Многопоточный пример на Geant4

Запуск и сборка примера на Geant4

Использование многопоточности в Geant4 позволяет существенно ускорить набор необходимой статистики. Использование персональных машин с 8-12 ядрами и высокой частотой процессора позволяет приблизится к скорости расчета, близкой к скорости расчета на одной ноде суперкомпьютера и забыть про необходимость набора необходимой статистики в течении недели, месяца.

Далее приведем пример, как преобразовать однопоточный пример Geant4, который приведен в статье полноценного рабочего примера по Geant4, в многопоточный.

В примере реализован простейший пример, где есть источник частиц, он направлен на мишень, за мишенью стоит детектор — чувствительный объем. Следует помнить, что c G4MTRunManager, пользовательские классы от G4VSensitiveDetector не работают.

Для регистрации частиц необходимо переопределять класс G4UserSteppingAction.

Миграция на многопоточность (multithreading)

Для начала надо изменить заголовочные файлы

Заменить в example3.cpp

#include "G4RunManager.hh"

На

// Подключаем заголовочные файлы
// RunManager - класс из ядра Geant4,
//должен быть включен обязательно для
//моделирования без распаралеливания
//а класс G4MTRunManager используется
//для запуска примера с распаралелливанием

// RunManager, класс из ядра Geant4,
//должен быть включен обязательно
#ifdef G4MULTITHREADED
#include "G4MTRunManager.hh"//Многопоточный
#else
#include "G4RunManager.hh"//Непараллельный, однопоточное моделирование
#endif

Здесь используется

#ifdef G4MULTITHREADED

который означает, что если Geant4 собран с поддержкой многопочности, то проект будет компилироваться с

#include "G4MTRunManager.hh"

иначе как обычно с

#include "G4RunManager.hh"

Далее в начале функции main ставим строчки

//Обеспечение песперебойной работы генератора случайных чисел в паралельном режиме
    CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine);
    CLHEP::HepRandom::setTheSeed(time(NULL));//Устанавливаем зерно для генератора
    //случайных чисел

что должно обеспечивать корректную работу генератору случайных чисел.

Далее меняем инициалицацию G4RunManager на G4MTRunManager

// Создание класс G4RunManager, он управляет течением программы и
// управляет событиями при запуске проекта
//Здесь использована директива ifdef, если при компилировании Geant4
//указывался флаг G4MULTITHREADED, то будет собран пример
//с многопоточным классом управления запуска моделирования G4MTRunManager
//если флаг не был указан, то Geant4 не поддерживает многопоточное
//моделирование и моделирование будет выполнятся в одном потоке с G4RunManager
#ifdef G4MULTITHREADED
  G4MTRunManager* runManager = new G4MTRunManager;
  //Выбор максимального числа потоков в выполняемой системе (на компьютере)
  // либо можно вручную задать число потоков
  runManager->SetNumberOfThreads( G4Threading::G4GetNumberOfCores());
#else
  G4RunManager* runManager = new G4RunManager;
#endif

Далее меняем файл

ExG4ActionInitialization01.hh

Добавляем функцию

virtual void BuildForMaster() const;

Далее в файле

ExG4ActionInitialization01.cpp

переписываем реализацию

//Создание главного потока
void ExG4ActionInitialization01::BuildForMaster() const {
    //Просто создаем и регистрируем поток
    RunAction* runAction = new RunAction;
    SetUserAction(runAction);
}

и

//Создание начального источника частиц и потока моделирования
void ExG4ActionInitialization01::Build() const
{
    SetUserAction(new ExG4PrimaryGeneratorAction01);//Задается начальный источник частиц
    // через обязательный класс ExG4PrimaryGeneratorAction01
    //Создание потока
    RunAction* runAction = new RunAction;
    //Регистрация потока
    SetUserAction(runAction);
    //В EventAction запускается 1 частица
    EventAction* eventAction = new EventAction(runAction);
    //Регистрируем
    SetUserAction(eventAction);
    //Так же создаем и регистрируем класс для получения происходящих событий
    //при моделировании
    SetUserAction(new SteppingAction(eventAction));
}}

