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

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

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

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

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

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

Схема уставки для мультипоточного проекта Geant4
Схема моделируемой установки, P — источник частиц, протоны, 1 — мишень, 2 — детектор из вольфрама.

Схема моделируемого эксперимента: пучёк протонов 50 ГэВ направляется на мишень 1, и после прохождения мишени измеряем выделенную энергию в детекторе 2.

Миграция на многопоточность (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

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

ActionInitialization.hh

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

virtual void BuildForMaster() const;

А в файле

ActionInitialization.cpp

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

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

и

//Создание начального источника частиц и потока моделирования
void ActionInitialization::Build() const
{
    SetUserAction(new PrimaryGeneratorAction);//Задается начальный источник частиц
    // через обязательный класс 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

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

void ActionInitialization::Build() const

строкой кода

SetUserAction(new SteppingAction(eventAction));

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

#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);
    }
}

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

Реализация класса 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);
    // В конце каждого события
    // отображаем прогресс моделирования, отправляем номер события
    // event->GetEventID()
    fRunAction->DisplayProgress(event->GetEventID()+1);
}

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

fRunAction->FillEnergy(energ);

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

//Функция накопления выделенной энергии в детекторе
void EventAction::addParticle(const G4Step* step)
{   //Здесь мы берем выделенную энергию в обьеме за текущий шаг и накапливаем
    energ+=step->GetTotalEnergyDeposit();
}

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

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

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

#ifndef B2RunAction_h
#define B2RunAction_h 1

#include "G4UserRunAction.hh"
#include "globals.hh"// Global Constants and typedefs, глобальные константы и определения
#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;

а пропускать std::

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

//Конструктор класса 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);
    //Получение экземпляра менеджера общих переменных (accumulables).
    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/GeV/eventsNumber<< endl;
        //если мы делим на x/GeV, то переводим из внутрених
        //единиц в ГэВы, если умножаем, то ГэВы переводим в
        //внутренние единицы Geant4
        // И делим на количество всех запущеных частиц и получаем
        // среднюю выделенную энергию в мишени на одну начальную
        // частицу.
        //Закрываем поток вывода
        file_energy_dep.close();
    }
    else {//Если у нас вычислительный поток, не главный, можно что нибудь
        //вывысти
        G4cout
                << G4endl
                << "--------------------End of Local Run------------------------";
    }
}

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

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

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