Далее необходимо создать три класса

  • RunAction — класс который соответствует одному потоку
  • EventAction — класс который соответствует запуску одной частицы
  • SteppingAction — класс в котором вызывается функция virtual void UserSteppingAction(const G4Step*) при каждом событии.

Начнем в обратном порядке

Регистрация частиц в G4UserSteppingAction

Создадим пользовательский класс SteppingAction наследуемый от G4UserSteppingAction

Не забываем, что он регистрируется в функции

ExG4ActionInitialization01::Build()

строкой кода

SetUserAction(new SteppingAction(eventAction));

Заголовочный файл выглядит так

#ifndef SteppingAction_h
#define SteppingAction_h 1
#include "G4UserSteppingAction.hh"
#include "globals.hh"

class EventAction;
class G4LogicalVolume;

/// \brief The SteppingAction class Пользовательский класс
///  для регистрации событий
///
class SteppingAction : public G4UserSteppingAction
{
public:
    ///Конструктор
    SteppingAction(EventAction* eventAction);
    virtual ~SteppingAction();
    ///method from the base class
    ///Метод из базового класса, который вызывается,
    ///когда происходит какое либо событие с
    ///частицей
    virtual void UserSteppingAction(const G4Step*);
private:
    ///Указатель на текущий запуск начальной частицы
    EventAction* fEventAction;
};
#endif

Самым важным является реализация функции, которая вызывается при каждом событии

//Метод который вызывается когда происходит какое либо событие
void SteppingAction::UserSteppingAction(const G4Step* step){
    //Создаем указатель на трек для удобства
    G4Track * track = step->GetTrack();
    //Узнаем физический объем в котором находится частица
    G4VPhysicalVolume* vel=track->GetVolume();
    //Название физического объема в котором мы будем регистрировать
    //частицы, он у нас будет играть роль чуствительного объема или
    //детектора
    G4String name="Detector";
    //Проверка условия, что частица находится в интересуещем нас объеме
    if(name==vel->GetName()){
        //Если находится то передаем все что связано с шагом
        //в EventAction
        fEventAction->addParticle(step);
        //Нефизично, но нам нужно только подсчитать энергию частиц попавших
        //в детектор, поэтому эту частицу уничтожаем что бы она
        //не зарегистрировалась несколько раз
        step->GetTrack()->SetTrackStatus(fStopAndKill);
    }
}

Идея регистрации проста, проверяем названия объема в котором происходит событие, если это интересующий нас объем, тогда регистрируем то что нам нужно. В данном примере мы также уничтожаем частицу, что не является физичным (не является естественным), но для наших целей это хороший способ подсчитать энергию попадающих частиц.

Реализация класса G4UserEventAction

Для управления запуском одной начальной частицы используется класс

#ifndef EventAction_h
#define EventAction_h 1
#include "G4UserEventAction.hh"
#include "globals.hh"
#include "G4Step.hh"
#include "globals.hh"
class RunAction;

/// \brief The EventAction class Класс отвечающий
///  за запуск 1 первичной частицы
///
class EventAction: public G4UserEventAction
{
public:
    ///Конструктор
    EventAction(RunAction* runAction);
    ~EventAction();
    ///Функция, выполняемая перед запуском частицы
    virtual void BeginOfEventAction(const G4Event* event);
    ///Функция, выполняемая после запуска частицы
    virtual void EndOfEventAction(const G4Event* event);
    ///В эту переменную будем накапливать энергию частиц
    /// пришедших в детектор
    G4double energ;
    ///Функция регистрации частицы в детекторе
    void addParticle(const G4Step* step);
private:
    ///Указатель на поток выполнения
    RunAction* fRunAction;
};
#endif

Последовательно пройдемся по реализации отдельных функций

Конструктору класса требует адрес экземпляра класса RunAction из которого и вызывается EventAction

/// \brief EventAction::EventAction Реализация конструктора
/// \param runAction поток моделирования
/// инициализируем переменную потока моделирования
/// в конструкторе таким образом fRunAction(runAction)
///
EventAction::EventAction(RunAction *runAction): G4UserEventAction(),
    fRunAction(runAction) {
    //Обнуляем, но можем это сделать и позже
    energ=0;
}

Для простоты мы просто считаем энергию всех частиц попавших на детектор, энергию мы записываем в переменную energ.

Перед запуском начальной частицы вызывается функция

//Выполняется перед запуском начальной частицы
void EventAction::BeginOfEventAction(const G4Event*)
{     //Обнуляем для каждого запуска накопленную энергию частиц в детекторе
    energ=0;
}

В ней мы обнуляем накопленную энергию. После окончания моделирования выполняется функция

//Запускается после моделирования прохождения частицы
void EventAction::EndOfEventAction(const G4Event* event)
{   //Передаем полученное значение накопленой энергии частиц в
    //поток моделирования
    fRunAction->FillEnergy(energ);
    // в конце каждого события
    // отображаем прогресс моделирования
    fRunAction->DisplayProgress(event->GetEventID()+1);
}

Здесь мы передаем значение накопленной энергии в класс RunAction

fRunAction->FillEnergy(energ);

Реализовываем функцию накопления энергии частиц

//Функция накапления энергии частиц
void EventAction::addParticle(const G4Step* step)
{   //Здесь мы берем энергию (кинетическую, без массы покоя) в ГэВ-ах
    energ+=step->GetTrack()->GetKineticEnergy()/GeV;
    //если мы делим на x/GeV, то переводим из внутрених
    //единиц в ГэВы, если умножаем, то ГэВы переводим в
    //внутренние единицы Geant4
}

Она вызывается в классе SteppingAction

Реализация класса G4UserRunAction

Теперь рассмотрим класс реализации одного потока вычисления RunAction

#ifndef B2RunAction_h
#define B2RunAction_h 1
#include "G4UserRunAction.hh"
#include "globals.hh"
#include "G4Accumulable.hh"

class G4Run;
class ActionInitilization;

/// Run action class или класс RunAction, выполняется в отдельном потоке,
/// по сути это как один RunManager в однопоточном режиме.
/// Количество частиц, которое выполняется в одном  RunAction,
/// вычисляется динамически в процессе расчета, оно не фиксированно

class RunAction : public G4UserRunAction
{
public:
    RunAction();
    ~RunAction();
    ///В этой функции выводится прогресс моделирования, необязательная
    void DisplayProgress(G4int);
    ///Здесь делаем необходимые приготовления перед
    ///запуском счета
    void BeginOfRunAction(const G4Run*);
    ///После окончания расчета с вспомогательных потоков передаем
    ///полученные данные на главный поток
    ///В главном потоке данные объединяем
    ///и выводим в файл
    void EndOfRunAction(const G4Run*);
    ///Объявляем функцию, которую будем использовать для сбора
    ///данных моделирования
    void FillEnergy(G4double);
private:
    ///Класс G4Accumulable используется для передачи и объединения
    ///данных моделирования в главном потоке. По умолчанию работает
    ///только с простыми типами, double, int и т.д. однако функциональность
    ///можно расширить за счет пользовательского класса.
    G4Accumulable<G4double> Energy;
    ///Сюда будем записывать колическово смоделированных частиц
    G4int eventsNumber;
    ///Переменная необходимая для вывода прогресса моделирования
    G4int printModulo;
    ///Переменная для записи энергии частиц
    G4double particlesEnergy;
};
#endif

Здесь следует обратить внимание на переменную

G4Accumulable<G4double> Energy;

С помощью переменных G4Accumulable происходит передача данных между различным вычислительными потоками. Данные должны накапливаться и каждый вычислительный поток должен передать данные в G4Accumulable, а главный забрать из него данные после окончания моделирования.

По умолчанию G4Accumulable работает только с простыми типами, как double, G4double, int, float и так далее. И не работает с массивами. Проблема в том, что для других элементов и экземпляров классов нужно определять функции merge(), reset(), по умолчанию они определены только для простых типов. Однако можно сделать свой пользовательский класс и определить в нем эти функции merge(), reset(), позже будет статья с примером такого класса. Также можно использовать встроенные в Geant4 встроенные механизмы отображения, например как G4AnalysisManager.

Далее перейдем к объяснению реализации.

Для упрощения кода используем пространство имен

using namespace CLHEP;
using namespace std;

это для того, что бы например при выводе в файл не писать так

file_energy_dep std::<< particlesEnergy std::<< std::endl;

Далее приведем конструктор

//Конструктор класса RunAction, обратите внимание, что переменные G4Accumulable
//как Energy(0) должны быть объявлены с неким начальным значением
RunAction::RunAction()
    : G4UserRunAction(),Energy(0)
{
    //Создаем экземпляр класса управления (менеджера) переменных G4Accumulable
    G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
    //В нем регистрируем все переменные G4Accumulable для того что бы
    //они коректно обрабатывались
    accumulableManager->RegisterAccumulable(Energy);
}

Переменные G4Accumulable нужно писать в конструкторе с инициализируемым числом по умолчанию Energy(0). А в самом конструкторе зарегистрировать переменную, иначе данные не будут передаваться между потоками.

Далее определяем функцию, которая выполняется перед запуском вычислительного потока

//Функция выполняется перед началом запуска потока моделирования
void RunAction::BeginOfRunAction(const G4Run* aRun)
{ 
    //Берем полное число частиц, которое будет запущено
    eventsNumber = aRun->GetNumberOfEventToBeProcessed();
    //Переменая используемая для вывода прогресса
    printModulo = eventsNumber/100;
    if (printModulo == 0)
        printModulo = 1;
    //Обнуляем переменную, в нее будем записывать энергию
    //частиц приходящих в детектор
    particlesEnergy = 0;

    //Сообщаем G4RunManager-у сохранить зерно для генератора случайных чисел
    G4RunManager::GetRunManager()->SetRandomNumberStore(false);

    G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
    //Сбрасываем значения накапливаемых переменных accumulables на их начальные значения
    accumulableManager->Reset();
}

Обратите внимание, здесь мы тоже обнуляем переменную в которую записывается энергия частиц попавших в детектор

particlesEnergy = 0;

И сбрасываем значения переменных на указанные по умолчанию

accumulableManager->Reset();

Реализовываем функцию накопления энергии частиц

//Функция, в которой суммируем энергиию частиц попавших в детектор
void RunAction::FillEnergy(G4double energy)
{
    particlesEnergy += energy ;
}

Далее функция выполняемая в конце работы каждого потока моделирования и главного потока

//Функция выполняется в конце выполнения потока
//Следует помнить, что есть работающие потоки, которые считаю
//И есть главный поток Master, в котором моделирование не происходит
//но он завершается последним и в нем собираются все данные G4AccumulableManager-ом
void RunAction::EndOfRunAction(const G4Run* )
{
    //Проверяем, если не главный процесс, то копируем значение энергии
    //частиц в G4Accumulable
    if (!IsMaster()) {
        Energy=particlesEnergy;
    };
    //Передаем данные в главный поток, слияние данных
    //Данные должны быть передадены до слияния, иначе просто не попадут
    //в главный поток
    G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
    accumulableManager->Merge();
    if (IsMaster()) {//Выполняем в главном потоке
        G4cout
                << G4endl
                << "--------------------End of Global Run-----------------------";
        G4cout
                << G4endl;
        //Из G4Accumulable берем значение
        particlesEnergy=Energy.GetValue();

        //Объявляем поток вывода в файл
        std::ofstream file_energy_dep("total_particles_energy.txt");
        //Выводим в файл
        file_energy_dep<< particlesEnergy<< endl;
        //Закрываем поток вывода
        file_energy_dep.close();
    }
    else {//Если у нас вычислительный поток, не главный, можно что нибудь
        //вывысти
        G4cout
                << G4endl
                << "--------------------End of Local Run------------------------";
    }
}

Код содержит комментарии, однако обращаю внимание, что вывод данных делаем в главном потоке, так как только в нем собранны данные со всех потоков. Также запись в файл делать следует только в главном потоке, так как это относительно долгая операция, а если все потоки начнут писать в файлы (как например 60 потоков), то программа зависнет надолго.

Код многопоточного примера Geant4

Весь представленный код в примере находится здесь. Данный пример был разработан для учебных целей, если вам необходимо смоделировать достаточно сложный детектор, получить его отклик, или хотите провести интересный научный эксперимент и у вас есть для этого свободные ресурсы (деньги), тогда мы, возможно, сможем наладить плодотворное сотрудничество! Пожалуйста, обращаетесь тогда к администрации сайта.