diff --git a/dependencies.py b/dependencies.py new file mode 100755 index 00000000..6354b133 --- /dev/null +++ b/dependencies.py @@ -0,0 +1,1178 @@ +#!/usr/bin/env python3 + +from __future__ import annotations + +import curses +import json +import os +from dataclasses import dataclass, field +from typing import Callable, List, Optional, Tuple + + +# ============================ +# colors: edit these +# ============================ + +# foreground colors (use curses.COLOR_* or -1 for default) +COLOR_TITLE_FG = curses.COLOR_BLUE +COLOR_TEXT_FG = curses.COLOR_WHITE +COLOR_SELECTED_FG = curses.COLOR_WHITE +COLOR_SELECTED_BG = curses.COLOR_BLACK +COLOR_HINT_FG = curses.COLOR_YELLOW +COLOR_OK_FG = curses.COLOR_GREEN +COLOR_ERR_FG = curses.COLOR_RED +COLOR_KEY_FG = curses.COLOR_MAGENTA +COLOR_DIM_FG = curses.COLOR_CYAN + +# pair IDs (must be unique small ints) +PAIR_TITLE = 1 +PAIR_TEXT = 2 +PAIR_SELECTED = 3 +PAIR_HINT = 4 +PAIR_OK = 5 +PAIR_ERR = 6 +PAIR_KEY = 7 +PAIR_DIM = 8 + + +KOKKOS_BACKENDS = ["cpu", "cuda", "hip", "sycl"] +ADIOS2_MPI_MODES = ["non-mpi", "mpi"] + +MESSAGE: str = "" + + +@dataclass +class Settings: + cluster: str = "(custom)" + write_modulefiles: bool = False + overwrite: bool = False + install_prefix: str = os.path.join(os.path.expanduser("~"), ".entity") + + apps: dict = field( + default_factory=lambda: {"Kokkos": False, "adios2": False, "nt2py": False} + ) + + # versions + kokkos_version: str = "5.0.1" + adios2_version: str = "2.11.0" + + # options + kokkos_backend: str = "cpu" + kokkos_arch: str = "" + extra_kokkos_flags: List[str] = field(default_factory=list) + adios2_mpi: str = "non-mpi" + extra_adios2_flags: List[str] = field(default_factory=list) + + module_loads: List[str] = field(default_factory=list) + + def from_json(self, json_str: str) -> None: + data = json.loads(json_str) + self.cluster = data.get("cluster", self.cluster) + self.write_modulefiles = data.get("write_modulefiles", self.write_modulefiles) + self.overwrite = data.get("overwrite", self.overwrite) + self.install_prefix = data.get("install_prefix", self.install_prefix) + self.apps = data.get("dependencies", self.apps) + versions = data.get("versions", {}) + self.kokkos_version = versions.get("Kokkos", self.kokkos_version) + self.adios2_version = versions.get("adios2", self.adios2_version) + options = data.get("options", {}) + self.kokkos_backend = options.get("kokkos_backend", self.kokkos_backend) + self.kokkos_arch = options.get("kokkos_arch", self.kokkos_arch) + self.adios2_mpi = options.get("adios2_mpi", self.adios2_mpi) + self.module_loads = data.get("module_loads", self.module_loads) + + def apps_summary(self) -> str: + chosen = [k for k, v in self.apps.items() if v] + return ", ".join(chosen) if chosen else "(none)" + + def to_json(self) -> str: + return json.dumps( + { + "cluster": self.cluster, + "write_modulefiles": self.write_modulefiles, + "overwrite": self.overwrite, + "install_prefix": self.install_prefix, + "dependencies": self.apps, + "versions": { + "Kokkos": self.kokkos_version, + "adios2": self.adios2_version, + }, + "options": { + "kokkos_backend": self.kokkos_backend, + "kokkos_arch": self.kokkos_arch, + "adios2_mpi": self.adios2_mpi, + }, + "module_loads": self.module_loads, + }, + indent=2, + ) + + +def unindent(script: str) -> str: + script_lines = script.splitlines() + min_indent = min( + (len(line) - len(line.lstrip()) for line in script_lines if line.strip()), + default=0, + ) + trimmed_lines = [line[min_indent:] for line in script_lines] + if trimmed_lines[0] == "": + trimmed_lines = trimmed_lines[1:] + if trimmed_lines[-1] == "": + trimmed_lines = trimmed_lines[:-1] + return "\n".join(trimmed_lines) + + +def InstallKokkosScriptModfile(settings: Settings) -> tuple[str, str]: + if settings.apps.get("Kokkos", False): + prefix = settings.install_prefix + version = settings.kokkos_version + backend = settings.kokkos_backend + arch = settings.kokkos_arch.strip() + modules = "\n".join( + [f"module load {module} && \\" for module in settings.module_loads] + ) + src_path = f"{prefix}/src/kokkos" + install_path = ( + f"{prefix}/kokkos/{version}/{backend}{f'_{arch}' if arch else ''}" + ) + if os.path.exists(install_path) and not settings.overwrite: + raise FileExistsError( + f"Kokkos install path {install_path} already exists and overwrite is disabled" + ) + + extra_flags = "-D ".join(settings.extra_kokkos_flags) + cxx_standard = 20 if tuple(map(int, version.split("."))) >= (5, 0, 0) else 17 + + if arch == "": + arch = "NATIVE" + arch = arch.upper() + + script = f""" + # Kokkos installation + {modules} + rm -rf {src_path} && \\ + git clone https://github.com/kokkos/kokkos.git {src_path} && \\ + cd {src_path} && \\ + git checkout {version} && \\ + cmake -B build \\ + -D CMAKE_CXX_STANDARD={cxx_standard} \\ + -D CMAKE_CXX_EXTENSIONS=OFF \\ + -D CMAKE_POSITION_INDEPENDENT_CODE=TRUE \\ + -D Kokkos_ARCH_{arch}=ON {f'-D Kokkos_ENABLE_{backend.upper()}=ON' if backend != 'cpu' else ''} \\ + -D CMAKE_INSTALL_PREFIX={install_path} {extra_flags} && \\ + cmake --build build -j $(nproc) && \\ + cmake --install build + """ + + modfile = f""" + #%Module1.0###################################################################### + ## + ## Kokkos @ {backend} @ {arch} modulefile + ## + ################################################################################# + proc ModulesHelp {{ }} {{ + puts stderr \"\\tKokkos @ {backend} @ {arch}\\n\" + }} + + module-whatis \"Sets up Kokkos @ {backend} @ {arch}\" + + conflict kokkos + {modules} + + set basedir {install_path} + prepend-path PATH $basedir/bin + setenv Kokkos_DIR $basedir + + setenv Kokkos_ARCH_{arch} ON + {f'setenv Kokkos_ENABLE_{backend.upper()} ON' if backend != 'cpu' else ''} + """ + + return (unindent(script), unindent(modfile)) + + else: + return ("""# skipping Kokkos install""", "") + + +def InstallAdios2Script(settings: Settings) -> tuple[str, str]: + if settings.apps.get("adios2", False): + prefix = settings.install_prefix + version = settings.adios2_version + mpi_mode = settings.adios2_mpi + modules = "\n".join( + [f"module load {module} && \\" for module in settings.module_loads] + ) + src_path = f"{prefix}/src/adios2" + install_path = f"{prefix}/adios2/{version}/{mpi_mode}" + if os.path.exists(install_path) and not settings.overwrite: + raise FileExistsError( + f"Adios2 install path {install_path} already exists and overwrite is disabled" + ) + + extra_flags = "-D ".join(settings.extra_adios2_flags) + cxx_standard = ( + 20 + if tuple(map(int, settings.kokkos_version.split("."))) >= (5, 0, 0) + else 17 + ) + + with_mpi = "ON" if mpi_mode == "mpi" else "OFF" + + script = f""" + # Adios2 installation + {modules} + rm -rf {src_path} && \\ + git clone https://github.com/ornladios/ADIOS2.git {src_path} && \\ + cd {src_path} && \\ + git checkout v{version} && \\ + cmake -B build \\ + -D CMAKE_CXX_STANDARD={cxx_standard} \\ + -D CMAKE_CXX_EXTENSIONS=OFF \\ + -D CMAKE_POSITION_INDEPENDENT_CODE=TRUE \\ + -D BUILD_SHARED_LIBS=ON \\ + -D ADIOS2_USE_Python=OFF \\ + -D ADIOS2_USE_Fortran=OFF \\ + -D ADIOS2_USE_ZeroMQ=OFF \\ + -D BUILD_TESTING=OFF \\ + -D ADIOS2_BUILD_EXAMPLES=OFF \\ + -D ADIOS2_USE_HDF5=OFF \\ + -D ADIOS2_USE_MPI={with_mpi} \\ + -D CMAKE_INSTALL_PREFIX={install_path} {extra_flags} && \\ + cmake --build build -j $(nproc) && \\ + cmake --install build + """ + + modfile = f""" + #%Module1.0###################################################################### + ## + ## ADIOS2 @ {mpi_mode} modulefile + ## + ################################################################################# + proc ModulesHelp {{ }} {{ + puts stderr \"\\tADIOS2 @ {mpi_mode}\\n\" + }} + + module-whatis \"Sets up ADIOS2 @ {mpi_mode}\" + + conflict adios2 + {modules} + + set basedir {install_path} + prepend-path PATH $basedir/bin + setenv ADIOS2_DIR $basedir + + setenv ADIOS2_USE_MPI {with_mpi} + """ + + return (unindent(script), unindent(modfile)) + + else: + return ("""# skipping Adios2 install""", "") + + +def InstallNt2pyScript(settings: Settings) -> str: + if settings.apps.get("nt2py", False): + prefix = settings.install_prefix + modules = "\n".join( + [f"module load {module} && \\" for module in settings.module_loads] + ) + install_path = f"{prefix}/.venv" + + script = f""" + # nt2py installation + {modules} + rm -rf {install_path} && \\ + python3 -m venv {install_path} && \\ + source {install_path}/bin/activate && \\ + pip install nt2py && \\ + deactivate + """ + return unindent(script) + else: + return """# skipping nt2py install""" + + +PRESETS = { + "rusty": {"module_loads": []}, + "stellar": {"module_loads": []}, + "perlmutter": {"module_loads": []}, + "frontier": {"module_loads": []}, + "aurora": {"module_loads": []}, +} + + +def apply_preset(s: Settings, name: str) -> None: + s.cluster = name + s.install_prefix = os.path.join(os.path.expanduser("~"), ".entity") + cluster_preset = PRESETS.get(name, {}) + s.apps["Kokkos"] = True + s.apps["adios2"] = True + s.apps["nt2py"] = True + s.module_loads = cluster_preset.get("module_loads", []) + s.extra_kokkos_flags = cluster_preset.get("extra_kokkos_flags", []) + s.extra_adios2_flags = cluster_preset.get("extra_adios2_flags", []) + + +def on_install_confirmed(settings: Settings) -> None: + global MESSAGE + os.makedirs(settings.install_prefix, exist_ok=True) + kokkos_script, kokkos_modfile = InstallKokkosScriptModfile(settings) + adios2_script, adios2_modfile = InstallAdios2Script(settings) + with open(os.path.join(settings.install_prefix, "install.sh"), "w") as f: + f.write("#!/usr/bin/env bash\n\n") + f.write(kokkos_script) + f.write("\n\n") + f.write(adios2_script) + f.write("\n") + if settings.write_modulefiles: + os.makedirs(os.path.join(settings.install_prefix, "modules"), exist_ok=True) + if kokkos_modfile != "": + kokkos_modfile_file = os.path.join( + settings.install_prefix, + "modules", + "kokkos", + settings.kokkos_backend + + ( + f"_{settings.kokkos_arch.strip()}" + if settings.kokkos_arch.strip() + else "" + ), + settings.kokkos_version, + ) + os.makedirs(os.path.dirname(kokkos_modfile_file), exist_ok=True) + if os.path.exists(kokkos_modfile_file) and not settings.overwrite: + raise FileExistsError( + f"modulefile {kokkos_modfile_file} already exists and overwrite is disabled" + ) + with open(kokkos_modfile_file, "w") as f: + f.write(kokkos_modfile) + if adios2_modfile != "": + adios2_modfile_file = os.path.join( + settings.install_prefix, + "modules", + "adios2", + settings.adios2_mpi, + settings.adios2_version, + ) + os.makedirs(os.path.dirname(adios2_modfile_file), exist_ok=True) + if os.path.exists(adios2_modfile_file) and not settings.overwrite: + raise FileExistsError( + f"modulefile {adios2_modfile_file} already exists and overwrite is disabled" + ) + with open(adios2_modfile_file, "w") as f: + f.write(adios2_modfile) + + os.chmod(os.path.join(settings.install_prefix, "install.sh"), 0o755) + MESSAGE = f"- installation script written to {os.path.join(settings.install_prefix, 'install.sh')}!\n" + MESSAGE += " please read and verify it before running.\n\n" + if settings.write_modulefiles: + MESSAGE += f"- module files have been written to {os.path.join(settings.install_prefix, 'modules')} directory.\n" + MESSAGE += f" add them to your .rc script as `module use --append {os.path.join(settings.install_prefix, 'modules')}`\n\n" + + if settings.apps.get("nt2py", False): + MESSAGE += ( + "- nt2py installed in a new virtual environment at " + f"{os.path.join(settings.install_prefix, '.venv')}.\n" + ) + MESSAGE += " activate it with `source {}/bin/activate`.\n\n".format( + os.path.join(settings.install_prefix, ".venv") + ) + + settings_json = os.path.join(settings.install_prefix, "settings.json") + with open(settings_json, "w") as f: + f.write(settings.to_json()) + return + + +@dataclass +class MenuItem: + label: str + hint: str = "" + right: Optional[Callable[[], str]] = None + on_enter: Optional[Callable[[], None]] = None + on_space: Optional[Callable[[], None]] = None + disabled: Optional[Callable[[], bool]] = None + + +class TuiExitInstall(Exception): + pass + + +class App: + def __init__(self, stdscr): + self.stdscr = stdscr + self.s = Settings() + if os.path.exists(os.path.join(self.s.install_prefix, "settings.json")): + with open(os.path.join(self.s.install_prefix, "settings.json"), "r") as f: + data = json.load(f) + self.s.from_json(json.dumps(data)) + + self.state = "mainmenu" + self.stack: List[Tuple[str, int]] = [] + self.selected = 0 + self.scroll = 0 + self.message = "use arrows or j/k" + + self.mod_sel = 0 + self.mod_scroll = 0 + + self._init_curses() + + def _init_curses(self) -> None: + curses.curs_set(0) + self.stdscr.keypad(True) + curses.noecho() + curses.cbreak() + + if curses.has_colors(): + curses.start_color() + curses.use_default_colors() + curses.init_pair(PAIR_TITLE, COLOR_TITLE_FG, -1) + curses.init_pair(PAIR_TEXT, COLOR_TEXT_FG, -1) + curses.init_pair(PAIR_SELECTED, COLOR_SELECTED_FG, COLOR_SELECTED_BG) + curses.init_pair(PAIR_HINT, COLOR_HINT_FG, -1) + curses.init_pair(PAIR_OK, COLOR_OK_FG, -1) + curses.init_pair(PAIR_ERR, COLOR_ERR_FG, -1) + curses.init_pair(PAIR_KEY, COLOR_KEY_FG, -1) + curses.init_pair(PAIR_DIM, COLOR_DIM_FG, -1) + + def cp(self, pair_id: int) -> int: + return curses.color_pair(pair_id) if curses.has_colors() else 0 + + # ----- formatting helpers ----- + + def checkbox(self, on: bool) -> str: + return "[x]" if on else "[ ]" + + def pill(self, on: bool) -> str: + return "[on]" if on else "[off]" + + def kokkos_right(self) -> str: + arch = self.s.kokkos_arch.strip() or "-" + return f"{self.s.kokkos_version} · {self.s.kokkos_backend} · {arch}" + + def adios2_right(self) -> str: + return f"{self.s.adios2_version} · {self.s.adios2_mpi}" + + # ----- nav stack ----- + + def push(self, st: str) -> None: + self.stack.append((self.state, self.selected)) + self.state = st + self.selected = 0 + self.scroll = 0 + self.message = "" + + def pop(self) -> None: + if self.stack: + self.state, self.selected = self.stack.pop() + else: + self.state, self.selected = "mainmenu", 0 + self.scroll = 0 + self.message = "" + + # ----- drawing ----- + + def add(self, y: int, x: int, s: str, attr: int = 0) -> None: + try: + self.stdscr.addstr(y, x, s, attr) + except curses.error: + pass + + def hline(self, y: int) -> None: + _, w = self.stdscr.getmaxyx() + try: + self.stdscr.hline(y, 0, curses.ACS_HLINE, max(0, w - 1)) + except curses.error: + pass + + def draw_keybar(self, y: int, x: int, pairs: List[Tuple[str, str]]) -> None: + cur_x = x + for key, action in pairs: + self.add(y, cur_x, key, self.cp(PAIR_KEY) | curses.A_BOLD) + cur_x += len(key) + self.add(y, cur_x, " ", self.cp(PAIR_DIM)) + cur_x += 1 + self.add(y, cur_x, action, self.cp(PAIR_HINT)) + cur_x += len(action) + self.add(y, cur_x, " ", self.cp(PAIR_DIM)) + cur_x += 3 + + def breadcrumb(self) -> str: + if self.state == "mainmenu": + return "mainmenu" + if self.state == "custom": + return "mainmenu › custom install" + if self.state == "dependencies": + return "mainmenu › custom install › dependencies" + if self.state == "versions": + return "mainmenu › custom install › versions" + if self.state == "options": + return "mainmenu › custom install › options" + if self.state == "cluster": + return "mainmenu › cluster-specific" + if self.state == "preset_applied": + return f"mainmenu › cluster-specific › {self.s.cluster}" + return "mainmenu" + + def draw_menu(self, title: str, prompt: str, items: List[MenuItem]) -> None: + self.stdscr.erase() + h, w = self.stdscr.getmaxyx() + + self.add(0, 2, title, self.cp(PAIR_TITLE) | curses.A_BOLD) + bc = self.breadcrumb() + self.add(0, max(2, w - 2 - len(bc)), bc, self.cp(PAIR_DIM)) + + self.draw_keybar( + 1, + 2, + [ + ("↑/↓/j/k", "move"), + ("enter", "select"), + ("space", "toggle/cycle"), + ("b", "back"), + ("q", "quit"), + ], + ) + self.hline(2) + + status1 = f"cluster: {self.s.cluster} write modulefiles: {self.pill(self.s.write_modulefiles)} module loads: {len(self.s.module_loads)}" + status2 = ( + f"prefix: {self.s.install_prefix} dependencies: {self.s.apps_summary()}" + ) + self.add(3, 2, status1[: w - 4], self.cp(PAIR_TEXT)) + self.add(4, 2, status2[: w - 4], self.cp(PAIR_TEXT)) + self.hline(5) + + self.add(6, 2, prompt[: w - 4], self.cp(PAIR_TEXT) | curses.A_BOLD) + + list_y = 8 + footer_h = 3 + view_h = max(1, h - list_y - footer_h) + n = len(items) + + if n == 0: + self.add(list_y, 2, "(empty)", self.cp(PAIR_HINT)) + else: + self.selected = max(0, min(self.selected, n - 1)) + + if self.selected < self.scroll: + self.scroll = self.selected + if self.selected >= self.scroll + view_h: + self.scroll = self.selected - view_h + 1 + self.scroll = max(0, min(self.scroll, max(0, n - view_h))) + + shown = items[self.scroll : self.scroll + view_h] + + for i, it in enumerate(shown): + idx = self.scroll + i + sel = idx == self.selected + dis = bool(it.disabled and it.disabled()) + + row_attr = ( + self.cp(PAIR_SELECTED) | curses.A_BOLD + if sel + else (self.cp(PAIR_DIM) if dis else self.cp(PAIR_TEXT)) + ) + self.add(list_y + i, 2, f" {it.label}"[: w - 4], row_attr) + + if it.right: + rt = (it.right() or "").strip() + if rt: + rt = rt[: max(0, w - 6)] + x = max(2, w - 2 - len(rt)) + rt_attr = ( + row_attr + if sel + else (self.cp(PAIR_HINT) if not dis else self.cp(PAIR_DIM)) + ) + self.add(list_y + i, x, rt, rt_attr) + + if sel and it.hint: + self.add( + list_y + i, + min(w - 4, 30), + f" {it.hint}"[: w - 4], + self.cp(PAIR_HINT), + ) + + self.hline(h - 3) + msg = self.message or "" + if msg: + is_err = msg.startswith("error") + attr = (self.cp(PAIR_ERR) if is_err else self.cp(PAIR_OK)) | curses.A_BOLD + self.add(h - 2, 2, msg[: w - 4], attr) + self.stdscr.refresh() + + # ----- modals ----- + + def input_box(self, title: str, prompt: str, initial: str) -> Optional[str]: + h, w = self.stdscr.getmaxyx() + win_h, win_w = 9, min(86, max(46, w - 6)) + top, left = max(0, (h - win_h) // 2), max(0, (w - win_w) // 2) + + win = curses.newwin(win_h, win_w, top, left) + win.keypad(True) + win.border() + + win.addstr(1, 2, title[: win_w - 4], self.cp(PAIR_TITLE) | curses.A_BOLD) + win.addstr(2, 2, prompt[: win_w - 4], self.cp(PAIR_TEXT)) + + buf = list(initial) + curses.curs_set(1) + + while True: + win.addstr(4, 2, " " * (win_w - 4), self.cp(PAIR_TEXT)) + text = "".join(buf) + if len(text) > win_w - 4: + text = text[-(win_w - 4) :] + win.addstr(4, 2, text, self.cp(PAIR_TEXT) | curses.A_BOLD) + win.addstr(6, 2, "enter=ok esc=cancel", self.cp(PAIR_DIM)) + win.refresh() + + ch = win.getch() + if ch == 27: + curses.curs_set(0) + return None + if ch in (curses.KEY_ENTER, 10, 13): + curses.curs_set(0) + return "".join(buf).strip() + if ch in (curses.KEY_BACKSPACE, 127, 8): + if buf: + buf.pop() + elif 32 <= ch <= 126: + buf.append(chr(ch)) + + def confirm_install(self) -> bool: + arch = self.s.kokkos_arch.strip() or "-" + lines = [ + f"cluster: {self.s.cluster}", + f"overwrite existing files: {self.pill(self.s.overwrite)}", + f"write modulefiles: {self.pill(self.s.write_modulefiles)}", + f"module loads: {len(self.s.module_loads)}", + f"prefix: {self.s.install_prefix}", + f"dependencies: {self.s.apps_summary()}", + f"kokkos: {self.s.kokkos_version} · {self.s.kokkos_backend} · {arch}", + f"adios2: {self.s.adios2_version} · {self.s.adios2_mpi}", + "", + "confirm install?", + ] + + h, w = self.stdscr.getmaxyx() + win_h, win_w = min(16, max(10, h - 6)), min(94, max(52, w - 6)) + top, left = max(0, (h - win_h) // 2), max(0, (w - win_w) // 2) + + win = curses.newwin(win_h, win_w, top, left) + win.keypad(True) + win.border() + win.addstr(1, 2, "confirm", self.cp(PAIR_TITLE) | curses.A_BOLD) + + y = 3 + for ln in lines[: win_h - 6]: + win.addstr(y, 2, ln[: win_w - 4], self.cp(PAIR_TEXT)) + y += 1 + + win.addstr(win_h - 3, 2, "y=yes n=no", self.cp(PAIR_DIM)) + win.refresh() + + while True: + ch = win.getch() + if ch in (ord("y"), ord("Y")): + return True + if ch in (ord("n"), ord("N"), 27): + return False + + # ----- helpers ----- + + def cycle(self, current: str, options: List[str]) -> str: + if current not in options: + return options[0] + i = options.index(current) + return options[(i + 1) % len(options)] + + # ----- module editor ----- + + def module_editor(self) -> None: + while True: + self.stdscr.erase() + h, w = self.stdscr.getmaxyx() + + self.add(0, 2, "module lines", self.cp(PAIR_TITLE) | curses.A_BOLD) + self.draw_keybar( + 1, + 2, + [ + ("↑/↓/j/k", "move"), + ("enter", "edit"), + ("a", "add"), + ("d", "delete"), + ("u/m", "reorder"), + ("b", "back"), + ], + ) + self.hline(2) + + self.add(3, 2, f"lines: {len(self.s.module_loads)}", self.cp(PAIR_TEXT)) + self.hline(4) + + lines = self.s.module_loads + n = len(lines) + list_y = 6 + view_h = max(1, h - list_y - 3) + + if n == 0: + self.add(list_y, 2, "(empty) press a to add", self.cp(PAIR_HINT)) + else: + self.mod_sel = max(0, min(self.mod_sel, n - 1)) + if self.mod_sel < self.mod_scroll: + self.mod_scroll = self.mod_sel + if self.mod_sel >= self.mod_scroll + view_h: + self.mod_scroll = self.mod_sel - view_h + 1 + self.mod_scroll = max(0, min(self.mod_scroll, max(0, n - view_h))) + + shown = lines[self.mod_scroll : self.mod_scroll + view_h] + for i, ln in enumerate(shown): + idx = self.mod_scroll + i + sel = idx == self.mod_sel + attr = ( + self.cp(PAIR_SELECTED) | curses.A_BOLD + if sel + else self.cp(PAIR_TEXT) + ) + self.add(list_y + i, 2, f" {ln}"[: w - 4], attr) + + self.hline(h - 3) + self.add( + h - 2, + 2, + "tip: example: cuda/12.9"[: w - 4], + self.cp(PAIR_HINT), + ) + self.stdscr.refresh() + + ch = self.stdscr.getch() + if ch in (ord("q"), ord("Q"), ord("b"), 8, 127): + return + + n = len(self.s.module_loads) + self.mod_sel = 0 if n == 0 else max(0, min(self.mod_sel, n - 1)) + + if ch in (curses.KEY_UP, ord("k"), ord("K")) and n: + self.mod_sel = (self.mod_sel - 1) % n + continue + if ch in (curses.KEY_DOWN, ord("j"), ord("J")) and n: + self.mod_sel = (self.mod_sel + 1) % n + continue + + if ch in (ord("a"), ord("A")): + val = self.input_box("add module line", "example: cuda/12.9", "") + if val: + self.s.module_loads.append(val) + self.mod_sel = len(self.s.module_loads) - 1 + continue + + if ch in (ord("d"), ord("D")): + if n == 0: + continue + val = self.input_box("delete line", "type 'delete' to confirm:", "") + if val == "delete": + del self.s.module_loads[self.mod_sel] + self.mod_sel = max( + 0, min(self.mod_sel, len(self.s.module_loads) - 1) + ) + continue + + if ch in (ord("u"), ord("U")): + if n >= 2 and self.mod_sel > 0: + i = self.mod_sel + self.s.module_loads[i - 1], self.s.module_loads[i] = ( + self.s.module_loads[i], + self.s.module_loads[i - 1], + ) + self.mod_sel -= 1 + continue + + if ch in (ord("m"), ord("M")): + if n >= 2 and self.mod_sel < n - 1: + i = self.mod_sel + self.s.module_loads[i + 1], self.s.module_loads[i] = ( + self.s.module_loads[i], + self.s.module_loads[i + 1], + ) + self.mod_sel += 1 + continue + + if ch in (curses.KEY_ENTER, 10, 13): + if n == 0: + continue + cur = self.s.module_loads[self.mod_sel] + val = self.input_box("edit module line", "edit the selected line:", cur) + if val is not None and val.strip(): + self.s.module_loads[self.mod_sel] = val.strip() + continue + + # ----- menus ----- + + def versions_menu(self) -> Tuple[str, str, List[MenuItem]]: + def edit_kokkos(): + val = self.input_box( + "kokkos version", "enter version/tag:", self.s.kokkos_version + ) + if val: + self.s.kokkos_version = val.strip() + + def edit_adios2(): + val = self.input_box( + "adios2 version", "enter version/tag:", self.s.adios2_version + ) + if val: + self.s.adios2_version = val.strip() + + return ( + "versions", + "set versions:", + [ + MenuItem( + "kokkos version", + "enter to edit", + right=lambda: self.s.kokkos_version, + on_enter=edit_kokkos, + ), + MenuItem( + "adios2 version", + "enter to edit", + right=lambda: self.s.adios2_version, + on_enter=edit_adios2, + ), + MenuItem("back", "return", on_enter=self.pop), + ], + ) + + def options_menu(self) -> Tuple[str, str, List[MenuItem]]: + def cycle_kokkos(): + self.s.kokkos_backend = self.cycle(self.s.kokkos_backend, KOKKOS_BACKENDS) + + def edit_kokkos_arch(): + val = self.input_box( + "kokkos arch", "enter arch text (free-form):", self.s.kokkos_arch + ) + if val is not None: + self.s.kokkos_arch = val.strip() + + def cycle_adios2(): + self.s.adios2_mpi = self.cycle(self.s.adios2_mpi, ADIOS2_MPI_MODES) + + return ( + "options", + "set build options:", + [ + MenuItem( + "kokkos backend", + "space cycles: cpu/cuda/hip/sycl", + right=lambda: self.s.kokkos_backend, + on_enter=cycle_kokkos, + on_space=cycle_kokkos, + disabled=lambda: not self.s.apps.get("Kokkos", False), + ), + MenuItem( + "kokkos arch", + "enter to edit (optional)", + right=lambda: (self.s.kokkos_arch.strip() or "-"), + on_enter=edit_kokkos_arch, + disabled=lambda: not self.s.apps.get("Kokkos", False), + ), + MenuItem( + "adios2 mpi", + "space cycles: non-mpi/mpi", + right=lambda: self.s.adios2_mpi, + on_enter=cycle_adios2, + on_space=cycle_adios2, + disabled=lambda: not self.s.apps.get("adios2", False), + ), + MenuItem("back", "return", on_enter=self.pop), + ], + ) + + def menu_main(self) -> Tuple[str, str, List[MenuItem]]: + return ( + "entity deps", + "main menu:", + [ + MenuItem( + "custom install", + "edit settings then install", + on_enter=lambda: self.push("custom"), + ), + MenuItem( + "cluster-specific", + "apply a cluster-specific preset (editable)", + on_enter=lambda: self.push("cluster"), + ), + MenuItem("exit", "", on_enter=lambda: setattr(self, "state", "exit")), + ], + ) + + def menu_custom(self) -> Tuple[str, str, List[MenuItem]]: + def toggle_write_modulefiles(): + self.s.write_modulefiles = not self.s.write_modulefiles + + def toggle_overwrite(): + self.s.overwrite = not self.s.overwrite + + def edit_prefix(): + val = self.input_box( + "install location", "enter install prefix:", self.s.install_prefix + ) + if val: + self.s.install_prefix = os.path.expanduser(val.strip()) + + def go_apps(): + self.push("dependencies") + + def go_versions(): + self.push("versions") + + def go_options(): + self.push("options") + + def do_install(): + if not self.confirm_install(): + self.message = "cancelled." + return + on_install_confirmed(self.s) + raise TuiExitInstall + + return ( + "custom install", + "settings:", + [ + MenuItem( + "overwrite existing files", + "whether to overwrite existing files", + right=lambda: "enabled" if self.s.overwrite else "disabled", + on_enter=toggle_overwrite, + on_space=toggle_overwrite, + ), + MenuItem( + "write modulefiles", + "whether to create module files", + right=lambda: "enabled" if self.s.write_modulefiles else "disabled", + on_enter=toggle_write_modulefiles, + on_space=toggle_write_modulefiles, + ), + MenuItem( + "module load lines", + "add/remove modules to load", + right=lambda: f"{len(self.s.module_loads)} entry(s)", + on_enter=self.module_editor, + ), + MenuItem( + "install location", + "root location where modules and dependencies are installed", + right=lambda: self.s.install_prefix, + on_enter=edit_prefix, + ), + MenuItem( + "dependencies to install", + "select which dependencies to install", + right=lambda: self.s.apps_summary(), + on_enter=go_apps, + ), + MenuItem( + "versions", + "edit dependency versions", + right=lambda: " · ".join( + [ + a + for (a, ae) in zip( + [ + f"kokkos {self.s.kokkos_version}", + f"adios2 {self.s.adios2_version}", + ], + [ + self.s.apps.get(app, False) + for app in ["Kokkos", "adios2"] + ], + ) + if ae + ] + ), + on_enter=go_versions, + ), + MenuItem( + "options", + "pick backends/architectures/mpi", + right=lambda: " · ".join( + [ + a + for (a, ae) in zip( + [ + f"kokkos {self.s.kokkos_backend}/{self.s.kokkos_arch.strip() or '-'}", + f"adios2 {self.s.adios2_mpi}", + ], + [ + self.s.apps.get(app, False) + for app in ["Kokkos", "adios2"] + ], + ) + if ae + ] + ), + on_enter=go_options, + ), + MenuItem("install", "", on_enter=do_install), + MenuItem("back", "", on_enter=self.pop), + ], + ) + + def menu_apps(self) -> Tuple[str, str, List[MenuItem]]: + def toggle(k: str): + self.s.apps[k] = not self.s.apps.get(k, False) + + return ( + "dependencies", + "select the dependencies:", + [ + MenuItem( + f"{self.checkbox(self.s.apps.get('Kokkos', False))} kokkos", + "", + on_enter=lambda: toggle("Kokkos"), + on_space=lambda: toggle("Kokkos"), + right=self.kokkos_right, + ), + MenuItem( + f"{self.checkbox(self.s.apps.get('adios2', False))} adios2", + "", + on_enter=lambda: toggle("adios2"), + on_space=lambda: toggle("adios2"), + right=self.adios2_right, + ), + MenuItem( + f"{self.checkbox(self.s.apps.get('nt2py', False))} nt2py", + "", + on_enter=lambda: toggle("nt2py"), + on_space=lambda: toggle("nt2py"), + ), + MenuItem("back", "", on_enter=self.pop), + ], + ) + + def menu_cluster(self) -> Tuple[str, str, List[MenuItem]]: + def choose(name: str): + apply_preset(self.s, name) + self.push("custom") + + return ( + "cluster-specific", + "pick a preset:", + [ + MenuItem("rusty", "apply preset", on_enter=lambda: choose("rusty")), + MenuItem("stellar", "apply preset", on_enter=lambda: choose("stellar")), + MenuItem( + "perlmutter", "apply preset", on_enter=lambda: choose("perlmutter") + ), + MenuItem( + "frontier", "apply preset", on_enter=lambda: choose("frontier") + ), + MenuItem("aurora", "apply preset", on_enter=lambda: choose("aurora")), + MenuItem("back", "", on_enter=self.pop), + ], + ) + + def get_menu(self) -> Tuple[str, str, List[MenuItem]]: + if self.state == "mainmenu": + return self.menu_main() + if self.state == "custom": + return self.menu_custom() + if self.state == "dependencies": + return self.menu_apps() + if self.state == "versions": + return self.versions_menu() + if self.state == "options": + return self.options_menu() + if self.state == "cluster": + return self.menu_cluster() + self.state = "mainmenu" + return self.menu_main() + + # ----- navigation ----- + + def is_disabled(self, it: MenuItem) -> bool: + return bool(it.disabled and it.disabled()) + + def move_sel(self, items: List[MenuItem], delta: int) -> None: + if not items: + return + n = len(items) + start = self.selected + for _ in range(n): + self.selected = (self.selected + delta) % n + if not self.is_disabled(items[self.selected]): + return + self.selected = start + + def activate(self, items: List[MenuItem], enter: bool) -> None: + if not items: + return + it = items[self.selected] + if self.is_disabled(it): + self.message = "error: option disabled." + return + fn = it.on_enter if enter else it.on_space + if fn: + fn() + + # ----- loop ----- + + def run(self) -> None: + while True: + if self.state == "exit": + return + + title, prompt, items = self.get_menu() + self.draw_menu(title, prompt, items) + + ch = self.stdscr.getch() + + if ch in (ord("q"), ord("Q")): + self.state = "exit" + continue + + if ch in (ord("b"), 8, 127): + self.pop() + continue + + if ch in (curses.KEY_UP, ord("k"), ord("K")): + self.move_sel(items, -1) + continue + + if ch in (curses.KEY_DOWN, ord("j"), ord("J")): + self.move_sel(items, +1) + continue + + if ch in (curses.KEY_ENTER, 10, 13): + self.activate(items, enter=True) + continue + + if ch == ord(" "): + self.activate(items, enter=False) + continue + + +def _wrapper_capture(stdscr) -> None: + app = App(stdscr) + try: + app.run() + except TuiExitInstall: + on_install_confirmed(app.s) + raise + + +if __name__ == "__main__": + try: + curses.wrapper(_wrapper_capture) + raise SystemExit(0) + except TuiExitInstall: + print(MESSAGE) + raise SystemExit(0) + except KeyboardInterrupt: + raise SystemExit(130) diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix index d8ae115c..fc51647b 100644 --- a/dev/nix/kokkos.nix +++ b/dev/nix/kokkos.nix @@ -17,15 +17,16 @@ let rocprim rocminfo rocm-smi + pkgs.llvmPackages_19.clang-tools ]; "CUDA" = with pkgs.cudaPackages; [ - llvmPackages_18.clang-tools + llvmPackages_19.clang-tools cudatoolkit cuda_cudart pkgs.gcc13 ]; "NONE" = [ - pkgs.llvmPackages_18.clang-tools + pkgs.llvmPackages_19.clang-tools pkgs.gcc13 ]; }; diff --git a/input.example.toml b/input.example.toml index 9a4aedee..7d0d73f1 100644 --- a/input.example.toml +++ b/input.example.toml @@ -195,6 +195,59 @@ # @from: `scales.larmor0` # @value: `1 / larmor0` +[radiation] + [radiation.drag] + [radiation.drag.synchrotron] + # Radiation reaction limit gamma-factor for synchrotron + # @type: float [> 0.0] + # @default: 1.0 + # @note: [required] if one of the species has `radiative_drag = "synchrotron"` + gamma_rad = "" + + [radiation.drag.compton] + # Radiation reaction limit gamma-factor for Compton drag + # @type: float [> 0.0] + # @default: 1.0 + # @note: [required] if one of the species has `radiative_drag = "compton"` + gamma_rad = "" + + [radiation.emission] + [radiation.emission.synchrotron] + # Gamma-factor of a particle emitting synchrotron photons at energy `m0 c^2` in fiducial magnetic field `B0` + # @type: float [> 1.0] + # @default: 10.0 + gamma_qed = "" + # Minimum photon energy for synchrotron emission (units of `m0 c^2`) + # @type: float [> 0.0] + # @default: 1e-4 + photon_energy_min = "" + # Weights for the emitted synchrotron photons + # @type: float [> 0.0] + # @default: 1.0 + photon_weight = "" + # Index of species for the emitted photon + # @type: ushort [> 0] + # @required + photon_species = "" + + [radiation.emission.compton] + # Gamma-factor of a particle emitting inverse Compton photons at energy `m0 c^2` in fiducial magnetic field `B0` + # @type: float [> 1.0] + # @default: 10.0 + gamma_qed = "" + # Minimum photon energy for inverse Compton emission (units of `m0 c^2`) + # @type: float [> 0.0] + # @default: 1e-4 + photon_energy_min = "" + # Weights for the emitted inverse Compton photons + # @type: float [> 0.0] + # @default: 1.0 + photon_weight = "" + # Index of species for the emitted photon + # @type: ushort [> 0] + # @required + photon_species = "" + [algorithms] # Number of current smoothing passes # @type: ushort [>= 0] @@ -224,10 +277,12 @@ # @type: bool # @default: true enable = "" - # Order of the particle shape function - # @type: int - # @default: 1 - order = "" + + # @inferred: + # - order + # @brief: order of the particle shape function + # @from: compile-time definition `shape_order` + # @type: ushort [0 -> 10] [algorithms.gr] # Stepsize for numerical differentiation in GR pusher @@ -250,16 +305,13 @@ # @note: When `larmor_max` == 0, the limit is disabled larmor_max = "" - [algorithms.synchrotron] - # Radiation reaction limit gamma-factor for synchrotron - # @type: float [> 0.0] - # @default: 1.0 - # @note: [required] if one of the species has `cooling = "synchrotron"` - gamma_rad = "" - # Stencil coefficients for the field solver [notation as in Blinne+ (2018)] # @note: Standard Yee solver: `delta_i = beta_ij = 0.0` [algorithms.fieldsolver] + # Enable the fieldsolver + # @type: bool + # @default: true + enable = "" # delta_x coefficient (for `F_{i +/- 3/2, j, k}`) # @type: float # @default: 0.0 @@ -366,8 +418,17 @@ # Radiation reaction to use for the species # @type: string # @default: "None" - # @enum: "None", "Synchrotron" - cooling = "" + # @enum: "None", "Synchrotron", "Compton" + # @note: Can also be coma-separated combination, e.g., "Synchrotron,Compton" + # @note: Relevant radiation.drag parameters should also be provided + radiative_drag = "" + # Particle emission policy for the species + # @type: string + # @default: "None" + # @enum: "None", "Synchrotron", "Compton", "StrongFieldPP" + # @note: Only one emission mechanism allowed + # @note: Appropriate radiation drag flag will be applied automatically (unless explicitly set to "None") + emission = "" # Parameters for specific problem generators and setups [setup] diff --git a/pgens/accretion/pgen.hpp b/pgens/accretion/pgen.hpp index 2266f4e1..b88d8dfc 100644 --- a/pgens/accretion/pgen.hpp +++ b/pgens/accretion/pgen.hpp @@ -5,13 +5,13 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "arch/traits.h" #include "utils/numeric.h" #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" #include "archetypes/spatial_dist.h" +#include "archetypes/traits.h" #include "framework/domain/metadomain.h" #include "kernels/particle_moments.hpp" @@ -43,7 +43,8 @@ namespace user { TWO * metric.spin() * g_00); } - Inline auto bx1(const coord_t& x_Ph) const -> real_t { // at ( i , j + HALF ) + Inline auto bx1(const coord_t& x_Ph) const + -> real_t { // at ( i , j + HALF ) coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; metric.template convert(x_Ph, xi); @@ -61,7 +62,8 @@ namespace user { } } - Inline auto bx2(const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j ) + Inline auto bx2(const coord_t& x_Ph) const + -> real_t { // at ( i + HALF , j ) coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; metric.template convert(x_Ph, xi); @@ -197,11 +199,17 @@ namespace user { template struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines { traits::compatible_with::value }; + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; static constexpr auto metrics { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value + }; + static constexpr auto dimensions { + arch::traits::pgen::compatible_with::value }; - static constexpr auto dimensions { traits::compatible_with::value }; // for easy access to variables in the child class using arch::ProblemGenerator::D; diff --git a/pgens/magnetosphere/pgen.hpp b/pgens/magnetosphere/pgen.hpp index eab55ffb..706c3501 100644 --- a/pgens/magnetosphere/pgen.hpp +++ b/pgens/magnetosphere/pgen.hpp @@ -5,10 +5,10 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "arch/traits.h" #include "utils/numeric.h" #include "archetypes/problem_generator.h" +#include "archetypes/traits.h" #include "framework/domain/metadomain.h" #include @@ -87,11 +87,15 @@ namespace user { template struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines { traits::compatible_with::value }; + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; static constexpr auto metrics { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value + }; + static constexpr auto dimensions { + arch::traits::pgen::compatible_with::value }; - static constexpr auto dimensions { traits::compatible_with::value }; // for easy access to variables in the child class using arch::ProblemGenerator::D; diff --git a/pgens/pgen.hpp b/pgens/pgen.hpp index 78d8e648..d2d091eb 100644 --- a/pgens/pgen.hpp +++ b/pgens/pgen.hpp @@ -4,10 +4,10 @@ #include "enums.h" #include "global.h" -#include "arch/traits.h" #include "utils/formatting.h" #include "archetypes/problem_generator.h" +#include "archetypes/traits.h" #include "framework/domain/metadomain.h" #include @@ -19,18 +19,18 @@ namespace user { struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator static constexpr auto engines { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value }; static constexpr auto metrics { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value }; static constexpr auto dimensions { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value }; // for easy access to variables in the child class diff --git a/pgens/reconnection/pgen.hpp b/pgens/reconnection/pgen.hpp index b446c11f..a615451b 100644 --- a/pgens/reconnection/pgen.hpp +++ b/pgens/reconnection/pgen.hpp @@ -5,13 +5,13 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "arch/traits.h" #include "utils/numeric.h" #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" #include "archetypes/spatial_dist.h" +#include "archetypes/traits.h" #include "archetypes/utils.h" #include "framework/domain/metadomain.h" @@ -141,10 +141,14 @@ namespace user { template struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines { traits::compatible_with::value }; - static constexpr auto metrics { traits::compatible_with::value }; + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto metrics { + arch::traits::pgen::compatible_with::value + }; static constexpr auto dimensions { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value }; // for easy access to variables in the child class diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp index c3904f02..b8289c9a 100644 --- a/pgens/shock/pgen.hpp +++ b/pgens/shock/pgen.hpp @@ -4,12 +4,12 @@ #include "enums.h" #include "global.h" -#include "arch/traits.h" #include "utils/error.h" #include "utils/numeric.h" #include "archetypes/field_setter.h" #include "archetypes/problem_generator.h" +#include "archetypes/traits.h" #include "archetypes/utils.h" #include "framework/domain/metadomain.h" @@ -69,10 +69,14 @@ namespace user { template struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines { traits::compatible_with::value }; - static constexpr auto metrics { traits::compatible_with::value }; + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto metrics { + arch::traits::pgen::compatible_with::value + }; static constexpr auto dimensions { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value }; // for easy access to variables in the child class @@ -117,8 +121,8 @@ namespace user { return init_flds; } - auto FixFieldsConst(const bc_in&, - const em& comp) const -> std::pair { + auto FixFieldsConst(const bc_in&, const em& comp) const + -> std::pair { if (comp == em::ex1) { return { init_flds.ex1({ ZERO }), true }; } else if (comp == em::ex2) { diff --git a/pgens/streaming/pgen.hpp b/pgens/streaming/pgen.hpp index dca6cc31..9c87d7d9 100644 --- a/pgens/streaming/pgen.hpp +++ b/pgens/streaming/pgen.hpp @@ -5,11 +5,11 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "arch/traits.h" #include "utils/error.h" #include "utils/numeric.h" #include "archetypes/problem_generator.h" +#include "archetypes/traits.h" #include "archetypes/utils.h" #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" @@ -54,18 +54,20 @@ namespace user { struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines = traits::compatible_with::value; - static constexpr auto metrics = traits::compatible_with::value; + static constexpr auto engines = + arch::traits::pgen::compatible_with::value; + static constexpr auto metrics = + arch::traits::pgen::compatible_with::value; static constexpr auto dimensions = - traits::compatible_with::value; + arch::traits::pgen::compatible_with::value; // for easy access to variables in the child class using arch::ProblemGenerator::D; using arch::ProblemGenerator::C; using arch::ProblemGenerator::params; - prmvec_t drifts_in_x, drifts_in_y, drifts_in_z; - prmvec_t densities, temperatures; + prmvec_t drifts_in_x, drifts_in_y, drifts_in_z; + prmvec_t densities, temperatures; // initial magnetic field real_t Btheta, Bphi, Bmag; InitFields init_flds; diff --git a/pgens/turbulence/pgen.hpp b/pgens/turbulence/pgen.hpp index e8001b09..c7579512 100644 --- a/pgens/turbulence/pgen.hpp +++ b/pgens/turbulence/pgen.hpp @@ -9,8 +9,8 @@ #include "utils/numeric.h" #include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" +#include "archetypes/traits.h" #include "archetypes/utils.h" #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" @@ -293,9 +293,12 @@ namespace user { struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines = traits::compatible_with::value; - static constexpr auto metrics = traits::compatible_with::value; - static constexpr auto dimensions = traits::compatible_with::value; + static constexpr auto engines = + arch::traits::pgen::compatible_with::value; + static constexpr auto metrics = + arch::traits::pgen::compatible_with::value; + static constexpr auto dimensions = + arch::traits::pgen::compatible_with::value; // for easy access to variables in the child class using arch::ProblemGenerator::D; diff --git a/pgens/wald/pgen.hpp b/pgens/wald/pgen.hpp index 71ee905e..292f6f7f 100644 --- a/pgens/wald/pgen.hpp +++ b/pgens/wald/pgen.hpp @@ -5,21 +5,16 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "arch/traits.h" #include "utils/comparators.h" #include "utils/error.h" #include "utils/formatting.h" -#include "utils/log.h" #include "utils/numeric.h" -#include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" -#include "framework/domain/domain.h" +#include "archetypes/traits.h" #include "framework/domain/metadomain.h" #include -#include enum InitFieldGeometry { Wald, @@ -64,7 +59,8 @@ namespace user { TWO * metric.spin() * g_00); } - Inline auto bx1(const coord_t& x_Ph) const -> real_t { // at ( i , j + HALF ) + Inline auto bx1(const coord_t& x_Ph) const + -> real_t { // at ( i , j + HALF ) coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; metric.template convert(x_Ph, xi); @@ -82,7 +78,8 @@ namespace user { } } - Inline auto bx2(const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j ) + Inline auto bx2(const coord_t& x_Ph) const + -> real_t { // at ( i + HALF , j ) coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; metric.template convert(x_Ph, xi); @@ -99,8 +96,8 @@ namespace user { } } - Inline auto bx3( - const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j + HALF ) + Inline auto bx3(const coord_t& x_Ph) const + -> real_t { // at ( i + HALF , j + HALF ) if (field_geometry == InitFieldGeometry::Wald) { coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; metric.template convert(x_Ph, xi); @@ -120,7 +117,8 @@ namespace user { } } - Inline auto dx1(const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j ) + Inline auto dx1(const coord_t& x_Ph) const + -> real_t { // at ( i + HALF , j ) if (field_geometry == InitFieldGeometry::Wald) { coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; metric.template convert(x_Ph, xi); @@ -158,7 +156,8 @@ namespace user { } } - Inline auto dx2(const coord_t& x_Ph) const -> real_t { // at ( i , j + HALF ) + Inline auto dx2(const coord_t& x_Ph) const + -> real_t { // at ( i , j + HALF ) if (field_geometry == InitFieldGeometry::Wald) { coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; metric.template convert(x_Ph, xi); @@ -232,11 +231,17 @@ namespace user { template struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines { traits::compatible_with::value }; + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; static constexpr auto metrics { - traits::compatible_with::value + arch::traits::pgen::compatible_with::value + }; + static constexpr auto dimensions { + arch::traits::pgen::compatible_with::value }; - static constexpr auto dimensions { traits::compatible_with::value }; // for easy access to variables in the child class using arch::ProblemGenerator::D; diff --git a/src/archetypes/energy_dist.h b/src/archetypes/energy_dist.h index f387922d..230e4a93 100644 --- a/src/archetypes/energy_dist.h +++ b/src/archetypes/energy_dist.h @@ -22,11 +22,12 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "arch/traits.h" #include "utils/comparators.h" #include "utils/error.h" #include "utils/numeric.h" +#include "metrics/traits.h" + #include #include @@ -34,7 +35,7 @@ namespace arch { using namespace ntt; template - requires traits::metric::HasD + requires metric::traits::HasD struct EnergyDistribution { static constexpr auto D = M::Dim; diff --git a/src/archetypes/field_setter.h b/src/archetypes/field_setter.h index 56a24af1..d4fab7a3 100644 --- a/src/archetypes/field_setter.h +++ b/src/archetypes/field_setter.h @@ -27,36 +27,34 @@ #include "arch/traits.h" #include "utils/numeric.h" +#include "metrics/traits.h" + #include namespace arch { using namespace ntt; template - requires traits::metric::HasD && - ((S == SimEngine::SRPIC && traits::metric::HasConvert && - traits::metric::HasTransform_i) || - (S == SimEngine::GRPIC && traits::metric::HasConvert_i)) + requires metric::traits::HasD && + ((S == SimEngine::SRPIC && metric::traits::HasConvert && + metric::traits::HasTransform_i) || + (S == SimEngine::GRPIC && metric::traits::HasConvert_i)) && + (S == SimEngine::SRPIC && + (::traits::fieldsetter::HasEx1 || + ::traits::fieldsetter::HasEx2 || + ::traits::fieldsetter::HasEx3 || + ::traits::fieldsetter::HasBx1 || + ::traits::fieldsetter::HasBx2 || + ::traits::fieldsetter::HasBx3) || + (S == SimEngine::GRPIC && + (::traits::fieldsetter::HasDx1 && + ::traits::fieldsetter::HasDx2 && + ::traits::fieldsetter::HasDx3) || + (::traits::fieldsetter::HasBx1 && + ::traits::fieldsetter::HasBx2 && + ::traits::fieldsetter::HasBx3))) class SetEMFields_kernel { static constexpr Dimension D = M::Dim; - static constexpr bool defines_ex1 = traits::has_method::value; - static constexpr bool defines_ex2 = traits::has_method::value; - static constexpr bool defines_ex3 = traits::has_method::value; - static constexpr bool defines_bx1 = traits::has_method::value; - static constexpr bool defines_bx2 = traits::has_method::value; - static constexpr bool defines_bx3 = traits::has_method::value; - static constexpr bool defines_dx1 = traits::has_method::value; - static constexpr bool defines_dx2 = traits::has_method::value; - static constexpr bool defines_dx3 = traits::has_method::value; - - static_assert(defines_ex1 || defines_ex2 || defines_ex3 || defines_bx1 || - defines_bx2 || defines_bx3 || defines_dx1 || defines_dx2 || - defines_dx3, - "No field initializer defined"); - static_assert((S != SimEngine::GRPIC) || - (defines_dx1 == defines_dx2 && defines_dx2 == defines_dx3 && - defines_bx1 == defines_bx2 && defines_bx2 == defines_bx3), - "In GR mode, all components must be defined or none"); ndfield_t EM; const I finit; @@ -75,37 +73,37 @@ namespace arch { const auto i1_ = COORD(i1); coord_t x_Phys { ZERO }; if constexpr (S == SimEngine::SRPIC) { - if constexpr (defines_ex1) { + if constexpr (::traits::fieldsetter::HasEx1) { metric.template convert({ i1_ + HALF }, x_Phys); EM(i1, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( { i1_ + HALF }, finit.ex1(x_Phys)); } - if constexpr (defines_ex2) { + if constexpr (::traits::fieldsetter::HasEx2) { metric.template convert({ i1_ }, x_Phys); EM(i1, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( { i1_ }, finit.ex2(x_Phys)); } - if constexpr (defines_ex3) { + if constexpr (::traits::fieldsetter::HasEx3) { metric.template convert({ i1_ }, x_Phys); EM(i1, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( { i1_ }, finit.ex3(x_Phys)); } - if constexpr (defines_bx1) { + if constexpr (::traits::fieldsetter::HasBx1) { metric.template convert({ i1_ }, x_Phys); EM(i1, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( { i1_ }, finit.bx1(x_Phys)); } - if constexpr (defines_bx2) { + if constexpr (::traits::fieldsetter::HasBx2) { metric.template convert({ i1_ + HALF }, x_Phys); EM(i1, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( { i1_ + HALF }, finit.bx2(x_Phys)); } - if constexpr (defines_bx3) { + if constexpr (::traits::fieldsetter::HasBx3) { metric.template convert({ i1_ + HALF }, x_Phys); EM(i1, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( { i1_ + HALF }, @@ -126,37 +124,37 @@ namespace arch { // srpic if constexpr (S == SimEngine::SRPIC) { coord_t x_Phys { ZERO }; - if constexpr (defines_ex1) { + if constexpr (::traits::fieldsetter::HasEx1) { metric.template convert({ i1_ + HALF, i2_ }, x_Phys); EM(i1, i2, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( { i1_ + HALF, i2_ }, finit.ex1(x_Phys)); } - if constexpr (defines_ex2) { + if constexpr (::traits::fieldsetter::HasEx2) { metric.template convert({ i1_, i2_ + HALF }, x_Phys); EM(i1, i2, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( { i1_, i2_ + HALF }, finit.ex2(x_Phys)); } - if constexpr (defines_ex3) { + if constexpr (::traits::fieldsetter::HasEx3) { metric.template convert({ i1_, i2_ }, x_Phys); EM(i1, i2, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( { i1_, i2_ }, finit.ex3(x_Phys)); } - if constexpr (defines_bx1) { + if constexpr (::traits::fieldsetter::HasBx1) { metric.template convert({ i1_, i2_ + HALF }, x_Phys); EM(i1, i2, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( { i1_, i2_ + HALF }, finit.bx1(x_Phys)); } - if constexpr (defines_bx2) { + if constexpr (::traits::fieldsetter::HasBx2) { metric.template convert({ i1_ + HALF, i2_ }, x_Phys); EM(i1, i2, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( { i1_ + HALF, i2_ }, finit.bx2(x_Phys)); } - if constexpr (defines_bx3) { + if constexpr (::traits::fieldsetter::HasBx3) { metric.template convert({ i1_ + HALF, i2_ + HALF }, x_Phys); EM(i1, i2, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( @@ -165,7 +163,9 @@ namespace arch { } } else if constexpr (S == SimEngine::GRPIC) { // grpic - if constexpr (defines_dx1 && defines_dx2 && defines_dx3) { + if constexpr (::traits::fieldsetter::HasDx1 && + ::traits::fieldsetter::HasDx2 && + ::traits::fieldsetter::HasDx3) { const real_t x1_0 { metric.template convert<1, Crd::Cd, Crd::Ph>(i1_) }; const real_t x1_H { metric.template convert<1, Crd::Cd, Crd::Ph>( i1_ + HALF) }; @@ -182,7 +182,9 @@ namespace arch { EM(i1, i2, em::dx3) = finit.dx3({ x1_0, x2_0 }); } } - if constexpr (defines_bx1 && defines_bx2 && defines_bx3) { + if constexpr (::traits::fieldsetter::HasBx1 && + ::traits::fieldsetter::HasBx2 && + ::traits::fieldsetter::HasBx3) { const real_t x1_0 { metric.template convert<1, Crd::Cd, Crd::Ph>(i1_) }; const real_t x1_H { metric.template convert<1, Crd::Cd, Crd::Ph>( i1_ + HALF) }; @@ -215,28 +217,28 @@ namespace arch { coord_t x_Phys { ZERO }; if constexpr (S == SimEngine::SRPIC) { // srpic - if constexpr (defines_ex1) { + if constexpr (::traits::fieldsetter::HasEx1) { metric.template convert({ i1_ + HALF, i2_, i3_ }, x_Phys); EM(i1, i2, i3, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( { i1_ + HALF, i2_, i3_ }, finit.ex1(x_Phys)); } - if constexpr (defines_ex2) { + if constexpr (::traits::fieldsetter::HasEx2) { metric.template convert({ i1_, i2_ + HALF, i3_ }, x_Phys); EM(i1, i2, i3, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( { i1_, i2_ + HALF, i3_ }, finit.ex2(x_Phys)); } - if constexpr (defines_ex3) { + if constexpr (::traits::fieldsetter::HasEx3) { metric.template convert({ i1_, i2_, i3_ + HALF }, x_Phys); EM(i1, i2, i3, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( { i1_, i2_, i3_ + HALF }, finit.ex3(x_Phys)); } - if constexpr (defines_bx1) { + if constexpr (::traits::fieldsetter::HasBx1) { metric.template convert( { i1_, i2_ + HALF, i3_ + HALF }, x_Phys); @@ -244,7 +246,7 @@ namespace arch { { i1_, i2_ + HALF, i3_ + HALF }, finit.bx1(x_Phys)); } - if constexpr (defines_bx2) { + if constexpr (::traits::fieldsetter::HasBx2) { metric.template convert( { i1_ + HALF, i2_, i3_ + HALF }, x_Phys); @@ -252,7 +254,7 @@ namespace arch { { i1_ + HALF, i2_, i3_ + HALF }, finit.bx2(x_Phys)); } - if constexpr (defines_bx3) { + if constexpr (::traits::fieldsetter::HasBx3) { metric.template convert( { i1_ + HALF, i2_ + HALF, i3_ }, x_Phys); @@ -272,7 +274,9 @@ namespace arch { const real_t x3_H { metric.template convert<3, Crd::Cd, Crd::Ph>( i3_ + HALF) }; - if constexpr (defines_dx1 && defines_dx2 && defines_dx3) { + if constexpr (::traits::fieldsetter::HasDx1 && + ::traits::fieldsetter::HasDx2 && + ::traits::fieldsetter::HasDx3) { { // dx1 EM(i1, i2, i3, em::dx1) = finit.dx1({ x1_H, x2_0, x3_0 }); } @@ -283,7 +287,9 @@ namespace arch { EM(i1, i2, i3, em::dx3) = finit.dx3({ x1_0, x2_0, x3_H }); } } - if constexpr (defines_bx1 && defines_bx2 && defines_bx3) { + if constexpr (::traits::fieldsetter::HasBx1 && + ::traits::fieldsetter::HasBx2 && + ::traits::fieldsetter::HasBx3) { { // bx1 EM(i1, i2, i3, em::bx1) = finit.bx1({ x1_0, x2_H, x3_H }); } diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index 1543e265..a5239a02 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -22,6 +22,8 @@ #include "utils/error.h" #include "utils/numeric.h" +#include "metrics/traits.h" + #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" @@ -53,7 +55,7 @@ namespace arch { * - array_t: maximum coordinates of the region in computational coords */ template - requires traits::metric::HasD && traits::metric::HasConvert + requires metric::traits::HasD && metric::traits::HasConvert auto DeduceRegion(const Domain& domain, const boundaries_t& box) -> std::tuple, array_t> { if (not domain.mesh.Intersects(box)) { @@ -108,7 +110,7 @@ namespace arch { * - array_t: maximum coordinates of the region in computational coords */ template - requires traits::metric::HasD + requires metric::traits::HasD auto ComputeNumInject(const SimulationParams& params, const Domain& domain, real_t number_density, @@ -200,7 +202,7 @@ namespace arch { * @tparam ED2 Energy distribution type for species 2 */ template - requires traits::metric::HasD && traits::energydist::IsValid && + requires metric::traits::HasD && traits::energydist::IsValid && traits::energydist::IsValid inline void InjectUniform(const SimulationParams& params, Domain& domain, @@ -309,7 +311,7 @@ namespace arch { * @tparam SD Spatial distribution type */ template - requires traits::metric::HasD && traits::energydist::IsValid && + requires metric::traits::HasD && traits::energydist::IsValid && traits::energydist::IsValid && traits::spatialdist::IsValid inline void InjectNonUniform(const SimulationParams& params, Domain& domain, diff --git a/src/archetypes/problem_generator.h b/src/archetypes/problem_generator.h index efcec90a..fc22428e 100644 --- a/src/archetypes/problem_generator.h +++ b/src/archetypes/problem_generator.h @@ -23,15 +23,15 @@ #include "enums.h" #include "global.h" -#include "arch/traits.h" +#include "metrics/traits.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" namespace arch { using namespace ntt; template - requires traits::metric::HasD + requires metric::traits::HasD and metric::traits::HasCoordType struct ProblemGenerator { static constexpr Dimension D { M::Dim }; static constexpr Coord C { M::CoordType }; diff --git a/src/archetypes/spatial_dist.h b/src/archetypes/spatial_dist.h index be42c3ba..bf8abe9c 100644 --- a/src/archetypes/spatial_dist.h +++ b/src/archetypes/spatial_dist.h @@ -23,11 +23,13 @@ #include "utils/error.h" #include "utils/numeric.h" +#include "metrics/traits.h" + namespace arch { using namespace ntt; template - requires traits::metric::HasD + requires metric::traits::HasD struct SpatialDistribution { static constexpr auto D = M::Dim; diff --git a/src/archetypes/traits.h b/src/archetypes/traits.h new file mode 100644 index 00000000..a0ca199f --- /dev/null +++ b/src/archetypes/traits.h @@ -0,0 +1,159 @@ +/** + * @file archetypes/traits.h + * @brief Defines a set of traits to check if archetype classes satisfy certain conditions + * @implements + * - arch::traits::energydist::IsValid<> - checks if energy distribution class has required operator() + * - arch::traits::spatialdist::IsValid<> - checks if spatial distribution class has required operator() + * - arch::traits::pgen::check_compatibility<> - checks if problem generator is compatible with given enums + * - arch::traits::pgen::compatible_with<> - defines compatible enums for problem generator + * - arch::traits::pgen::HasD<> - checks if problem generator has Dim static member + * - arch::traits::pgen::HasInitFlds<> - checks if problem generator has init_flds member + * - arch::traits::pgen::HasInitPrtls<> - checks if problem generator has InitPrtls method + * - arch::traits::pgen::HasExtForce<> - checks if problem generator has ext_force member + * - arch::traits::pgen::HasExtCurrent<> - checks if problem generator has ext_current member + * - arch::traits::pgen::HasAtmFields<> - checks if problem generator has AtmFields method + * - arch::traits::pgen::HasMatchFields<> - checks if problem generator has MatchFields method + * - arch::traits::pgen::HasMatchFieldsInX1<> - checks if problem generator has MatchFieldsInX1 method + * - arch::traits::pgen::HasMatchFieldsInX2<> - checks if problem generator has MatchFieldsInX2 method + * - arch::traits::pgen::HasMatchFieldsInX3<> - checks if problem generator has MatchFieldsInX3 method + * - arch::traits::pgen::HasFixFieldsConst<> - checks if problem generator has FixFieldsConst method + * - arch::traits::pgen::HasCustomPostStep<> - checks if problem generator has CustomPostStep method + * - arch::traits::pgen::HasCustomFieldOutput<> - checks if problem generator has CustomFieldOutput method + * - arch::traits::pgen::HasCustomStatOutput<> - checks if problem generator has CustomStat method + * @namespaces: + * - arch::traits:: + */ +#ifndef ARCHETYPES_TRAITS_H +#define ARCHETYPES_TRAITS_H + +#include "global.h" + +#include "arch/kokkos_aliases.h" + +namespace arch { + namespace traits { + + namespace energydist { + + template + concept IsValid = requires(const ED& edist, + const coord_t& x_Ph, + vec_t& v) { + { edist(x_Ph, v) } -> std::same_as; + }; + + } // namespace energydist + + namespace spatialdist { + + template + concept IsValid = requires(const SD& sdist, const coord_t& x_Ph) { + { sdist(x_Ph) } -> std::convertible_to; + }; + + } // namespace spatialdist + + namespace pgen { + + // checking compat for the problem generator + engine + template + struct check_compatibility { + template + static constexpr bool value(std::integer_sequence) { + return ((Is == N) || ...); + } + }; + + template + struct compatible_with { + static constexpr auto value = std::integer_sequence {}; + }; + + template + concept HasD = requires { + { PG::D } -> std::convertible_to; + }; + + template + concept HasInitFlds = requires(const PG& pgen) { pgen.init_flds; }; + + template + concept HasInitPrtls = requires(PG& pgen, D& domain) { + { pgen.InitPrtls(domain) } -> std::same_as; + }; + + template + concept HasExtForce = requires(const PG& pgen) { pgen.ext_force; }; + + template + concept HasExtCurrent = requires(const PG& pgen) { pgen.ext_current; }; + + template + concept HasAtmFields = requires(const PG& pgen, simtime_t time) { + pgen.AtmFields(time); + }; + + template + concept HasMatchFields = requires(const PG& pgen, simtime_t time) { + pgen.MatchFields(time); + }; + + template + concept HasMatchFieldsInX1 = requires(const PG& pgen, simtime_t time) { + pgen.MatchFieldsInX1(time); + }; + + template + concept HasMatchFieldsInX2 = requires(const PG& pgen, simtime_t time) { + pgen.MatchFieldsInX2(time); + }; + + template + concept HasMatchFieldsInX3 = requires(const PG& pgen, simtime_t time) { + pgen.MatchFieldsInX3(time); + }; + + template + concept HasFixFieldsConst = requires(const PG& pgen, + const bc_in& bc, + const ntt::em& comp) { + { + pgen.FixFieldsConst(bc, comp) + } -> std::convertible_to>; + }; + + template + concept HasCustomPostStep = requires(PG& pgen, + timestep_t s, + simtime_t t, + D& domain) { + { pgen.CustomPostStep(s, t, domain) } -> std::same_as; + }; + + template + concept HasCustomFieldOutput = requires(PG& pgen, + const std::string& name, + ndfield_t& buff, + index_t idx, + timestep_t step, + simtime_t time, + const D& dom) { + { + pgen.CustomFieldOutput(name, buff, idx, step, time, dom) + } -> std::same_as; + }; + + template + concept HasCustomStatOutput = requires(PG& pgen, + const std::string& name, + timestep_t s, + simtime_t t, + const D& dom) { + { pgen.CustomStat(name, s, t, dom) } -> std::convertible_to; + }; + + } // namespace pgen + } // namespace traits +} // namespace arch + +#endif // ARCHETYPES_TRAITS_H diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h index b3eb24ae..bbb4f8fb 100644 --- a/src/archetypes/utils.h +++ b/src/archetypes/utils.h @@ -17,7 +17,7 @@ #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "framework/domain/domain.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include diff --git a/src/engines/CMakeLists.txt b/src/engines/CMakeLists.txt index 4cef1863..61295e43 100644 --- a/src/engines/CMakeLists.txt +++ b/src/engines/CMakeLists.txt @@ -4,10 +4,6 @@ # # @sources: # -# * engine_printer.cpp -# * engine_init.cpp -# * engine_run.cpp -# # @includes: # # * ../ @@ -31,10 +27,7 @@ # * mpi [optional] # ------------------------------ -set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}) -set(SOURCES ${SRC_DIR}/engine_printer.cpp ${SRC_DIR}/engine_init.cpp - ${SRC_DIR}/engine_run.cpp) -add_library(ntt_engines ${SOURCES}) +add_library(ntt_engines INTERFACE) set(libs ntt_global ntt_framework ntt_metrics ntt_archetypes ntt_kernels ntt_pgen) @@ -42,10 +35,9 @@ if(${output}) list(APPEND libs ntt_output) endif() add_dependencies(ntt_engines ${libs}) -target_link_libraries(ntt_engines PUBLIC ${libs}) -target_compile_definitions(ntt_engines PRIVATE PGEN=\"${PGEN}\") +target_link_libraries(ntt_engines INTERFACE ${libs}) +target_compile_definitions(ntt_engines INTERFACE PGEN=\"${PGEN}\") target_include_directories( ntt_engines - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../ INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/../) diff --git a/src/engines/engine.hpp b/src/engines/engine.hpp index 3117b082..910d0cc0 100644 --- a/src/engines/engine.hpp +++ b/src/engines/engine.hpp @@ -3,9 +3,6 @@ * @brief Base simulation class which just initializes the metadomain * @implements * - ntt::Engine<> - * @cpp: - * - engine_init.cpp - * - engine_printer.cpp * @namespaces: * - ntt:: * @macros: @@ -21,19 +18,41 @@ #include "enums.h" #include "global.h" -#include "arch/traits.h" -#include "utils/error.h" +#include "arch/directions.h" +#include "arch/mpi_aliases.h" +#include "utils/colors.h" +#include "utils/diag.h" +#include "utils/formatting.h" #include "utils/timer.h" #include "utils/toml.h" +#include "metrics/traits.h" + +#include "archetypes/field_setter.h" +#include "archetypes/traits.h" #include "framework/containers/species.h" +#include "framework/domain/domain.h" #include "framework/domain/metadomain.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "pgen.hpp" +#if defined(CUDA_ENABLED) + #include +#elif defined(HIP_ENABLED) + #include +#endif + +#if defined(OUTPUT_ENABLED) + #include +#endif + #include +#include +#include +#include + #if defined(OUTPUT_ENABLED) #include #include @@ -49,7 +68,13 @@ namespace ntt { template - concept IsCompatibleWithEngine = traits::metric::HasD; + concept IsCompatibleWithEngine = + metric::traits::HasD and + arch::traits::pgen::check_compatibility::value(user::PGen::engines) and + arch::traits::pgen::check_compatibility::value( + user::PGen::metrics) and + arch::traits::pgen::check_compatibility::value( + user::PGen::dimensions); template requires IsCompatibleWithEngine @@ -78,12 +103,6 @@ namespace ntt { timestep_t step; public: - static constexpr bool pgen_is_ok { - traits::check_compatibility::value(user::PGen::engines) and - traits::check_compatibility::value(user::PGen::metrics) and - traits::check_compatibility::value(user::PGen::dimensions) - }; - static constexpr Dimension D { M::Dim }; Engine(const SimulationParams& params) @@ -109,9 +128,7 @@ namespace ntt { , start_step { m_params.get("checkpoint.start_step") } , start_time { m_params.get("checkpoint.start_time") } , time { start_time } - , step { start_step } { - raise::ErrorIf(not pgen_is_ok, "Problem generator is not compatible with the picked engine/metric/dimension", HERE); - } + , step { start_step } {} ~Engine() = default; @@ -123,6 +140,652 @@ namespace ntt { void run(); }; + template + requires IsCompatibleWithEngine + void Engine::init() { + m_metadomain.InitStatsWriter(m_params, is_resuming); +#if defined(OUTPUT_ENABLED) + m_metadomain.InitWriter(&m_adios, m_params); + m_metadomain.InitCheckpointWriter(&m_adios, m_params); +#endif + logger::Checkpoint("Initializing Engine", HERE); + if (not is_resuming) { + // start a new simulation with initial conditions + logger::Checkpoint("Loading initial conditions", HERE); + if constexpr (arch::traits::pgen::HasInitFlds>) { + logger::Checkpoint("Initializing fields from problem generator", HERE); + m_metadomain.runOnLocalDomains([&](auto& loc_dom) { + Kokkos::parallel_for( + "InitFields", + loc_dom.mesh.rangeActiveCells(), + arch::SetEMFields_kernel { + loc_dom.fields.em, + m_pgen.init_flds, + loc_dom.mesh.metric }); + }); + } + if constexpr ( + arch::traits::pgen::HasInitPrtls, Domain>) { + logger::Checkpoint("Initializing particles from problem generator", HERE); + m_metadomain.runOnLocalDomains([&](auto& loc_dom) { + m_pgen.InitPrtls(loc_dom); + }); + } + } else { +#if defined(OUTPUT_ENABLED) + // read simulation data from the checkpoint + raise::ErrorIf( + m_params.template get("checkpoint.start_step") == 0, + "Resuming simulation from a checkpoint requires a valid start_step", + HERE); + logger::Checkpoint("Resuming simulation from a checkpoint", HERE); + m_metadomain.ContinueFromCheckpoint(&m_adios, m_params); +#else + raise::Error( + "Resuming simulation from a checkpoint requires -D output=ON", + HERE); +#endif + } + print_report(); + } + + namespace { + void add_header(std::string& report, + const std::vector& lines, + const std::vector& colors) { + report += fmt::format("%s╔%s╗%s\n", + color::BRIGHT_BLACK, + fmt::repeat("═", 58).c_str(), + color::RESET); + for (auto i { 0u }; i < lines.size(); ++i) { + report += fmt::format("%s║%s %s%s%s%s%s║%s\n", + color::BRIGHT_BLACK, + color::RESET, + colors[i], + lines[i].c_str(), + color::RESET, + fmt::repeat(" ", 57 - lines[i].size()).c_str(), + color::BRIGHT_BLACK, + color::RESET); + } + report += fmt::format("%s╚%s╝%s\n", + color::BRIGHT_BLACK, + fmt::repeat("═", 58).c_str(), + color::RESET); + } + + void add_category(std::string& report, unsigned short indent, const char* name) { + report += fmt::format("%s%s%s%s\n", + std::string(indent, ' ').c_str(), + color::BLUE, + name, + color::RESET); + } + + void add_subcategory(std::string& report, unsigned short indent, const char* name) { + report += fmt::format("%s%s-%s %s:\n", + std::string(indent, ' ').c_str(), + color::BRIGHT_BLACK, + color::RESET, + name); + } + + void add_label(std::string& report, unsigned short indent, const char* label) { + report += fmt::format("%s%s\n", std::string(indent, ' ').c_str(), label); + } + + template + void add_param(std::string& report, + unsigned short indent, + const char* name, + const char* format, + Args... args) { + report += fmt::format("%s%s-%s %s: %s%s%s\n", + std::string(indent, ' ').c_str(), + color::BRIGHT_BLACK, + color::RESET, + name, + color::BRIGHT_YELLOW, + fmt::format(format, args...).c_str(), + color::RESET); + } + + template + void add_unlabeled_param(std::string& report, + unsigned short indent, + const char* name, + const char* format, + Args... args) { + report += fmt::format("%s%s: %s%s%s\n", + std::string(indent, ' ').c_str(), + name, + color::BRIGHT_YELLOW, + fmt::format(format, args...).c_str(), + color::RESET); + } + + auto bytes_to_human_readable(std::size_t bytes) + -> std::pair { + const std::vector units { "B", "KB", "MB", "GB", "TB" }; + idx_t unit_idx = 0; + auto size = static_cast(bytes); + while ((size >= 1024.0) and (unit_idx < units.size() - 1)) { + size /= 1024.0; + ++unit_idx; + } + return { size, units[unit_idx] }; + } + } // namespace + + template + requires IsCompatibleWithEngine + void Engine::print_report() const { + const auto colored_stdout = m_params.template get( + "diagnostics.colored_stdout"); + std::string report = ""; + CallOnce( + [&](auto& metadomain, auto& params) { +#if defined(MPI_ENABLED) + int mpi_v, mpi_subv; + MPI_Get_version(&mpi_v, &mpi_subv); + const std::string mpi_version = fmt::format("%d.%d", mpi_v, mpi_subv); +#else // not MPI_ENABLED + const std::string mpi_version = "OFF"; +#endif // MPI_ENABLED + + const auto entity_version = "Entity v" + std::string(ENTITY_VERSION); + const auto hash = std::string(ENTITY_GIT_HASH); + const auto pgen = std::string(PGEN); + const auto nspec = metadomain.species_params().size(); + const auto precision = (sizeof(real_t) == 4) ? "single" : "double"; + +#if defined(__clang__) + const std::string ccx = "Clang/LLVM " __clang_version__; +#elif defined(__ICC) || defined(__INTEL_COMPILER) + const std::string ccx = "Intel ICC/ICPC " __VERSION__; +#elif defined(__GNUC__) || defined(__GNUG__) + const std::string ccx = "GNU GCC/G++ " __VERSION__; +#elif defined(__HP_cc) || defined(__HP_aCC) + const std::string ccx = "Hewlett-Packard C/aC++ " __HP_aCC; +#elif defined(__IBMC__) || defined(__IBMCPP__) + const std::string ccx = "IBM XL C/C++ " __IBMCPP__; +#elif defined(_MSC_VER) + const std::string ccx = "Microsoft Visual Studio " _MSC_VER; +#else + const std::string ccx = "Unknown compiler"; +#endif + std::string cpp_standard; + if (__cplusplus == 202101L) { + cpp_standard = "C++23"; + } else if (__cplusplus == 202002L) { + cpp_standard = "C++20"; + } else if (__cplusplus == 201703L) { + cpp_standard = "C++17"; + } else if (__cplusplus == 201402L) { + cpp_standard = "C++14"; + } else if (__cplusplus == 201103L) { + cpp_standard = "C++11"; + } else if (__cplusplus == 199711L) { + cpp_standard = "C++98"; + } else { + cpp_standard = "pre-standard " + std::to_string(__cplusplus); + } + +#if defined(CUDA_ENABLED) + int cuda_v; + cudaRuntimeGetVersion(&cuda_v); + const auto major { cuda_v / 1000 }; + const auto minor { cuda_v % 1000 / 10 }; + const auto patch { cuda_v % 10 }; + const auto cuda_version = fmt::format("%d.%d.%d", major, minor, patch); +#elif defined(HIP_ENABLED) + int hip_v; + auto status = hipDriverGetVersion(&hip_v); + raise::ErrorIf(status != hipSuccess, + "hipDriverGetVersion failed with error code %d", + HERE); + const auto major { hip_v / 10000000 }; + const auto minor { (hip_v % 10000000) / 100000 }; + const auto patch { hip_v % 100000 }; + const auto hip_version = fmt::format("%d.%d.%d", major, minor, patch); +#endif + + const auto kokkos_version = fmt::format("%d.%d.%d", + KOKKOS_VERSION / 10000, + KOKKOS_VERSION / 100 % 100, + KOKKOS_VERSION % 100); + +#if defined(OUTPUT_ENABLED) + const std::string adios2_version = fmt::format("%d.%d.%d", + ADIOS2_VERSION / 10000, + ADIOS2_VERSION / 100 % 100, + ADIOS2_VERSION % 100); +#else // not OUTPUT_ENABLED + const std::string adios2_version = "OFF"; +#endif + +#if defined(DEBUG) + const std::string dbg = "ON"; +#else // not DEBUG + const std::string dbg = "OFF"; +#endif + + report += "\n\n"; + add_header(report, { entity_version }, { color::BRIGHT_GREEN }); + report += "\n"; + + /* + * Backend + */ + add_category(report, 4, "Backend"); + add_param(report, 4, "Build hash", "%s", hash.c_str()); + add_param(report, 4, "CXX", "%s [%s]", ccx.c_str(), cpp_standard.c_str()); +#if defined(CUDA_ENABLED) + add_param(report, 4, "CUDA", "%s", cuda_version.c_str()); +#elif defined(HIP_VERSION) + add_param(report, 4, "HIP", "%s", hip_version.c_str()); +#endif + add_param(report, 4, "MPI", "%s", mpi_version.c_str()); +#if defined(MPI_ENABLED) && defined(DEVICE_ENABLED) + #if defined(GPU_AWARE_MPI) + const std::string gpu_aware_mpi = "ON"; + #else + const std::string gpu_aware_mpi = "OFF"; + #endif + add_param(report, 4, "GPU-aware MPI", "%s", gpu_aware_mpi.c_str()); +#endif + add_param(report, 4, "Kokkos", "%s", kokkos_version.c_str()); + add_param(report, 4, "ADIOS2", "%s", adios2_version.c_str()); + add_param(report, 4, "Precision", "%s", precision); + add_param(report, 4, "Debug", "%s", dbg.c_str()); + report += "\n"; + + /* + * Compilation flags + */ + add_category(report, 4, "Compilation flags"); +#if defined(SINGLE_PRECISION) + add_param(report, 4, "SINGLE_PRECISION", "%s", "ON"); +#else + add_param(report, 4, "SINGLE_PRECISION", "%s", "OFF"); +#endif + +#if defined(OUTPUT_ENABLED) + add_param(report, 4, "OUTPUT_ENABLED", "%s", "ON"); +#else + add_param(report, 4, "OUTPUT_ENABLED", "%s", "OFF"); +#endif + +#if defined(DEBUG) + add_param(report, 4, "DEBUG", "%s", "ON"); +#else + add_param(report, 4, "DEBUG", "%s", "OFF"); +#endif + +#if defined(CUDA_ENABLED) + add_param(report, 4, "CUDA_ENABLED", "%s", "ON"); +#else + add_param(report, 4, "CUDA_ENABLED", "%s", "OFF"); +#endif + +#if defined(HIP_ENABLED) + add_param(report, 4, "HIP_ENABLED", "%s", "ON"); +#else + add_param(report, 4, "HIP_ENABLED", "%s", "OFF"); +#endif + +#if defined(DEVICE_ENABLED) + add_param(report, 4, "DEVICE_ENABLED", "%s", "ON"); +#else + add_param(report, 4, "DEVICE_ENABLED", "%s", "OFF"); +#endif + +#if defined(MPI_ENABLED) + add_param(report, 4, "MPI_ENABLED", "%s", "ON"); +#else + add_param(report, 4, "MPI_ENABLED", "%s", "OFF"); +#endif + +#if defined(GPU_AWARE_MPI) + add_param(report, 4, "GPU_AWARE_MPI", "%s", "ON"); +#else + add_param(report, 4, "GPU_AWARE_MPI", "%s", "OFF"); +#endif + report += "\n"; + + /* + * Simulation configs + */ + add_category(report, 4, "Configuration"); + add_param(report, + 4, + "Name", + "%s", + params.template get("simulation.name").c_str()); + add_param(report, 4, "Problem generator", "%s", pgen.c_str()); + add_param(report, 4, "Engine", "%s", SimEngine(S).to_string()); + add_param(report, 4, "Metric", "%s", Metric(M::MetricType).to_string()); +#if SHAPE_ORDER == 0 + add_param(report, 4, "Deposit", "%s", "zigzag"); +#else + add_param(report, 4, "Deposit", "%s", "esirkepov"); + add_param(report, 4, "Interpolation order", "%i", SHAPE_ORDER); +#endif + add_param(report, 4, "Timestep [dt]", "%.3e", dt); + add_param(report, 4, "Runtime", "%.3e [%d steps]", runtime, max_steps); + report += "\n"; + add_category(report, 4, "Global domain"); + add_param(report, + 4, + "Resolution", + "%s", + params.template stringize("grid.resolution").c_str()); + add_param(report, + 4, + "Extent", + "%s", + params.template stringize("grid.extent").c_str()); + add_param(report, + 4, + "Fiducial cell size [dx0]", + "%.3e", + params.template get("scales.dx0")); + add_subcategory(report, 4, "Boundary conditions"); + add_param( + report, + 6, + "Fields", + "%s", + params.template stringize("grid.boundaries.fields").c_str()); + add_param( + report, + 6, + "Particles", + "%s", + params.template stringize("grid.boundaries.particles").c_str()); + add_param(report, + 4, + "Domain decomposition", + "%s [%d total]", + fmt::formatVector(m_metadomain.ndomains_per_dim()).c_str(), + m_metadomain.ndomains()); + report += "\n"; + add_category(report, 4, "Fiducial parameters"); + add_param(report, + 4, + "Particles per cell [ppc0]", + "%.1f", + params.template get("particles.ppc0")); + add_param(report, + 4, + "Larmor radius [larmor0]", + "%.3e [%.3f dx0]", + params.template get("scales.larmor0"), + params.template get("scales.larmor0") / + params.template get("scales.dx0")); + add_param(report, + 4, + "Larmor frequency [omegaB0 * dt]", + "%.3e", + params.template get("scales.omegaB0") * + params.template get("algorithms.timestep.dt")); + add_param(report, + 4, + "Skin depth [skindepth0]", + "%.3e [%.3f dx0]", + params.template get("scales.skindepth0"), + params.template get("scales.skindepth0") / + params.template get("scales.dx0")); + add_param(report, + 4, + "Plasma frequency [omp0 * dt]", + "%.3e", + params.template get("algorithms.timestep.dt") / + params.template get("scales.skindepth0")); + add_param(report, + 4, + "Magnetization [sigma0]", + "%.3e", + params.template get("scales.sigma0")); + + if (nspec > 0) { + report += "\n"; + add_category(report, 4, "Particles"); + } + for (const auto& species : metadomain.species_params()) { + add_subcategory(report, + 4, + fmt::format("Species #%d", species.index()).c_str()); + add_param(report, 6, "Label", "%s", species.label().c_str()); + add_param(report, 6, "Mass", "%.1f", species.mass()); + add_param(report, 6, "Charge", "%.1f", species.charge()); + add_param(report, 6, "Max #", "%d [per domain]", species.maxnpart()); + add_param(report, + 6, + "Pusher", + "%s", + ParticlePusher::to_string(species.pusher()).c_str()); + add_param( + report, + 6, + "Radiative drag", + "%s", + RadiativeDrag::to_string(species.radiative_drag_flags()).c_str()); + add_param(report, 6, "# of real-value payloads", "%d", species.npld_r()); + add_param(report, 6, "# of integer-value payloads", "%d", species.npld_i()); + } + report.pop_back(); + }, + m_metadomain, + m_params); + info::Print(report, colored_stdout); + + report = "\n"; + CallOnce([&]() { + add_category(report, 4, "Domains"); + report.pop_back(); + }); + info::Print(report, colored_stdout); + + for (unsigned int idx { 0 }; idx < m_metadomain.ndomains(); ++idx) { + auto is_local = false; + for (const auto& lidx : m_metadomain.l_subdomain_indices()) { + is_local |= (idx == lidx); + } + if (is_local) { + report = ""; + const auto& domain = m_metadomain.subdomain(idx); + add_subcategory(report, + 4, + fmt::format("Domain #%d", domain.index()).c_str()); +#if defined(MPI_ENABLED) + add_param(report, 6, "Rank", "%d", domain.mpi_rank()); +#endif + add_param(report, + 6, + "Resolution", + "%s", + fmt::formatVector(domain.mesh.n_active()).c_str()); + add_param(report, + 6, + "Extent", + "%s", + fmt::formatVector(domain.mesh.extent()).c_str()); + add_subcategory(report, 6, "Boundary conditions"); + + add_label( + report, + 8 + 2 + 2 * M::Dim, + fmt::format("%-10s %-10s %-10s", "[flds]", "[prtl]", "[neighbor]").c_str()); + for (auto& direction : dir::Directions::all) { + const auto flds_bc = domain.mesh.flds_bc_in(direction); + const auto prtl_bc = domain.mesh.prtl_bc_in(direction); + bool has_sync = false; + auto neighbor_idx = domain.neighbor_idx_in(direction); + if (flds_bc == FldsBC::SYNC || prtl_bc == PrtlBC::SYNC) { + has_sync = true; + } + add_unlabeled_param(report, + 8, + direction.to_string().c_str(), + "%-10s %-10s %-10s", + flds_bc.to_string(), + prtl_bc.to_string(), + has_sync ? std::to_string(neighbor_idx).c_str() + : "."); + } + add_subcategory(report, 6, "Memory footprint"); + auto flds_footprint = domain.fields.memory_footprint(); + auto [flds_size, flds_unit] = bytes_to_human_readable(flds_footprint); + add_param(report, 8, "Fields", "%.2f %s", flds_size, flds_unit.c_str()); + if (domain.species.size() > 0) { + add_subcategory(report, 8, "Particles"); + } + for (auto& species : domain.species) { + const auto str = fmt::format("Species #%d (%s)", + species.index(), + species.label().c_str()); + auto [size, unit] = bytes_to_human_readable(species.memory_footprint()); + add_param(report, 10, str.c_str(), "%.2f %s", size, unit.c_str()); + } + report.pop_back(); + if (idx == m_metadomain.ndomains() - 1) { + report += "\n\n"; + } + info::Print(report, colored_stdout, true, false); + } +#if defined(MPI_ENABLED) + MPI_Barrier(MPI_COMM_WORLD); +#endif + } + } + + template + requires IsCompatibleWithEngine + void Engine::run() { + init(); + + auto timers = timer::Timers { + { "FieldSolver", + "CurrentFiltering", "CurrentDeposit", + "ParticlePusher", "FieldBoundaries", + "ParticleBoundaries", "Communications", + "Injector", "Custom", + "PrtlClear", "Output", + "Checkpoint" }, + []() { + Kokkos::fence(); + }, + m_params.get("diagnostics.blocking_timers") + }; + const auto diag_interval = m_params.template get( + "diagnostics.interval"); + + auto time_history = pbar::DurationHistory { 1000 }; + const auto clear_interval = m_params.template get( + "particles.clear_interval"); + + // main algorithm loop + while (step < max_steps) { + // run the engine-dependent algorithm step + m_metadomain.runOnLocalDomains([&timers, this](auto& dom) { + step_forward(timers, dom); + }); + // poststep (if defined) + if constexpr ( + arch::traits::pgen::HasCustomPostStep>) { + timers.start("Custom"); + m_metadomain.runOnLocalDomains([&timers, this](auto& dom) { + m_pgen.CustomPostStep(step, time, dom); + }); + timers.stop("Custom"); + } + auto print_prtl_clear = (clear_interval > 0 and + step % clear_interval == 0 and step > 0); + + // advance time & step + time += dt; + ++step; + + auto print_output = false; + auto print_checkpoint = false; +#if defined(OUTPUT_ENABLED) + timers.start("Output"); + if constexpr ( + arch::traits::pgen::HasCustomFieldOutput>) { + auto lambda_custom_field_output = [&](const std::string& name, + ndfield_t& buff, + index_t idx, + timestep_t step, + simtime_t time, + const Domain& dom) { + m_pgen.CustomFieldOutput(name, buff, idx, step, time, dom); + }; + print_output &= m_metadomain.Write(m_params, + step, + step - 1, + time, + time - dt, + lambda_custom_field_output); + } else { + print_output &= m_metadomain.Write(m_params, step, step - 1, time, time - dt); + } + if constexpr ( + arch::traits::pgen::HasCustomStatOutput>) { + auto lambda_custom_stat = [&](const std::string& name, + timestep_t step, + simtime_t time, + const Domain& dom) -> real_t { + return m_pgen.CustomStat(name, step, time, dom); + }; + print_output &= m_metadomain.WriteStats(m_params, + step, + step - 1, + time, + time - dt, + lambda_custom_stat); + } else { + print_output &= m_metadomain.WriteStats(m_params, + step, + step - 1, + time, + time - dt); + } + timers.stop("Output"); + + timers.start("Checkpoint"); + print_checkpoint = m_metadomain.WriteCheckpoint(m_params, + step, + step - 1, + time, + time - dt); + timers.stop("Checkpoint"); +#endif + + // advance time_history + time_history.tick(); + // print timestep report + if (diag_interval > 0 and step % diag_interval == 0) { + diag::printDiagnostics( + step - 1, + max_steps, + time - dt, + dt, + timers, + time_history, + m_metadomain.l_ncells(), + m_metadomain.species_labels(), + m_metadomain.l_npart_perspec(), + m_metadomain.l_maxnpart_perspec(), + print_prtl_clear, + print_output, + print_checkpoint, + m_params.get("diagnostics.colored_stdout")); + } + timers.resetAll(); + } + } + } // namespace ntt #endif // ENGINES_ENGINE_H diff --git a/src/engines/engine_init.cpp b/src/engines/engine_init.cpp deleted file mode 100644 index e822ba3f..00000000 --- a/src/engines/engine_init.cpp +++ /dev/null @@ -1,73 +0,0 @@ -#include "enums.h" -#include "global.h" - -#include "arch/traits.h" - -#include "archetypes/field_setter.h" -#include "framework/specialization_registry.h" - -#include "engines/engine.hpp" - -#include - -#include - -namespace ntt { - - template - requires IsCompatibleWithEngine - void Engine::init() { - if constexpr (pgen_is_ok) { - m_metadomain.InitStatsWriter(m_params, is_resuming); -#if defined(OUTPUT_ENABLED) - m_metadomain.InitWriter(&m_adios, m_params); - m_metadomain.InitCheckpointWriter(&m_adios, m_params); -#endif - logger::Checkpoint("Initializing Engine", HERE); - if (not is_resuming) { - // start a new simulation with initial conditions - logger::Checkpoint("Loading initial conditions", HERE); - if constexpr (traits::pgen::HasInitFlds>) { - logger::Checkpoint("Initializing fields from problem generator", HERE); - m_metadomain.runOnLocalDomains([&](auto& loc_dom) { - Kokkos::parallel_for( - "InitFields", - loc_dom.mesh.rangeActiveCells(), - arch::SetEMFields_kernel { - loc_dom.fields.em, - m_pgen.init_flds, - loc_dom.mesh.metric }); - }); - } - if constexpr (traits::pgen::HasInitPrtls, Domain>) { - logger::Checkpoint("Initializing particles from problem generator", HERE); - m_metadomain.runOnLocalDomains([&](auto& loc_dom) { - m_pgen.InitPrtls(loc_dom); - }); - } - } else { -#if defined(OUTPUT_ENABLED) - // read simulation data from the checkpoint - raise::ErrorIf( - m_params.template get("checkpoint.start_step") == 0, - "Resuming simulation from a checkpoint requires a valid start_step", - HERE); - logger::Checkpoint("Resuming simulation from a checkpoint", HERE); - m_metadomain.ContinueFromCheckpoint(&m_adios, m_params); -#else - raise::Error( - "Resuming simulation from a checkpoint requires -D output=ON", - HERE); -#endif - } - } - print_report(); - } - -#define ENGINE_INIT(S, M, D) template class Engine>; - - NTT_FOREACH_SPECIALIZATION(ENGINE_INIT) - -#undef ENGINE_INIT - -} // namespace ntt diff --git a/src/engines/engine_printer.cpp b/src/engines/engine_printer.cpp deleted file mode 100644 index 34df741c..00000000 --- a/src/engines/engine_printer.cpp +++ /dev/null @@ -1,501 +0,0 @@ -#include "enums.h" -#include "global.h" - -#include "arch/directions.h" -#include "arch/mpi_aliases.h" -#include "utils/colors.h" -#include "utils/formatting.h" - -#include "framework/specialization_registry.h" - -#include "engines/engine.hpp" - -#if defined(CUDA_ENABLED) - #include -#elif defined(HIP_ENABLED) - #include -#endif - -#if defined(OUTPUT_ENABLED) - #include -#endif - -#include -#include -#include - -namespace ntt { - - namespace { - void add_header(std::string& report, - const std::vector& lines, - const std::vector& colors) { - report += fmt::format("%s╔%s╗%s\n", - color::BRIGHT_BLACK, - fmt::repeat("═", 58).c_str(), - color::RESET); - for (auto i { 0u }; i < lines.size(); ++i) { - report += fmt::format("%s║%s %s%s%s%s%s║%s\n", - color::BRIGHT_BLACK, - color::RESET, - colors[i], - lines[i].c_str(), - color::RESET, - fmt::repeat(" ", 57 - lines[i].size()).c_str(), - color::BRIGHT_BLACK, - color::RESET); - } - report += fmt::format("%s╚%s╝%s\n", - color::BRIGHT_BLACK, - fmt::repeat("═", 58).c_str(), - color::RESET); - } - - void add_category(std::string& report, unsigned short indent, const char* name) { - report += fmt::format("%s%s%s%s\n", - std::string(indent, ' ').c_str(), - color::BLUE, - name, - color::RESET); - } - - void add_subcategory(std::string& report, unsigned short indent, const char* name) { - report += fmt::format("%s%s-%s %s:\n", - std::string(indent, ' ').c_str(), - color::BRIGHT_BLACK, - color::RESET, - name); - } - - void add_label(std::string& report, unsigned short indent, const char* label) { - report += fmt::format("%s%s\n", std::string(indent, ' ').c_str(), label); - } - - template - void add_param(std::string& report, - unsigned short indent, - const char* name, - const char* format, - Args... args) { - report += fmt::format("%s%s-%s %s: %s%s%s\n", - std::string(indent, ' ').c_str(), - color::BRIGHT_BLACK, - color::RESET, - name, - color::BRIGHT_YELLOW, - fmt::format(format, args...).c_str(), - color::RESET); - } - - template - void add_unlabeled_param(std::string& report, - unsigned short indent, - const char* name, - const char* format, - Args... args) { - report += fmt::format("%s%s: %s%s%s\n", - std::string(indent, ' ').c_str(), - name, - color::BRIGHT_YELLOW, - fmt::format(format, args...).c_str(), - color::RESET); - } - - auto bytes_to_human_readable( - std::size_t bytes) -> std::pair { - const std::vector units { "B", "KB", "MB", "GB", "TB" }; - idx_t unit_idx = 0; - auto size = static_cast(bytes); - while ((size >= 1024.0) and (unit_idx < units.size() - 1)) { - size /= 1024.0; - ++unit_idx; - } - return { size, units[unit_idx] }; - } - } // namespace - - template - requires IsCompatibleWithEngine - void Engine::print_report() const { - const auto colored_stdout = m_params.template get( - "diagnostics.colored_stdout"); - std::string report = ""; - CallOnce( - [&](auto& metadomain, auto& params) { -#if defined(MPI_ENABLED) - int mpi_v, mpi_subv; - MPI_Get_version(&mpi_v, &mpi_subv); - const std::string mpi_version = fmt::format("%d.%d", mpi_v, mpi_subv); -#else // not MPI_ENABLED - const std::string mpi_version = "OFF"; -#endif // MPI_ENABLED - - const auto entity_version = "Entity v" + std::string(ENTITY_VERSION); - const auto hash = std::string(ENTITY_GIT_HASH); - const auto pgen = std::string(PGEN); - const auto nspec = metadomain.species_params().size(); - const auto precision = (sizeof(real_t) == 4) ? "single" : "double"; - -#if defined(__clang__) - const std::string ccx = "Clang/LLVM " __clang_version__; -#elif defined(__ICC) || defined(__INTEL_COMPILER) - const std::string ccx = "Intel ICC/ICPC " __VERSION__; -#elif defined(__GNUC__) || defined(__GNUG__) - const std::string ccx = "GNU GCC/G++ " __VERSION__; -#elif defined(__HP_cc) || defined(__HP_aCC) - const std::string ccx = "Hewlett-Packard C/aC++ " __HP_aCC; -#elif defined(__IBMC__) || defined(__IBMCPP__) - const std::string ccx = "IBM XL C/C++ " __IBMCPP__; -#elif defined(_MSC_VER) - const std::string ccx = "Microsoft Visual Studio " _MSC_VER; -#else - const std::string ccx = "Unknown compiler"; -#endif - std::string cpp_standard; - if (__cplusplus == 202101L) { - cpp_standard = "C++23"; - } else if (__cplusplus == 202002L) { - cpp_standard = "C++20"; - } else if (__cplusplus == 201703L) { - cpp_standard = "C++17"; - } else if (__cplusplus == 201402L) { - cpp_standard = "C++14"; - } else if (__cplusplus == 201103L) { - cpp_standard = "C++11"; - } else if (__cplusplus == 199711L) { - cpp_standard = "C++98"; - } else { - cpp_standard = "pre-standard " + std::to_string(__cplusplus); - } - -#if defined(CUDA_ENABLED) - int cuda_v; - cudaRuntimeGetVersion(&cuda_v); - const auto major { cuda_v / 1000 }; - const auto minor { cuda_v % 1000 / 10 }; - const auto patch { cuda_v % 10 }; - const auto cuda_version = fmt::format("%d.%d.%d", major, minor, patch); -#elif defined(HIP_ENABLED) - int hip_v; - auto status = hipDriverGetVersion(&hip_v); - raise::ErrorIf(status != hipSuccess, - "hipDriverGetVersion failed with error code %d", - HERE); - const auto major { hip_v / 10000000 }; - const auto minor { (hip_v % 10000000) / 100000 }; - const auto patch { hip_v % 100000 }; - const auto hip_version = fmt::format("%d.%d.%d", major, minor, patch); -#endif - - const auto kokkos_version = fmt::format("%d.%d.%d", - KOKKOS_VERSION / 10000, - KOKKOS_VERSION / 100 % 100, - KOKKOS_VERSION % 100); - -#if defined(OUTPUT_ENABLED) - const std::string adios2_version = fmt::format("%d.%d.%d", - ADIOS2_VERSION / 10000, - ADIOS2_VERSION / 100 % 100, - ADIOS2_VERSION % 100); -#else // not OUTPUT_ENABLED - const std::string adios2_version = "OFF"; -#endif - -#if defined(DEBUG) - const std::string dbg = "ON"; -#else // not DEBUG - const std::string dbg = "OFF"; -#endif - - report += "\n\n"; - add_header(report, { entity_version }, { color::BRIGHT_GREEN }); - report += "\n"; - - /* - * Backend - */ - add_category(report, 4, "Backend"); - add_param(report, 4, "Build hash", "%s", hash.c_str()); - add_param(report, 4, "CXX", "%s [%s]", ccx.c_str(), cpp_standard.c_str()); -#if defined(CUDA_ENABLED) - add_param(report, 4, "CUDA", "%s", cuda_version.c_str()); -#elif defined(HIP_VERSION) - add_param(report, 4, "HIP", "%s", hip_version.c_str()); -#endif - add_param(report, 4, "MPI", "%s", mpi_version.c_str()); -#if defined(MPI_ENABLED) && defined(DEVICE_ENABLED) - #if defined(GPU_AWARE_MPI) - const std::string gpu_aware_mpi = "ON"; - #else - const std::string gpu_aware_mpi = "OFF"; - #endif - add_param(report, 4, "GPU-aware MPI", "%s", gpu_aware_mpi.c_str()); -#endif - add_param(report, 4, "Kokkos", "%s", kokkos_version.c_str()); - add_param(report, 4, "ADIOS2", "%s", adios2_version.c_str()); - add_param(report, 4, "Precision", "%s", precision); - add_param(report, 4, "Debug", "%s", dbg.c_str()); - report += "\n"; - - /* - * Compilation flags - */ - add_category(report, 4, "Compilation flags"); -#if defined(SINGLE_PRECISION) - add_param(report, 4, "SINGLE_PRECISION", "%s", "ON"); -#else - add_param(report, 4, "SINGLE_PRECISION", "%s", "OFF"); -#endif - -#if defined(OUTPUT_ENABLED) - add_param(report, 4, "OUTPUT_ENABLED", "%s", "ON"); -#else - add_param(report, 4, "OUTPUT_ENABLED", "%s", "OFF"); -#endif - -#if defined(DEBUG) - add_param(report, 4, "DEBUG", "%s", "ON"); -#else - add_param(report, 4, "DEBUG", "%s", "OFF"); -#endif - -#if defined(CUDA_ENABLED) - add_param(report, 4, "CUDA_ENABLED", "%s", "ON"); -#else - add_param(report, 4, "CUDA_ENABLED", "%s", "OFF"); -#endif - -#if defined(HIP_ENABLED) - add_param(report, 4, "HIP_ENABLED", "%s", "ON"); -#else - add_param(report, 4, "HIP_ENABLED", "%s", "OFF"); -#endif - -#if defined(DEVICE_ENABLED) - add_param(report, 4, "DEVICE_ENABLED", "%s", "ON"); -#else - add_param(report, 4, "DEVICE_ENABLED", "%s", "OFF"); -#endif - -#if defined(MPI_ENABLED) - add_param(report, 4, "MPI_ENABLED", "%s", "ON"); -#else - add_param(report, 4, "MPI_ENABLED", "%s", "OFF"); -#endif - -#if defined(GPU_AWARE_MPI) - add_param(report, 4, "GPU_AWARE_MPI", "%s", "ON"); -#else - add_param(report, 4, "GPU_AWARE_MPI", "%s", "OFF"); -#endif - report += "\n"; - - /* - * Simulation configs - */ - add_category(report, 4, "Configuration"); - add_param(report, - 4, - "Name", - "%s", - params.template get("simulation.name").c_str()); - add_param(report, 4, "Problem generator", "%s", pgen.c_str()); - add_param(report, 4, "Engine", "%s", SimEngine(S).to_string()); - add_param(report, 4, "Metric", "%s", Metric(M::MetricType).to_string()); -#if SHAPE_ORDER == 0 - add_param(report, 4, "Deposit", "%s", "zigzag"); -#else - add_param(report, 4, "Deposit", "%s", "esirkepov"); - add_param(report, 4, "Interpolation order", "%i", SHAPE_ORDER); -#endif - add_param(report, 4, "Timestep [dt]", "%.3e", dt); - add_param(report, 4, "Runtime", "%.3e [%d steps]", runtime, max_steps); - report += "\n"; - add_category(report, 4, "Global domain"); - add_param(report, - 4, - "Resolution", - "%s", - params.template stringize("grid.resolution").c_str()); - add_param(report, - 4, - "Extent", - "%s", - params.template stringize("grid.extent").c_str()); - add_param(report, - 4, - "Fiducial cell size [dx0]", - "%.3e", - params.template get("scales.dx0")); - add_subcategory(report, 4, "Boundary conditions"); - add_param( - report, - 6, - "Fields", - "%s", - params.template stringize("grid.boundaries.fields").c_str()); - add_param( - report, - 6, - "Particles", - "%s", - params.template stringize("grid.boundaries.particles").c_str()); - add_param(report, - 4, - "Domain decomposition", - "%s [%d total]", - fmt::formatVector(m_metadomain.ndomains_per_dim()).c_str(), - m_metadomain.ndomains()); - report += "\n"; - add_category(report, 4, "Fiducial parameters"); - add_param(report, - 4, - "Particles per cell [ppc0]", - "%.1f", - params.template get("particles.ppc0")); - add_param(report, - 4, - "Larmor radius [larmor0]", - "%.3e [%.3f dx0]", - params.template get("scales.larmor0"), - params.template get("scales.larmor0") / - params.template get("scales.dx0")); - add_param(report, - 4, - "Larmor frequency [omegaB0 * dt]", - "%.3e", - params.template get("scales.omegaB0") * - params.template get("algorithms.timestep.dt")); - add_param(report, - 4, - "Skin depth [skindepth0]", - "%.3e [%.3f dx0]", - params.template get("scales.skindepth0"), - params.template get("scales.skindepth0") / - params.template get("scales.dx0")); - add_param(report, - 4, - "Plasma frequency [omp0 * dt]", - "%.3e", - params.template get("algorithms.timestep.dt") / - params.template get("scales.skindepth0")); - add_param(report, - 4, - "Magnetization [sigma0]", - "%.3e", - params.template get("scales.sigma0")); - - if (nspec > 0) { - report += "\n"; - add_category(report, 4, "Particles"); - } - for (const auto& species : metadomain.species_params()) { - add_subcategory(report, - 4, - fmt::format("Species #%d", species.index()).c_str()); - add_param(report, 6, "Label", "%s", species.label().c_str()); - add_param(report, 6, "Mass", "%.1f", species.mass()); - add_param(report, 6, "Charge", "%.1f", species.charge()); - add_param(report, 6, "Max #", "%d [per domain]", species.maxnpart()); - add_param(report, 6, "Pusher", "%s", species.pusher().to_string()); - if (species.mass() != 0.0) { - add_param(report, 6, "GCA", "%s", species.use_gca() ? "ON" : "OFF"); - } - add_param(report, 6, "Cooling", "%s", species.cooling().to_string()); - add_param(report, 6, "# of real-value payloads", "%d", species.npld_r()); - add_param(report, 6, "# of integer-value payloads", "%d", species.npld_i()); - } - report.pop_back(); - }, - m_metadomain, - m_params); - info::Print(report, colored_stdout); - - report = "\n"; - CallOnce([&]() { - add_category(report, 4, "Domains"); - report.pop_back(); - }); - info::Print(report, colored_stdout); - - for (unsigned int idx { 0 }; idx < m_metadomain.ndomains(); ++idx) { - auto is_local = false; - for (const auto& lidx : m_metadomain.l_subdomain_indices()) { - is_local |= (idx == lidx); - } - if (is_local) { - report = ""; - const auto& domain = m_metadomain.subdomain(idx); - add_subcategory(report, - 4, - fmt::format("Domain #%d", domain.index()).c_str()); -#if defined(MPI_ENABLED) - add_param(report, 6, "Rank", "%d", domain.mpi_rank()); -#endif - add_param(report, - 6, - "Resolution", - "%s", - fmt::formatVector(domain.mesh.n_active()).c_str()); - add_param(report, - 6, - "Extent", - "%s", - fmt::formatVector(domain.mesh.extent()).c_str()); - add_subcategory(report, 6, "Boundary conditions"); - - add_label( - report, - 8 + 2 + 2 * M::Dim, - fmt::format("%-10s %-10s %-10s", "[flds]", "[prtl]", "[neighbor]").c_str()); - for (auto& direction : dir::Directions::all) { - const auto flds_bc = domain.mesh.flds_bc_in(direction); - const auto prtl_bc = domain.mesh.prtl_bc_in(direction); - bool has_sync = false; - auto neighbor_idx = domain.neighbor_idx_in(direction); - if (flds_bc == FldsBC::SYNC || prtl_bc == PrtlBC::SYNC) { - has_sync = true; - } - add_unlabeled_param(report, - 8, - direction.to_string().c_str(), - "%-10s %-10s %-10s", - flds_bc.to_string(), - prtl_bc.to_string(), - has_sync ? std::to_string(neighbor_idx).c_str() - : "."); - } - add_subcategory(report, 6, "Memory footprint"); - auto flds_footprint = domain.fields.memory_footprint(); - auto [flds_size, flds_unit] = bytes_to_human_readable(flds_footprint); - add_param(report, 8, "Fields", "%.2f %s", flds_size, flds_unit.c_str()); - if (domain.species.size() > 0) { - add_subcategory(report, 8, "Particles"); - } - for (auto& species : domain.species) { - const auto str = fmt::format("Species #%d (%s)", - species.index(), - species.label().c_str()); - auto [size, unit] = bytes_to_human_readable(species.memory_footprint()); - add_param(report, 10, str.c_str(), "%.2f %s", size, unit.c_str()); - } - report.pop_back(); - if (idx == m_metadomain.ndomains() - 1) { - report += "\n\n"; - } - info::Print(report, colored_stdout, true, false); - } -#if defined(MPI_ENABLED) - MPI_Barrier(MPI_COMM_WORLD); -#endif - } - } - -#define ENGINE_PRINTER(S, M, D) \ - template void Engine>::print_report() const; - - NTT_FOREACH_SPECIALIZATION(ENGINE_PRINTER) - -#undef ENGINE_PRINTER - -} // namespace ntt diff --git a/src/engines/engine_run.cpp b/src/engines/engine_run.cpp deleted file mode 100644 index 7d81fe4b..00000000 --- a/src/engines/engine_run.cpp +++ /dev/null @@ -1,147 +0,0 @@ -#include "enums.h" - -#include "arch/traits.h" -#include "utils/diag.h" - -#include "framework/domain/domain.h" -#include "framework/specialization_registry.h" - -#include "engines/engine.hpp" - -namespace ntt { - - template - requires IsCompatibleWithEngine - void Engine::run() { - if constexpr (pgen_is_ok) { - init(); - - auto timers = timer::Timers { - { "FieldSolver", - "CurrentFiltering", "CurrentDeposit", - "ParticlePusher", "FieldBoundaries", - "ParticleBoundaries", "Communications", - "Injector", "Custom", - "PrtlClear", "Output", - "Checkpoint" }, - []() { - Kokkos::fence(); - }, - m_params.get("diagnostics.blocking_timers") - }; - const auto diag_interval = m_params.get( - "diagnostics.interval"); - - auto time_history = pbar::DurationHistory { 1000 }; - const auto clear_interval = m_params.template get( - "particles.clear_interval"); - - // main algorithm loop - while (step < max_steps) { - // run the engine-dependent algorithm step - m_metadomain.runOnLocalDomains([&timers, this](auto& dom) { - step_forward(timers, dom); - }); - // poststep (if defined) - if constexpr ( - traits::pgen::HasCustomPostStep>) { - timers.start("Custom"); - m_metadomain.runOnLocalDomains([&timers, this](auto& dom) { - m_pgen.CustomPostStep(step, time, dom); - }); - timers.stop("Custom"); - } - auto print_prtl_clear = (clear_interval > 0 and - step % clear_interval == 0 and step > 0); - - // advance time & step - time += dt; - ++step; - - auto print_output = false; - auto print_checkpoint = false; -#if defined(OUTPUT_ENABLED) - timers.start("Output"); - if constexpr ( - traits::pgen::HasCustomFieldOutput, M::Dim>) { - auto lambda_custom_field_output = [&](const std::string& name, - ndfield_t& buff, - index_t idx, - timestep_t step, - simtime_t time, - const Domain& dom) { - m_pgen.CustomFieldOutput(name, buff, idx, step, time, dom); - }; - print_output &= m_metadomain.Write(m_params, - step, - step - 1, - time, - time - dt, - lambda_custom_field_output); - } else { - print_output &= m_metadomain.Write(m_params, step, step - 1, time, time - dt); - } - if constexpr ( - traits::pgen::HasCustomStatOutput>) { - auto lambda_custom_stat = [&](const std::string& name, - timestep_t step, - simtime_t time, - const Domain& dom) -> real_t { - return m_pgen.CustomStat(name, step, time, dom); - }; - print_output &= m_metadomain.WriteStats(m_params, - step, - step - 1, - time, - time - dt, - lambda_custom_stat); - } else { - print_output &= m_metadomain.WriteStats(m_params, - step, - step - 1, - time, - time - dt); - } - timers.stop("Output"); - - timers.start("Checkpoint"); - print_checkpoint = m_metadomain.WriteCheckpoint(m_params, - step, - step - 1, - time, - time - dt); - timers.stop("Checkpoint"); -#endif - - // advance time_history - time_history.tick(); - // print timestep report - if (diag_interval > 0 and step % diag_interval == 0) { - diag::printDiagnostics( - step - 1, - max_steps, - time - dt, - dt, - timers, - time_history, - m_metadomain.l_ncells(), - m_metadomain.species_labels(), - m_metadomain.l_npart_perspec(), - m_metadomain.l_maxnpart_perspec(), - print_prtl_clear, - print_output, - print_checkpoint, - m_params.get("diagnostics.colored_stdout")); - } - timers.resetAll(); - } - } - } - -#define ENGINE_RUN(S, M, D) template void Engine>::run(); - - NTT_FOREACH_SPECIALIZATION(ENGINE_RUN) - -#undef ENGINE_RUN - -} // namespace ntt diff --git a/src/engines/engine_traits.h b/src/engines/engine_traits.h deleted file mode 100644 index d26b0add..00000000 --- a/src/engines/engine_traits.h +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef ENGINES_ENGINE_TRAITS_H -#define ENGINES_ENGINE_TRAITS_H - -#include "enums.h" - -#include "engines/grpic.hpp" -#include "engines/srpic.hpp" - -namespace ntt { - - template - struct EngineSelector; - - template <> - struct EngineSelector { - template - using type = SRPICEngine; - }; - - template <> - struct EngineSelector { - template - using type = GRPICEngine; - }; - -} // namespace ntt - -#endif // ENGINES_ENGINE_TRAITS_H diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index 275c70b5..8a5a5ebe 100644 --- a/src/engines/grpic.hpp +++ b/src/engines/grpic.hpp @@ -22,7 +22,7 @@ #include "utils/toml.h" #include "framework/domain/domain.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "engines/engine.hpp" #include "kernels/ampere_gr.hpp" @@ -71,8 +71,6 @@ namespace ntt { using base_t = Engine; using pgen_t = user::PGen; using domain_t = Domain; - // constexprs - using base_t::pgen_is_ok; // contents using base_t::m_metadomain; using base_t::m_params; @@ -1082,36 +1080,40 @@ namespace ntt { "algorithms.gr.pusher_eps"); const auto niter = m_params.template get( "algorithms.gr.pusher_niter"); - // clang-format off - if (species.pusher() == PrtlPusher::PHOTON) { - auto range_policy = Kokkos::RangePolicy( - 0, - species.npart()); + if (species.pusher() == ParticlePusher::PHOTON) { + auto range_policy = + Kokkos::RangePolicy( + 0, + species.npart()); - Kokkos::parallel_for( - "ParticlePusher", - range_policy, - kernel::gr::Pusher_kernel( - domain.fields.em, - domain.fields.em0, - species.i1, species.i2, species.i3, - species.i1_prev, species.i2_prev, species.i3_prev, - species.dx1, species.dx2, species.dx3, - species.dx1_prev, species.dx2_prev, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.phi, species.tag, - domain.mesh.metric, - coeff, dt, - domain.mesh.n_active(in::x1), - domain.mesh.n_active(in::x2), - domain.mesh.n_active(in::x3), - eps, niter, - domain.mesh.prtl_bc() - )); - } else if (species.pusher() == PrtlPusher::BORIS) { - auto range_policy = Kokkos::RangePolicy( - 0, - species.npart()); + // clang-format off + Kokkos::parallel_for( + "ParticlePusher", + range_policy, + kernel::gr::Pusher_kernel( + domain.fields.em, + domain.fields.em0, + species.i1, species.i2, species.i3, + species.i1_prev, species.i2_prev, species.i3_prev, + species.dx1, species.dx2, species.dx3, + species.dx1_prev, species.dx2_prev, species.dx3_prev, + species.ux1, species.ux2, species.ux3, + species.phi, species.tag, + domain.mesh.metric, + coeff, dt, + domain.mesh.n_active(in::x1), + domain.mesh.n_active(in::x2), + domain.mesh.n_active(in::x3), + eps, niter, + domain.mesh.prtl_bc() + )); + // clang-format on + } else if (species.pusher() == ParticlePusher::BORIS) { + auto range_policy = + Kokkos::RangePolicy( + 0, + species.npart()); + // clang-format off Kokkos::parallel_for( "ParticlePusher", range_policy, @@ -1132,12 +1134,12 @@ namespace ntt { eps, niter, domain.mesh.prtl_bc() )); - } else if (species.pusher() == PrtlPusher::NONE) { + // clang-format on + } else if (species.pusher() == ParticlePusher::NONE) { // do nothing } else { raise::Error("not implemented", HERE); } - // clang-format on } } }; diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp index fc95e6cf..006a3412 100644 --- a/src/engines/srpic.hpp +++ b/src/engines/srpic.hpp @@ -23,11 +23,13 @@ #include "utils/timer.h" #include "utils/toml.h" +#include "metrics/traits.h" + #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/spatial_dist.h" #include "framework/domain/domain.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "engines/engine.hpp" #include "kernels/ampere_mink.hpp" @@ -50,15 +52,13 @@ namespace ntt { template requires IsCompatibleWithEngine && - traits::metric::HasH_ij && traits::metric::HasConvert_i && - traits::metric::HasSqrtH_ij + metric::traits::HasH_ij && metric::traits::HasConvert_i && + metric::traits::HasSqrtH_ij class SRPICEngine : public Engine { using base_t = Engine; using pgen_t = user::PGen; using domain_t = Domain; - // constexprs - using base_t::pgen_is_ok; // contents using base_t::m_metadomain; using base_t::m_params; @@ -312,7 +312,7 @@ namespace ntt { } } for (auto& species : domain.species) { - if ((species.pusher() == PrtlPusher::NONE) or (species.npart() == 0)) { + if ((species.pusher() == ParticlePusher::NONE) or (species.npart() == 0)) { continue; } species.set_unsorted(); @@ -322,60 +322,69 @@ namespace ntt { species.label().c_str(), species.npart()), HERE); - const auto q_ovr_m = species.mass() > ZERO - ? species.charge() / species.mass() - : ZERO; - // coeff = q / m (dt / 2) omegaB0 - const auto coeff = q_ovr_m * HALF * dt * - m_params.template get("scales.omegaB0"); - PrtlPusher::type pusher; - if (species.pusher() == PrtlPusher::PHOTON) { - pusher = PrtlPusher::PHOTON; - } else if (species.pusher() == PrtlPusher::BORIS) { - pusher = PrtlPusher::BORIS; - } else if (species.pusher() == PrtlPusher::VAY) { - pusher = PrtlPusher::VAY; - } else { - raise::Fatal("Invalid particle pusher", HERE); + + kernel::sr::PusherParams pusher_params {}; + pusher_params.pusher_flags = species.pusher(); + pusher_params.radiative_drag_flags = species.radiative_drag_flags(); + pusher_params.mass = species.mass(); + pusher_params.charge = species.charge(); + pusher_params.time = time; + pusher_params.dt = dt; + pusher_params.omegaB0 = m_params.template get("scales.omegaB0"); + pusher_params.ni1 = domain.mesh.n_active(in::x1); + pusher_params.ni2 = domain.mesh.n_active(in::x2); + pusher_params.ni3 = domain.mesh.n_active(in::x3); + pusher_params.boundaries = domain.mesh.prtl_bc(); + + if (species.pusher() & ParticlePusher::GCA) { + pusher_params.gca_params.set( + "larmor_max", + m_params.template get("algorithms.gca.larmor_max")); + pusher_params.gca_params.set( + "e_ovr_b_max", + m_params.template get("algorithms.gca.e_ovr_b_max")); + } + + if (species.radiative_drag_flags() & RadiativeDrag::SYNCHROTRON) { + pusher_params.radiative_drag_params.set( + "synchrotron_gamma_rad", + m_params.template get( + "radiation.drag.synchrotron.gamma_rad")); + } + + if (species.radiative_drag_flags() & RadiativeDrag::COMPTON) { + pusher_params.radiative_drag_params.set( + "compton_gamma_rad", + m_params.template get("radiation.drag.compton.gamma_rad")); } - const auto cooling = species.cooling(); - - // coefficients to be forwarded to the dispatcher - // gca - const auto has_gca = species.use_gca(); - const auto gca_larmor_max = has_gca ? m_params.template get( - "algorithms.gca.larmor_max") - : ZERO; - const auto gca_eovrb_max = has_gca ? m_params.template get( - "algorithms.gca.e_ovr_b_max") - : ZERO; - // cooling - const auto has_synchrotron = (cooling == Cooling::SYNCHROTRON); - const auto has_compton = (cooling == Cooling::COMPTON); - const auto sync_grad = has_synchrotron - ? m_params.template get( - "algorithms.synchrotron.gamma_rad") - : ZERO; - const auto sync_coeff = has_synchrotron - ? (real_t)(0.1) * dt * - m_params.template get( - "scales.omegaB0") / - (SQR(sync_grad) * species.mass()) - : ZERO; - const auto comp_grad = has_compton ? m_params.template get( - "algorithms.compton.gamma_rad") - : ZERO; - const auto comp_coeff = has_compton ? (real_t)(0.1) * dt * - m_params.template get( - "scales.omegaB0") / - (SQR(comp_grad) * species.mass()) - : ZERO; + + kernel::sr::PusherArrays pusher_arrays {}; + pusher_arrays.sp = species.index(); + pusher_arrays.i1 = species.i1; + pusher_arrays.i2 = species.i2; + pusher_arrays.i3 = species.i3; + pusher_arrays.i1_prev = species.i1_prev; + pusher_arrays.i2_prev = species.i2_prev; + pusher_arrays.i3_prev = species.i3_prev; + pusher_arrays.dx1 = species.dx1; + pusher_arrays.dx2 = species.dx2; + pusher_arrays.dx3 = species.dx3; + pusher_arrays.dx1_prev = species.dx1_prev; + pusher_arrays.dx2_prev = species.dx2_prev; + pusher_arrays.dx3_prev = species.dx3_prev; + pusher_arrays.ux1 = species.ux1; + pusher_arrays.ux2 = species.ux2; + pusher_arrays.ux3 = species.ux3; + pusher_arrays.phi = species.phi; + pusher_arrays.tag = species.tag; + // toggle to indicate whether pgen defines the external force - bool has_extforce = false; - if constexpr (traits::pgen::HasExtForce) { + bool has_extforce = false; + if constexpr (arch::traits::pgen::HasExtForce) { has_extforce = true; // toggle to indicate whether the ext force applies to current species - if (traits::has_member::value) { + if ( + ::traits::has_member<::traits::species_t, decltype(pgen_t::ext_force)>::value) { has_extforce &= std::find(m_pgen.ext_force.species.begin(), m_pgen.ext_force.species.end(), species.index()) != @@ -383,69 +392,32 @@ namespace ntt { } } - kernel::sr::CoolingTags cooling_tags = 0; - if (cooling == Cooling::SYNCHROTRON) { - cooling_tags = kernel::sr::Cooling::Synchrotron; - } - if (cooling == Cooling::COMPTON) { - cooling_tags = kernel::sr::Cooling::Compton; - } - // clang-format off + pusher_params.ext_force = has_extforce; + if (not has_atmosphere and not has_extforce) { - Kokkos::parallel_for( - "ParticlePusher", - species.rangeActiveParticles(), - kernel::sr::Pusher_kernel( - pusher, has_gca, false, - cooling_tags, - domain.fields.em, - species.index(), - species.i1, species.i2, species.i3, - species.i1_prev, species.i2_prev, species.i3_prev, - species.dx1, species.dx2, species.dx3, - species.dx1_prev, species.dx2_prev, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.phi, species.tag, - domain.mesh.metric, - time, coeff, dt, - domain.mesh.n_active(in::x1), - domain.mesh.n_active(in::x2), - domain.mesh.n_active(in::x3), - domain.mesh.prtl_bc(), - gca_larmor_max, gca_eovrb_max, sync_coeff, comp_coeff - )); + Kokkos::parallel_for("ParticlePusher", + species.rangeActiveParticles(), + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric)); } else if (has_atmosphere and not has_extforce) { const auto force = kernel::sr::Force { - {gx1, gx2, gx3}, + { gx1, gx2, gx3 }, x_surf, ds - }; + }; Kokkos::parallel_for( "ParticlePusher", species.rangeActiveParticles(), - kernel::sr::Pusher_kernel( - pusher, has_gca, false, - cooling_tags, - domain.fields.em, - species.index(), - species.i1, species.i2, species.i3, - species.i1_prev, species.i2_prev, species.i3_prev, - species.dx1, species.dx2, species.dx3, - species.dx1_prev, species.dx2_prev, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.phi, species.tag, - domain.mesh.metric, - force, - time, coeff, dt, - domain.mesh.n_active(in::x1), - domain.mesh.n_active(in::x2), - domain.mesh.n_active(in::x3), - domain.mesh.prtl_bc(), - gca_larmor_max, gca_eovrb_max, sync_coeff, comp_coeff - )); + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric, + force)); } else if (not has_atmosphere and has_extforce) { - if constexpr (traits::pgen::HasExtForce) { + if constexpr (arch::traits::pgen::HasExtForce) { const auto force = kernel::sr::Force { m_pgen.ext_force @@ -453,63 +425,35 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", species.rangeActiveParticles(), - kernel::sr::Pusher_kernel( - pusher, has_gca, true, - cooling_tags, - domain.fields.em, - species.index(), - species.i1, species.i2, species.i3, - species.i1_prev, species.i2_prev, species.i3_prev, - species.dx1, species.dx2, species.dx3, - species.dx1_prev, species.dx2_prev, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.phi, species.tag, - domain.mesh.metric, - force, - time, coeff, dt, - domain.mesh.n_active(in::x1), - domain.mesh.n_active(in::x2), - domain.mesh.n_active(in::x3), - domain.mesh.prtl_bc(), - gca_larmor_max, gca_eovrb_max, sync_coeff, comp_coeff - )); + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric, + force)); } else { raise::Error("External force not implemented", HERE); } } else { // has_atmosphere and has_extforce - if constexpr (traits::pgen::HasExtForce) { + if constexpr (arch::traits::pgen::HasExtForce) { const auto force = kernel::sr::Force { - m_pgen.ext_force, {gx1, gx2, gx3}, x_surf, ds - }; + m_pgen.ext_force, + { gx1, gx2, gx3 }, + x_surf, + ds + }; Kokkos::parallel_for( "ParticlePusher", species.rangeActiveParticles(), - kernel::sr::Pusher_kernel( - pusher, has_gca, true, - cooling_tags, - domain.fields.em, - species.index(), - species.i1, species.i2, species.i3, - species.i1_prev, species.i2_prev, species.i3_prev, - species.dx1, species.dx2, species.dx3, - species.dx1_prev, species.dx2_prev, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.phi, species.tag, - domain.mesh.metric, - force, - time, coeff, dt, - domain.mesh.n_active(in::x1), - domain.mesh.n_active(in::x2), - domain.mesh.n_active(in::x3), - domain.mesh.prtl_bc(), - gca_larmor_max, gca_eovrb_max, sync_coeff, comp_coeff - )); + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric, + force)); } else { raise::Error("External force not implemented", HERE); - } + } } - // clang-format on } } @@ -544,10 +488,9 @@ namespace ntt { void CurrentsDeposit(domain_t& domain) { auto scatter_cur = Kokkos::Experimental::create_scatter_view( domain.fields.cur); - auto shape_order = m_params.template get("algorithms.deposit.order"); for (auto& species : domain.species) { - if ((species.pusher() == PrtlPusher::NONE) or (species.npart() == 0) or - cmp::AlmostZero_host(species.charge())) { + if ((species.pusher() == ParticlePusher::NONE) or + (species.npart() == 0) or cmp::AlmostZero_host(species.charge())) { continue; } logger::Checkpoint( @@ -573,7 +516,7 @@ namespace ntt { const auto V0 = m_params.template get("scales.V0"); const auto ppc0 = m_params.template get("particles.ppc0"); const auto coeff = -dt * q0 / (B0 * V0); - if constexpr (traits::pgen::HasExtCurrent) { + if constexpr (arch::traits::pgen::HasExtCurrent) { const std::vector xmin { domain.mesh.extent(in::x1).first, domain.mesh.extent(in::x2).first, domain.mesh.extent(in::x3).first }; @@ -722,7 +665,7 @@ namespace ntt { } if (dim == in::x1) { - if constexpr (traits::pgen::HasMatchFields) { + if constexpr (arch::traits::pgen::HasMatchFields) { auto match_fields = m_pgen.MatchFields(time); call_match_fields(domain.fields.em, domain.mesh.flds_bc(), @@ -733,7 +676,7 @@ namespace ntt { tags, range_min, range_max); - } else if constexpr (traits::pgen::HasMatchFieldsInX1) { + } else if constexpr (arch::traits::pgen::HasMatchFieldsInX1) { auto match_fields = m_pgen.MatchFieldsInX1(time); call_match_fields(domain.fields.em, domain.mesh.flds_bc(), @@ -747,7 +690,7 @@ namespace ntt { } } else if (dim == in::x2) { if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - if constexpr (traits::pgen::HasMatchFields) { + if constexpr (arch::traits::pgen::HasMatchFields) { auto match_fields = m_pgen.MatchFields(time); call_match_fields(domain.fields.em, domain.mesh.flds_bc(), @@ -758,7 +701,7 @@ namespace ntt { tags, range_min, range_max); - } else if constexpr (traits::pgen::HasMatchFieldsInX2) { + } else if constexpr (arch::traits::pgen::HasMatchFieldsInX2) { auto match_fields = m_pgen.MatchFieldsInX2(time); call_match_fields(domain.fields.em, domain.mesh.flds_bc(), @@ -775,7 +718,7 @@ namespace ntt { } } else if (dim == in::x3) { if constexpr (M::Dim == Dim::_3D) { - if constexpr (traits::pgen::HasMatchFields) { + if constexpr (arch::traits::pgen::HasMatchFields) { auto match_fields = m_pgen.MatchFields(time); call_match_fields(domain.fields.em, domain.mesh.flds_bc(), @@ -786,7 +729,7 @@ namespace ntt { tags, range_min, range_max); - } else if constexpr (traits::pgen::HasMatchFieldsInX3) { + } else if constexpr (arch::traits::pgen::HasMatchFieldsInX3) { auto match_fields = m_pgen.MatchFieldsInX3(time); call_match_fields(domain.fields.em, domain.mesh.flds_bc(), @@ -896,7 +839,7 @@ namespace ntt { if (tags & BC::B) { comps.push_back(normal_b_comp); } - if constexpr (traits::pgen::HasFixFieldsConst) { + if constexpr (arch::traits::pgen::HasFixFieldsConst) { for (const auto& comp : comps) { auto value = ZERO; bool shouldset = false; @@ -1068,7 +1011,7 @@ namespace ntt { /** * atmosphere field boundaries */ - if constexpr (traits::pgen::HasAtmFields) { + if constexpr (arch::traits::pgen::HasAtmFields) { const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction); const auto dd = static_cast(dim); boundaries_t box; diff --git a/src/engines/traits.h b/src/engines/traits.h new file mode 100644 index 00000000..f6466117 --- /dev/null +++ b/src/engines/traits.h @@ -0,0 +1,26 @@ +/** + * @file engine/traits.h + * @brief Defines a set of traits to check if an engine class satisfies certain conditions + * @implements + * - ntt::traits::engine::HasRun<> - checks if an engine has a run() method + * @namespaces: + * - ntt::traits::engine:: + */ +#ifndef ENGINES_TRAITS_H +#define ENGINES_TRAITS_H + +#include + +namespace ntt { + + namespace traits { + namespace engine { + template + concept HasRun = requires(E& engine) { + { engine.run() } -> std::same_as; + }; + } // namespace engine + } // namespace traits +} // namespace ntt + +#endif // ENGINES_TRAITS_H diff --git a/src/entity.cpp b/src/entity.cpp index 8d88e265..969a2637 100644 --- a/src/entity.cpp +++ b/src/entity.cpp @@ -3,19 +3,39 @@ #include "arch/traits.h" #include "utils/error.h" -#include "engines/engine_traits.h" +#include "archetypes/traits.h" #include "framework/simulation.h" #include "framework/specialization_registry.h" +#include "engines/grpic.hpp" +#include "engines/srpic.hpp" #include "pgen.hpp" - + #include +namespace ntt { + template + struct EngineSelector; + + template <> + struct EngineSelector { + template + using type = SRPICEngine; + }; + + template <> + struct EngineSelector { + template + using type = GRPICEngine; + }; +} // namespace ntt + template class M, Dimension D> static constexpr bool should_compile { - traits::check_compatibility::value(user::PGen>::engines) && - traits::check_compatibility::MetricType>::value(user::PGen>::metrics) && - traits::check_compatibility::value(user::PGen>::dimensions) + arch::traits::pgen::check_compatibility::value(user::PGen>::engines) && + arch::traits::pgen::check_compatibility::MetricType>::value( + user::PGen>::metrics) && + arch::traits::pgen::check_compatibility::value(user::PGen>::dimensions) }; template class M, Dimension D> @@ -26,7 +46,7 @@ void dispatch_engine(ntt::Simulation& sim) { sim.run::template type, M, D>(); } else { static_assert( - traits::always_false>::value, + ::traits::always_false>::value, "Unsupported engine"); } } diff --git a/src/framework/CMakeLists.txt b/src/framework/CMakeLists.txt index 01496787..907b04f2 100644 --- a/src/framework/CMakeLists.txt +++ b/src/framework/CMakeLists.txt @@ -4,7 +4,12 @@ # # @sources: # -# * parameters.cpp +# * parameters/parameters.cpp +# * parameters/particles.cpp +# * parameters/grid.cpp +# * parameters/output.cpp +# * parameters/algorithms.cpp +# * parameters/extra.cpp # * simulation.cpp # * domain/grid.cpp # * domain/metadomain.cpp @@ -36,8 +41,13 @@ set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}) set(SOURCES - ${SRC_DIR}/parameters.cpp ${SRC_DIR}/simulation.cpp + ${SRC_DIR}/parameters/parameters.cpp + ${SRC_DIR}/parameters/particles.cpp + ${SRC_DIR}/parameters/grid.cpp + ${SRC_DIR}/parameters/output.cpp + ${SRC_DIR}/parameters/algorithms.cpp + ${SRC_DIR}/parameters/extra.cpp ${SRC_DIR}/domain/grid.cpp ${SRC_DIR}/domain/metadomain.cpp ${SRC_DIR}/domain/communications.cpp diff --git a/src/framework/containers/particles.cpp b/src/framework/containers/particles.cpp index e9e51622..b1b77497 100644 --- a/src/framework/containers/particles.cpp +++ b/src/framework/containers/particles.cpp @@ -16,26 +16,26 @@ namespace ntt { template - Particles::Particles(spidx_t index, - const std::string& label, - float m, - float ch, - npart_t maxnpart, - const PrtlPusher& pusher, - bool use_tracking, - bool use_gca, - const Cooling& cooling, - unsigned short npld_r, - unsigned short npld_i) + Particles::Particles(spidx_t index, + const std::string& label, + float m, + float ch, + npart_t maxnpart, + ParticlePusherFlags particle_pusher_flags, + bool use_tracking, + RadiativeDragFlags radiative_drag_flags, + EmissionTypeFlag emission_policy_flag, + unsigned short npld_r, + unsigned short npld_i) : ParticleSpecies(index, label, m, ch, maxnpart, - pusher, + particle_pusher_flags, use_tracking, - use_gca, - cooling, + radiative_drag_flags, + emission_policy_flag, npld_r, npld_i) { diff --git a/src/framework/containers/particles.h b/src/framework/containers/particles.h index 144ca611..f265425e 100644 --- a/src/framework/containers/particles.h +++ b/src/framework/containers/particles.h @@ -83,24 +83,24 @@ namespace ntt { * @param m The mass of the species * @param ch The charge of the species * @param maxnpart The maximum number of allocated particles for the species - * @param pusher The pusher assigned for the species + * @param particle_pusher_flags The pusher(s) assigned for the species * @param use_tracking Use particle tracking for the species - * @param use_gca Use hybrid GCA pusher for the species - * @param cooling The cooling mechanism assigned for the species + * @param radiative_drag_flags The radiative drag mechanism(s) assigned for the species + * @param emission_policy_flag The emission policy assigned for the species * @param npld_r The number of real-valued payloads for the species * @param npld_i The number of integer-valued payloads for the species */ - Particles(spidx_t index, - const std::string& label, - float m, - float ch, - npart_t maxnpart, - const PrtlPusher& pusher, - bool use_gca, - bool use_tracking, - const Cooling& cooling, - unsigned short npld_r = 0, - unsigned short npld_i = 0); + Particles(spidx_t index, + const std::string& label, + float m, + float ch, + npart_t maxnpart, + ParticlePusherFlags particle_pusher_flags, + bool use_tracking, + RadiativeDragFlags radiative_drag_flags, + EmissionTypeFlag emission_policy_flag, + unsigned short npld_r, + unsigned short npld_i); /** * @brief Constructor for the particle container @@ -115,8 +115,8 @@ namespace ntt { spec.maxnpart(), spec.pusher(), spec.use_tracking(), - spec.use_gca(), - spec.cooling(), + spec.radiative_drag_flags(), + spec.emission_policy_flag(), spec.npld_r(), spec.npld_i()) {} diff --git a/src/framework/containers/species.h b/src/framework/containers/species.h index baf02487..5f374249 100644 --- a/src/framework/containers/species.h +++ b/src/framework/containers/species.h @@ -31,16 +31,16 @@ namespace ntt { npart_t m_maxnpart; // Pusher assigned for the species - const PrtlPusher m_pusher; + const ParticlePusherFlags m_particle_pusher_flags; // Use particle tracking for the species const bool m_use_tracking; - // Use byrid gca pusher for the species - const bool m_use_gca; + // Radiative drag mechanism(s) assigned for the species + const RadiativeDragFlags m_radiative_drag_flags; - // Cooling drag mechanism assigned for the species - const Cooling m_cooling; + // Emission policy assigned for the species + const EmissionTypeFlag m_emission_policy_flag; // Number of payloads for the species const unsigned short m_npld_r; @@ -53,10 +53,10 @@ namespace ntt { , m_mass { 0.0 } , m_charge { 0.0 } , m_maxnpart { 0 } - , m_pusher { PrtlPusher::INVALID } + , m_particle_pusher_flags { ParticlePusher::NONE } , m_use_tracking { false } - , m_use_gca { false } - , m_cooling { Cooling::INVALID } + , m_radiative_drag_flags { RadiativeDrag::NONE } + , m_emission_policy_flag { EmissionType::NONE } , m_npld_r { 0 } , m_npld_i { 0 } {} @@ -68,33 +68,33 @@ namespace ntt { * @param m The mass of the species. * @param ch The charge of the species. * @param maxnpart The maximum number of allocated particles for the species. - * @param pusher The pusher assigned for the species. + * @param particle_pusher_flags The pusher(s) assigned for the species. * @param use_tracking Use particle tracking for the species. - * @param use_gca Use hybrid GCA pusher for the species. - * @param cooling The cooling mechanism assigned for the species. + * @param radiative_drag_flags The radiative drag mechanism(s) assigned for the species. + * @param emission_policy_flag The emission policy assigned for the species. * @param npld_r The number of real-valued payloads for the species * @param npld_i The number of integer-valued payloads for the species */ - ParticleSpecies(spidx_t index, - const std::string& label, - float m, - float ch, - npart_t maxnpart, - const PrtlPusher& pusher, - bool use_tracking, - bool use_gca, - const Cooling& cooling, - unsigned short npld_r = 0, - unsigned short npld_i = 0) + ParticleSpecies(spidx_t index, + const std::string& label, + float m, + float ch, + npart_t maxnpart, + ParticlePusherFlags particle_pusher_flags, + bool use_tracking, + RadiativeDragFlags radiative_drag_flags, + EmissionTypeFlag emission_policy_flag, + unsigned short npld_r, + unsigned short npld_i) : m_index { index } , m_label { std::move(label) } , m_mass { m } , m_charge { ch } , m_maxnpart { maxnpart } - , m_pusher { pusher } + , m_particle_pusher_flags { particle_pusher_flags } , m_use_tracking { use_tracking } - , m_use_gca { use_gca } - , m_cooling { cooling } + , m_radiative_drag_flags { radiative_drag_flags } + , m_emission_policy_flag { emission_policy_flag } , m_npld_r { npld_r } , m_npld_i { npld_i } { if (use_tracking) { @@ -144,8 +144,8 @@ namespace ntt { } [[nodiscard]] - auto pusher() const -> PrtlPusher { - return m_pusher; + auto pusher() const -> ParticlePusherFlags { + return m_particle_pusher_flags; } [[nodiscard]] @@ -154,13 +154,13 @@ namespace ntt { } [[nodiscard]] - auto use_gca() const -> bool { - return m_use_gca; + auto radiative_drag_flags() const -> RadiativeDragFlags { + return m_radiative_drag_flags; } [[nodiscard]] - auto cooling() const -> Cooling { - return m_cooling; + auto emission_policy_flag() const -> EmissionTypeFlag { + return m_emission_policy_flag; } [[nodiscard]] diff --git a/src/framework/domain/checkpoint.cpp b/src/framework/domain/checkpoint.cpp index 6e10eda2..693cafc1 100644 --- a/src/framework/domain/checkpoint.cpp +++ b/src/framework/domain/checkpoint.cpp @@ -8,7 +8,7 @@ #include "utils/log.h" #include "framework/domain/metadomain.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "framework/specialization_registry.h" namespace ntt { diff --git a/src/framework/domain/domain.h b/src/framework/domain/domain.h index c4546110..a7da455d 100644 --- a/src/framework/domain/domain.h +++ b/src/framework/domain/domain.h @@ -43,10 +43,11 @@ #include "global.h" #include "arch/directions.h" -#include "arch/traits.h" #include "utils/formatting.h" #include "utils/numeric.h" +#include "metrics/traits.h" + #include "framework/containers/fields.h" #include "framework/containers/particles.h" #include "framework/containers/species.h" @@ -60,7 +61,7 @@ namespace ntt { template - requires traits::metric::HasD + requires metric::traits::HasD struct Domain { static constexpr Dimension D { M::Dim }; @@ -148,7 +149,8 @@ namespace ntt { } /* setters -------------------------------------------------------------- */ - auto set_neighbor_idx(const dir::direction_t& dir, unsigned int idx) -> void { + auto set_neighbor_idx(const dir::direction_t& dir, unsigned int idx) + -> void { m_neighbor_idx[dir] = idx; } @@ -166,8 +168,8 @@ namespace ntt { }; template - inline auto operator<<(std::ostream& os, - const Domain& domain) -> std::ostream& { + inline auto operator<<(std::ostream& os, const Domain& domain) + -> std::ostream& { os << "Domain #" << domain.index(); #if defined(MPI_ENABLED) os << " [MPI rank: " << domain.mpi_rank() << "]"; diff --git a/src/framework/domain/mesh.h b/src/framework/domain/mesh.h index 3cfd9117..b820eb89 100644 --- a/src/framework/domain/mesh.h +++ b/src/framework/domain/mesh.h @@ -17,11 +17,12 @@ #include "global.h" #include "arch/directions.h" -#include "arch/traits.h" #include "utils/comparators.h" #include "utils/error.h" #include "utils/numeric.h" +#include "metrics/traits.h" + #include "framework/domain/grid.h" #include @@ -32,7 +33,7 @@ namespace ntt { template - requires traits::metric::HasD && traits::metric::HasConvert_i + requires metric::traits::HasD && metric::traits::HasConvert_i struct Mesh : public Grid { static constexpr bool is_mesh { true }; static constexpr Dimension D { M::Dim }; @@ -132,9 +133,8 @@ namespace ntt { * @note indices are already shifted by N_GHOSTS (i.e. they start at N_GHOSTS not 0) */ [[nodiscard]] - auto ExtentToRange( - boundaries_t box, - boundaries_t incl_ghosts) const -> boundaries_t { + auto ExtentToRange(boundaries_t box, boundaries_t incl_ghosts) const + -> boundaries_t { raise::ErrorIf(box.size() != M::Dim, "Invalid box dimension", HERE); raise::ErrorIf(incl_ghosts.size() != M::Dim, "Invalid incl_ghosts dimension", diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h index bd6ec864..3319cafd 100644 --- a/src/framework/domain/metadomain.h +++ b/src/framework/domain/metadomain.h @@ -19,12 +19,13 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "arch/traits.h" + +#include "metrics/traits.h" #include "framework/containers/species.h" #include "framework/domain/domain.h" #include "framework/domain/mesh.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "output/stats.h" #if defined(MPI_ENABLED) @@ -48,9 +49,9 @@ namespace ntt { template - concept IsCompatibleWithMetadomain = traits::metric::HasD && - traits::metric::HasConvert && - traits::metric::HasTotVolume; + concept IsCompatibleWithMetadomain = metric::traits::HasD && + metric::traits::HasConvert && + metric::traits::HasTotVolume; template requires IsCompatibleWithMetadomain diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index b4ff2123..5cced2d9 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -9,7 +9,7 @@ #include "framework/containers/particles.h" #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "framework/specialization_registry.h" #include "kernels/divergences.hpp" diff --git a/src/framework/domain/stats.cpp b/src/framework/domain/stats.cpp index 83d65956..811b3511 100644 --- a/src/framework/domain/stats.cpp +++ b/src/framework/domain/stats.cpp @@ -9,7 +9,7 @@ #include "framework/containers/particles.h" #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "framework/specialization_registry.h" #include "kernels/reduced_stats.hpp" diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp deleted file mode 100644 index f7ef82fd..00000000 --- a/src/framework/parameters.cpp +++ /dev/null @@ -1,1121 +0,0 @@ -#include "framework/parameters.h" - -#include "defaults.h" -#include "enums.h" -#include "global.h" - -#include "utils/error.h" -#include "utils/formatting.h" -#include "utils/log.h" -#include "utils/numeric.h" -#include "utils/toml.h" - -#include "metrics/kerr_schild.h" -#include "metrics/kerr_schild_0.h" -#include "metrics/minkowski.h" -#include "metrics/qkerr_schild.h" -#include "metrics/qspherical.h" -#include "metrics/spherical.h" - -#include "framework/containers/species.h" - -#if defined(MPI_ENABLED) - #include -#endif - -#include -#include -#include -#include -#include -#include - -namespace ntt { - - template - auto get_dx0_V0( - const std::vector& resolution, - const boundaries_t& extent, - const std::map& params) -> std::pair { - const auto metric = M(resolution, extent, params); - const auto dx0 = metric.dxMin(); - coord_t x_corner { ZERO }; - for (auto d { 0u }; d < M::Dim; ++d) { - x_corner[d] = HALF; - } - const auto V0 = metric.sqrt_det_h(x_corner); - return { dx0, V0 }; - } - - /* - * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - * Parameters that must not be changed during the checkpoint restart - * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - */ - void SimulationParams::setImmutableParams(const toml::value& toml_data) { - /* [simulation] --------------------------------------------------------- */ - const auto engine_enum = SimEngine::pick( - fmt::toLower(toml::find(toml_data, "simulation", "engine")).c_str()); - set("simulation.engine", engine_enum); - - int default_ndomains = 1; -#if defined(MPI_ENABLED) - raise::ErrorIf(MPI_Comm_size(MPI_COMM_WORLD, &default_ndomains) != MPI_SUCCESS, - "MPI_Comm_size failed", - HERE); -#endif - const auto ndoms = toml::find_or(toml_data, - "simulation", - "domain", - "number", - default_ndomains); - set("simulation.domain.number", (unsigned int)ndoms); - - auto decomposition = toml::find_or>( - toml_data, - "simulation", - "domain", - "decomposition", - std::vector { -1, -1, -1 }); - promiseToDefine("simulation.domain.decomposition"); - - /* [grid] --------------------------------------------------------------- */ - const auto res = toml::find>(toml_data, - "grid", - "resolution"); - raise::ErrorIf(res.size() < 1 || res.size() > 3, - "invalid `grid.resolution`", - HERE); - set("grid.resolution", res); - const auto dim = static_cast(res.size()); - set("grid.dim", dim); - - if (decomposition.size() > dim) { - decomposition.erase(decomposition.begin() + (std::size_t)(dim), - decomposition.end()); - } - raise::ErrorIf(decomposition.size() != dim, - "invalid `simulation.domain.decomposition`", - HERE); - set("simulation.domain.decomposition", decomposition); - - auto extent = toml::find>>(toml_data, - "grid", - "extent"); - raise::ErrorIf(extent.size() < 1 || extent.size() > 3, - "invalid `grid.extent`", - HERE); - promiseToDefine("grid.extent"); - - /* [grid.metric] -------------------------------------------------------- */ - const auto metric_enum = Metric::pick( - fmt::toLower(toml::find(toml_data, "grid", "metric", "metric")) - .c_str()); - promiseToDefine("grid.metric.metric"); - std::string coord; - if (metric_enum == Metric::Minkowski) { - raise::ErrorIf(engine_enum != SimEngine::SRPIC, - "minkowski metric is only supported for SRPIC", - HERE); - coord = "cart"; - } else if (metric_enum == Metric::QKerr_Schild or - metric_enum == Metric::QSpherical) { - // quasi-spherical geometry - raise::ErrorIf(dim == Dim::_1D, - "not enough dimensions for qspherical geometry", - HERE); - raise::ErrorIf(dim == Dim::_3D, - "3D not implemented for qspherical geometry", - HERE); - coord = "qsph"; - set("grid.metric.qsph_r0", - toml::find_or(toml_data, "grid", "metric", "qsph_r0", defaults::qsph::r0)); - set("grid.metric.qsph_h", - toml::find_or(toml_data, "grid", "metric", "qsph_h", defaults::qsph::h)); - } else { - // spherical geometry - raise::ErrorIf(dim == Dim::_1D, - "not enough dimensions for spherical geometry", - HERE); - raise::ErrorIf(dim == Dim::_3D, - "3D not implemented for spherical geometry", - HERE); - coord = "sph"; - } - if ((engine_enum == SimEngine::GRPIC) && - (metric_enum != Metric::Kerr_Schild_0)) { - const auto ks_a = toml::find_or(toml_data, - "grid", - "metric", - "ks_a", - defaults::ks::a); - set("grid.metric.ks_a", ks_a); - set("grid.metric.ks_rh", ONE + math::sqrt(ONE - SQR(ks_a))); - } - const auto coord_enum = Coord::pick(coord.c_str()); - set("grid.metric.coord", coord_enum); - - /* [scales] ------------------------------------------------------------- */ - const auto larmor0 = toml::find(toml_data, "scales", "larmor0"); - const auto skindepth0 = toml::find(toml_data, "scales", "skindepth0"); - raise::ErrorIf(larmor0 <= ZERO || skindepth0 <= ZERO, - "larmor0 and skindepth0 must be positive", - HERE); - set("scales.larmor0", larmor0); - set("scales.skindepth0", skindepth0); - promiseToDefine("scales.dx0"); - promiseToDefine("scales.V0"); - promiseToDefine("scales.n0"); - promiseToDefine("scales.q0"); - set("scales.sigma0", SQR(skindepth0 / larmor0)); - set("scales.B0", ONE / larmor0); - set("scales.omegaB0", ONE / larmor0); - - /* [particles] ---------------------------------------------------------- */ - const auto ppc0 = toml::find(toml_data, "particles", "ppc0"); - set("particles.ppc0", ppc0); - raise::ErrorIf(ppc0 <= 0.0, "ppc0 must be positive", HERE); - set("particles.use_weights", - toml::find_or(toml_data, "particles", "use_weights", false)); - - /* [particles.species] -------------------------------------------------- */ - std::vector species; - const auto species_tab = toml::find_or(toml_data, - "particles", - "species", - toml::array {}); - set("particles.nspec", species_tab.size()); - - spidx_t idx = 1; - for (const auto& sp : species_tab) { - const auto label = toml::find_or(sp, - "label", - "s" + std::to_string(idx)); - const auto mass = toml::find(sp, "mass"); - const auto charge = toml::find(sp, "charge"); - raise::ErrorIf((charge != 0.0f) && (mass == 0.0f), - "mass of the charged species must be non-zero", - HERE); - const auto is_massless = (mass == 0.0f) && (charge == 0.0f); - const auto def_pusher = (is_massless ? defaults::ph_pusher - : defaults::em_pusher); - const auto maxnpart_real = toml::find(sp, "maxnpart"); - const auto maxnpart = static_cast(maxnpart_real); - auto pusher = toml::find_or(sp, "pusher", std::string(def_pusher)); - const auto npayloads_real = toml::find_or(sp, - "n_payloads_real", - static_cast(0)); - const auto use_tracking = toml::find_or(sp, "tracking", false); - auto npayloads_int = toml::find_or(sp, - "n_payloads_int", - static_cast(0)); - if (use_tracking) { -#if !defined(MPI_ENABLED) - npayloads_int += 1; -#else - npayloads_int += 2; -#endif - } - const auto cooling = toml::find_or(sp, "cooling", std::string("None")); - raise::ErrorIf((fmt::toLower(cooling) != "none") && is_massless, - "cooling is only applicable to massive particles", - HERE); - raise::ErrorIf((fmt::toLower(pusher) == "photon") && !is_massless, - "photon pusher is only applicable to massless particles", - HERE); - bool use_gca = false; - if (pusher.find(',') != std::string::npos) { - raise::ErrorIf(fmt::toLower(pusher.substr(pusher.find(',') + 1, - pusher.size())) != "gca", - "invalid pusher syntax", - HERE); - use_gca = true; - pusher = pusher.substr(0, pusher.find(',')); - } - const auto pusher_enum = PrtlPusher::pick(pusher.c_str()); - const auto cooling_enum = Cooling::pick(cooling.c_str()); - if (use_gca) { - raise::ErrorIf(engine_enum != SimEngine::SRPIC, - "GCA pushers are only supported for SRPIC", - HERE); - promiseToDefine("algorithms.gca.e_ovr_b_max"); - promiseToDefine("algorithms.gca.larmor_max"); - } - if (cooling_enum == Cooling::SYNCHROTRON) { - raise::ErrorIf(engine_enum != SimEngine::SRPIC, - "Synchrotron cooling is only supported for SRPIC", - HERE); - promiseToDefine("algorithms.synchrotron.gamma_rad"); - } - - if (cooling_enum == Cooling::COMPTON) { - raise::ErrorIf(engine_enum != SimEngine::SRPIC, - "Inverse Compton cooling is only supported for SRPIC", - HERE); - promiseToDefine("algorithms.compton.gamma_rad"); - } - - species.emplace_back(ParticleSpecies(idx, - label, - mass, - charge, - maxnpart, - pusher_enum, - use_tracking, - use_gca, - cooling_enum, - npayloads_real, - npayloads_int)); - idx += 1; - } - set("particles.species", species); - - /* inferred variables --------------------------------------------------- */ - // extent - if (extent.size() > dim) { - extent.erase(extent.begin() + (std::size_t)(dim), extent.end()); - } - raise::ErrorIf(extent[0].size() != 2, "invalid `grid.extent[0]`", HERE); - if (coord_enum != Coord::Cart) { - raise::ErrorIf(extent.size() > 1, - "invalid `grid.extent` for non-cartesian geometry", - HERE); - extent.push_back({ ZERO, constant::PI }); - if (dim == Dim::_3D) { - extent.push_back({ ZERO, TWO * constant::PI }); - } - } - raise::ErrorIf(extent.size() != dim, "invalid inferred `grid.extent`", HERE); - boundaries_t extent_pairwise; - for (auto d { 0u }; d < (dim_t)dim; ++d) { - raise::ErrorIf(extent[d].size() != 2, - fmt::format("invalid inferred `grid.extent[%d]`", d), - HERE); - extent_pairwise.push_back({ extent[d][0], extent[d][1] }); - } - set("grid.extent", extent_pairwise); - - // metric, dx0, V0, n0, q0 - { - boundaries_t ext; - for (const auto& e : extent) { - ext.push_back({ e[0], e[1] }); - } - std::map params; - if (coord_enum == Coord::Qsph) { - params["r0"] = get("grid.metric.qsph_r0"); - params["h"] = get("grid.metric.qsph_h"); - } - if ((engine_enum == SimEngine::GRPIC) && - (metric_enum != Metric::Kerr_Schild_0)) { - params["a"] = get("grid.metric.ks_a"); - } - set("grid.metric.params", params); - - std::pair dx0_V0; - if (metric_enum == Metric::Minkowski) { - if (dim == Dim::_1D) { - dx0_V0 = get_dx0_V0>(res, ext, params); - } else if (dim == Dim::_2D) { - dx0_V0 = get_dx0_V0>(res, ext, params); - } else { - dx0_V0 = get_dx0_V0>(res, ext, params); - } - } else if (metric_enum == Metric::Spherical) { - dx0_V0 = get_dx0_V0>(res, ext, params); - } else if (metric_enum == Metric::QSpherical) { - dx0_V0 = get_dx0_V0>(res, ext, params); - } else if (metric_enum == Metric::Kerr_Schild) { - dx0_V0 = get_dx0_V0>(res, ext, params); - } else if (metric_enum == Metric::Kerr_Schild_0) { - dx0_V0 = get_dx0_V0>(res, ext, params); - } else if (metric_enum == Metric::QKerr_Schild) { - dx0_V0 = get_dx0_V0>(res, ext, params); - } - auto [dx0, V0] = dx0_V0; - set("scales.dx0", dx0); - set("scales.V0", V0); - set("scales.n0", ppc0 / V0); - set("scales.q0", V0 / (ppc0 * SQR(skindepth0))); - - set("grid.metric.metric", metric_enum); - } - } - - /* - * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - * Parameters that may be changed during the checkpoint restart - * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - */ - void SimulationParams::setMutableParams(const toml::value& toml_data) { - const auto engine_enum = get("simulation.engine"); - const auto coord_enum = get("grid.metric.coord"); - const auto dim = get("grid.dim"); - const auto extent_pairwise = get>("grid.extent"); - - /* [simulation] --------------------------------------------------------- */ - set("simulation.name", - toml::find(toml_data, "simulation", "name")); - set("simulation.runtime", - toml::find(toml_data, "simulation", "runtime")); - - /* [grid.boundaraies] --------------------------------------------------- */ - auto flds_bc = toml::find>>( - toml_data, - "grid", - "boundaries", - "fields"); - { - raise::ErrorIf(flds_bc.size() < 1 || flds_bc.size() > 3, - "invalid `grid.boundaries.fields`", - HERE); - promiseToDefine("grid.boundaries.fields"); - auto atm_defined = false; - for (const auto& bcs : flds_bc) { - for (const auto& bc : bcs) { - if (fmt::toLower(bc) == "match") { - promiseToDefine("grid.boundaries.match.ds"); - } - if (fmt::toLower(bc) == "atmosphere") { - raise::ErrorIf(atm_defined, - "ATMOSPHERE is only allowed in one direction", - HERE); - atm_defined = true; - promiseToDefine("grid.boundaries.atmosphere.temperature"); - promiseToDefine("grid.boundaries.atmosphere.density"); - promiseToDefine("grid.boundaries.atmosphere.height"); - promiseToDefine("grid.boundaries.atmosphere.ds"); - promiseToDefine("grid.boundaries.atmosphere.species"); - promiseToDefine("grid.boundaries.atmosphere.g"); - } - } - } - } - - auto prtl_bc = toml::find>>( - toml_data, - "grid", - "boundaries", - "particles"); - { - raise::ErrorIf(prtl_bc.size() < 1 || prtl_bc.size() > 3, - "invalid `grid.boundaries.particles`", - HERE); - promiseToDefine("grid.boundaries.particles"); - auto atm_defined = false; - for (const auto& bcs : prtl_bc) { - for (const auto& bc : bcs) { - if (fmt::toLower(bc) == "absorb") { - promiseToDefine("grid.boundaries.absorb.ds"); - } - if (fmt::toLower(bc) == "atmosphere") { - raise::ErrorIf(atm_defined, - "ATMOSPHERE is only allowed in one direction", - HERE); - atm_defined = true; - promiseToDefine("grid.boundaries.atmosphere.temperature"); - promiseToDefine("grid.boundaries.atmosphere.density"); - promiseToDefine("grid.boundaries.atmosphere.height"); - promiseToDefine("grid.boundaries.atmosphere.ds"); - promiseToDefine("grid.boundaries.atmosphere.species"); - promiseToDefine("grid.boundaries.atmosphere.g"); - } - } - } - } - - /* [algorithms] --------------------------------------------------------- */ - set("algorithms.current_filters", - toml::find_or(toml_data, - "algorithms", - "current_filters", - defaults::current_filters)); - - /* [algorithms.deposit] ------------------------------------------------- */ - set("algorithms.deposit.enable", - toml::find_or(toml_data, "algorithms", "deposit", "enable", true)); - set("algorithms.deposit.order", - toml::find_or(toml_data, "algorithms", "deposit", "order", 1)); - - /* [algorithms.fieldsolver] --------------------------------------------- */ - set("algorithms.fieldsolver.enable", - toml::find_or(toml_data, "algorithms", "fieldsolver", "enable", true)); - - set("algorithms.fieldsolver.delta_x", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "delta_x", - defaults::fieldsolver::delta_x)); - set("algorithms.fieldsolver.delta_y", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "delta_y", - defaults::fieldsolver::delta_y)); - set("algorithms.fieldsolver.delta_z", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "delta_z", - defaults::fieldsolver::delta_z)); - set("algorithms.fieldsolver.beta_xy", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "beta_xy", - defaults::fieldsolver::beta_xy)); - set("algorithms.fieldsolver.beta_yx", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "beta_yx", - defaults::fieldsolver::beta_yx)); - set("algorithms.fieldsolver.beta_xz", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "beta_xz", - defaults::fieldsolver::beta_xz)); - set("algorithms.fieldsolver.beta_zx", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "beta_zx", - defaults::fieldsolver::beta_zx)); - set("algorithms.fieldsolver.beta_yz", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "beta_yz", - defaults::fieldsolver::beta_yz)); - set("algorithms.fieldsolver.beta_zy", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "beta_zy", - defaults::fieldsolver::beta_zy)); - /* [algorithms.timestep] ------------------------------------------------ */ - set("algorithms.timestep.CFL", - toml::find_or(toml_data, "algorithms", "timestep", "CFL", defaults::cfl)); - set("algorithms.timestep.dt", - get("algorithms.timestep.CFL") * get("scales.dx0")); - set("algorithms.timestep.correction", - toml::find_or(toml_data, - "algorithms", - "timestep", - "correction", - defaults::correction)); - - /* [algorithms.gr] ------------------------------------------------------ */ - if (engine_enum == SimEngine::GRPIC) { - set("algorithms.gr.pusher_eps", - toml::find_or(toml_data, - "algorithms", - "gr", - "pusher_eps", - defaults::gr::pusher_eps)); - set("algorithms.gr.pusher_niter", - toml::find_or(toml_data, - "algorithms", - "gr", - "pusher_niter", - defaults::gr::pusher_niter)); - } - /* [particles] ---------------------------------------------------------- */ - set("particles.clear_interval", - toml::find_or(toml_data, "particles", "clear_interval", defaults::clear_interval)); - const auto species_tab = toml::find_or(toml_data, - "particles", - "species", - toml::array {}); - std::vector species = get>( - "particles.species"); - raise::ErrorIf(species_tab.size() != species.size(), - "number of species changed after restart", - HERE); - - std::vector new_species; - - spidx_t idxM1 = 0; - for (const auto& sp : species_tab) { - const auto maxnpart_real = toml::find(sp, "maxnpart"); - const auto maxnpart = static_cast(maxnpart_real); - const auto particle_species = species[idxM1]; - new_species.emplace_back(particle_species.index(), - particle_species.label(), - particle_species.mass(), - particle_species.charge(), - maxnpart, - particle_species.pusher(), - particle_species.use_tracking(), - particle_species.use_gca(), - particle_species.cooling(), - particle_species.npld_r(), - particle_species.npld_i()); - idxM1++; - } - set("particles.species", new_species); - - /* [output] ------------------------------------------------------------- */ - // fields - set("output.format", - toml::find_or(toml_data, "output", "format", defaults::output::format)); - set("output.interval", - toml::find_or(toml_data, "output", "interval", defaults::output::interval)); - set("output.interval_time", - toml::find_or(toml_data, "output", "interval_time", -1.0)); - set("output.separate_files", - toml::find_or(toml_data, "output", "separate_files", true)); - - promiseToDefine("output.fields.enable"); - promiseToDefine("output.fields.interval"); - promiseToDefine("output.fields.interval_time"); - promiseToDefine("output.particles.enable"); - promiseToDefine("output.particles.interval"); - promiseToDefine("output.particles.interval_time"); - promiseToDefine("output.spectra.enable"); - promiseToDefine("output.spectra.interval"); - promiseToDefine("output.spectra.interval_time"); - promiseToDefine("output.stats.enable"); - promiseToDefine("output.stats.interval"); - promiseToDefine("output.stats.interval_time"); - - const auto flds_out = toml::find_or(toml_data, - "output", - "fields", - "quantities", - std::vector {}); - const auto custom_flds_out = toml::find_or(toml_data, - "output", - "fields", - "custom", - std::vector {}); - if (flds_out.size() == 0) { - raise::Warning("No fields output specified", HERE); - } - set("output.fields.quantities", flds_out); - set("output.fields.custom", custom_flds_out); - set("output.fields.mom_smooth", - toml::find_or(toml_data, - "output", - "fields", - "mom_smooth", - defaults::output::mom_smooth)); - std::vector field_dwn; - try { - auto field_dwn_ = toml::find>(toml_data, - "output", - "fields", - "downsampling"); - for (auto i = 0u; i < field_dwn_.size(); ++i) { - field_dwn.push_back(field_dwn_[i]); - } - } catch (...) { - try { - auto field_dwn_ = toml::find(toml_data, - "output", - "fields", - "downsampling"); - for (auto i = 0u; i < dim; ++i) { - field_dwn.push_back(field_dwn_); - } - } catch (...) { - for (auto i = 0u; i < dim; ++i) { - field_dwn.push_back(1u); - } - } - } - raise::ErrorIf(field_dwn.size() > 3, "invalid `output.fields.downsampling`", HERE); - if (field_dwn.size() > dim) { - field_dwn.erase(field_dwn.begin() + (std::size_t)(dim), field_dwn.end()); - } - for (const auto& dwn : field_dwn) { - raise::ErrorIf(dwn == 0, "downsampling factor must be nonzero", HERE); - } - set("output.fields.downsampling", field_dwn); - - // particles - auto all_specs = std::vector {}; - const auto nspec = get("particles.nspec"); - for (auto i = 0u; i < nspec; ++i) { - all_specs.push_back(static_cast(i + 1)); - } - const auto prtl_out = toml::find_or(toml_data, - "output", - "particles", - "species", - all_specs); - set("output.particles.species", prtl_out); - set("output.particles.stride", - toml::find_or(toml_data, - "output", - "particles", - "stride", - defaults::output::prtl_stride)); - - // spectra - set("output.spectra.e_min", - toml::find_or(toml_data, "output", "spectra", "e_min", defaults::output::spec_emin)); - set("output.spectra.e_max", - toml::find_or(toml_data, "output", "spectra", "e_max", defaults::output::spec_emax)); - set("output.spectra.log_bins", - toml::find_or(toml_data, - "output", - "spectra", - "log_bins", - defaults::output::spec_log)); - set("output.spectra.n_bins", - toml::find_or(toml_data, - "output", - "spectra", - "n_bins", - defaults::output::spec_nbins)); - - // stats - set("output.stats.quantities", - toml::find_or(toml_data, - "output", - "stats", - "quantities", - defaults::output::stats_quantities)); - set("output.stats.custom", - toml::find_or(toml_data, - "output", - "stats", - "custom", - std::vector {})); - - // intervals - for (const auto& type : { "fields", "particles", "spectra", "stats" }) { - const auto q_int = toml::find_or(toml_data, - "output", - std::string(type), - "interval", - 0); - const auto q_int_time = toml::find_or(toml_data, - "output", - std::string(type), - "interval_time", - -1.0); - set("output." + std::string(type) + ".enable", - toml::find_or(toml_data, "output", std::string(type), "enable", true)); - if ((q_int == 0) and (q_int_time == -1.0)) { - set("output." + std::string(type) + ".interval", - get("output.interval")); - set("output." + std::string(type) + ".interval_time", - get("output.interval_time")); - } else { - set("output." + std::string(type) + ".interval", q_int); - set("output." + std::string(type) + ".interval_time", q_int_time); - } - } - - /* [output.debug] ------------------------------------------------------- */ - set("output.debug.as_is", - toml::find_or(toml_data, "output", "debug", "as_is", false)); - const auto output_ghosts = toml::find_or(toml_data, - "output", - "debug", - "ghosts", - false); - set("output.debug.ghosts", output_ghosts); - if (output_ghosts) { - for (const auto& dwn : field_dwn) { - raise::ErrorIf( - dwn != 1, - "full resolution required when outputting with ghost cells", - HERE); - } - } - - /* [checkpoint] --------------------------------------------------------- */ - set("checkpoint.interval", - toml::find_or(toml_data, - "checkpoint", - "interval", - defaults::checkpoint::interval)); - set("checkpoint.interval_time", - toml::find_or(toml_data, "checkpoint", "interval_time", -1.0)); - set("checkpoint.keep", - toml::find_or(toml_data, "checkpoint", "keep", defaults::checkpoint::keep)); - auto walltime_str = toml::find_or(toml_data, - "checkpoint", - "walltime", - defaults::checkpoint::walltime); - if (walltime_str.empty()) { - walltime_str = defaults::checkpoint::walltime; - } - set("checkpoint.walltime", walltime_str); - - const auto checkpoint_write_path = toml::find_or( - toml_data, - "checkpoint", - "write_path", - fmt::format(defaults::checkpoint::write_path.c_str(), - get("simulation.name").c_str())); - set("checkpoint.write_path", checkpoint_write_path); - set("checkpoint.read_path", - toml::find_or(toml_data, "checkpoint", "read_path", checkpoint_write_path)); - - /* [diagnostics] -------------------------------------------------------- */ - set("diagnostics.interval", - toml::find_or(toml_data, "diagnostics", "interval", defaults::diag::interval)); - set("diagnostics.blocking_timers", - toml::find_or(toml_data, "diagnostics", "blocking_timers", false)); - set("diagnostics.colored_stdout", - toml::find_or(toml_data, "diagnostics", "colored_stdout", false)); - set("diagnostics.log_level", - toml::find_or(toml_data, "diagnostics", "log_level", defaults::diag::log_level)); - - /* inferred variables --------------------------------------------------- */ - // fields/particle boundaries - std::vector> flds_bc_enum; - std::vector> prtl_bc_enum; - if (coord_enum == Coord::Cart) { - raise::ErrorIf(flds_bc.size() != (std::size_t)dim, - "invalid `grid.boundaries.fields`", - HERE); - raise::ErrorIf(prtl_bc.size() != (std::size_t)dim, - "invalid `grid.boundaries.particles`", - HERE); - for (auto d { 0u }; d < (dim_t)dim; ++d) { - flds_bc_enum.push_back({}); - prtl_bc_enum.push_back({}); - const auto fbc = flds_bc[d]; - const auto pbc = prtl_bc[d]; - raise::ErrorIf(fbc.size() < 1 || fbc.size() > 2, - "invalid `grid.boundaries.fields`", - HERE); - raise::ErrorIf(pbc.size() < 1 || pbc.size() > 2, - "invalid `grid.boundaries.particles`", - HERE); - auto fbc_enum = FldsBC::pick(fmt::toLower(fbc[0]).c_str()); - auto pbc_enum = PrtlBC::pick(fmt::toLower(pbc[0]).c_str()); - if (fbc.size() == 1) { - raise::ErrorIf(fbc_enum != FldsBC::PERIODIC, - "invalid `grid.boundaries.fields`", - HERE); - flds_bc_enum.back().push_back(FldsBC(FldsBC::PERIODIC)); - flds_bc_enum.back().push_back(FldsBC(FldsBC::PERIODIC)); - } else { - raise::ErrorIf(fbc_enum == FldsBC::PERIODIC, - "invalid `grid.boundaries.fields`", - HERE); - flds_bc_enum.back().push_back(fbc_enum); - auto fbc_enum = FldsBC::pick(fmt::toLower(fbc[1]).c_str()); - raise::ErrorIf(fbc_enum == FldsBC::PERIODIC, - "invalid `grid.boundaries.fields`", - HERE); - flds_bc_enum.back().push_back(fbc_enum); - } - if (pbc.size() == 1) { - raise::ErrorIf(pbc_enum != PrtlBC::PERIODIC, - "invalid `grid.boundaries.particles`", - HERE); - prtl_bc_enum.back().push_back(PrtlBC(PrtlBC::PERIODIC)); - prtl_bc_enum.back().push_back(PrtlBC(PrtlBC::PERIODIC)); - } else { - raise::ErrorIf(pbc_enum == PrtlBC::PERIODIC, - "invalid `grid.boundaries.particles`", - HERE); - prtl_bc_enum.back().push_back(pbc_enum); - auto pbc_enum = PrtlBC::pick(fmt::toLower(pbc[1]).c_str()); - raise::ErrorIf(pbc_enum == PrtlBC::PERIODIC, - "invalid `grid.boundaries.particles`", - HERE); - prtl_bc_enum.back().push_back(pbc_enum); - } - } - } else { - raise::ErrorIf(flds_bc.size() > 1, "invalid `grid.boundaries.fields`", HERE); - raise::ErrorIf(prtl_bc.size() > 1, "invalid `grid.boundaries.particles`", HERE); - if (engine_enum == SimEngine::SRPIC) { - raise::ErrorIf(flds_bc[0].size() != 2, - "invalid `grid.boundaries.fields`", - HERE); - flds_bc_enum.push_back( - { FldsBC::pick(fmt::toLower(flds_bc[0][0]).c_str()), - FldsBC::pick(fmt::toLower(flds_bc[0][1]).c_str()) }); - flds_bc_enum.push_back({ FldsBC::AXIS, FldsBC::AXIS }); - if (dim == Dim::_3D) { - flds_bc_enum.push_back({ FldsBC::PERIODIC, FldsBC::PERIODIC }); - } - raise::ErrorIf(prtl_bc[0].size() != 2, - "invalid `grid.boundaries.particles`", - HERE); - prtl_bc_enum.push_back( - { PrtlBC::pick(fmt::toLower(prtl_bc[0][0]).c_str()), - PrtlBC::pick(fmt::toLower(prtl_bc[0][1]).c_str()) }); - prtl_bc_enum.push_back({ PrtlBC::AXIS, PrtlBC::AXIS }); - if (dim == Dim::_3D) { - prtl_bc_enum.push_back({ PrtlBC::PERIODIC, PrtlBC::PERIODIC }); - } - } else { - raise::ErrorIf(flds_bc[0].size() != 1, - "invalid `grid.boundaries.fields`", - HERE); - raise::ErrorIf(prtl_bc[0].size() != 1, - "invalid `grid.boundaries.particles`", - HERE); - flds_bc_enum.push_back( - { FldsBC::HORIZON, FldsBC::pick(fmt::toLower(flds_bc[0][0]).c_str()) }); - flds_bc_enum.push_back({ FldsBC::AXIS, FldsBC::AXIS }); - if (dim == Dim::_3D) { - flds_bc_enum.push_back({ FldsBC::PERIODIC, FldsBC::PERIODIC }); - } - prtl_bc_enum.push_back( - { PrtlBC::HORIZON, PrtlBC::pick(fmt::toLower(prtl_bc[0][0]).c_str()) }); - prtl_bc_enum.push_back({ PrtlBC::AXIS, PrtlBC::AXIS }); - if (dim == Dim::_3D) { - prtl_bc_enum.push_back({ PrtlBC::PERIODIC, PrtlBC::PERIODIC }); - } - } - } - - raise::ErrorIf(flds_bc_enum.size() != (std::size_t)dim, - "invalid inferred `grid.boundaries.fields`", - HERE); - raise::ErrorIf(prtl_bc_enum.size() != (std::size_t)dim, - "invalid inferred `grid.boundaries.particles`", - HERE); - boundaries_t flds_bc_pairwise; - boundaries_t prtl_bc_pairwise; - for (auto d { 0u }; d < (dim_t)dim; ++d) { - raise::ErrorIf( - flds_bc_enum[d].size() != 2, - fmt::format("invalid inferred `grid.boundaries.fields[%d]`", d), - HERE); - raise::ErrorIf( - prtl_bc_enum[d].size() != 2, - fmt::format("invalid inferred `grid.boundaries.particles[%d]`", d), - HERE); - flds_bc_pairwise.push_back({ flds_bc_enum[d][0], flds_bc_enum[d][1] }); - prtl_bc_pairwise.push_back({ prtl_bc_enum[d][0], prtl_bc_enum[d][1] }); - } - set("grid.boundaries.fields", flds_bc_pairwise); - set("grid.boundaries.particles", prtl_bc_pairwise); - - if (isPromised("grid.boundaries.match.ds")) { - if (coord_enum == Coord::Cart) { - auto min_extent = std::numeric_limits::max(); - for (const auto& e : extent_pairwise) { - min_extent = std::min(min_extent, e.second - e.first); - } - const auto default_ds = min_extent * defaults::bc::match::ds_frac; - boundaries_t ds_array; - try { - auto ds = toml::find(toml_data, "grid", "boundaries", "match", "ds"); - for (auto d = 0u; d < dim; ++d) { - ds_array.push_back({ ds, ds }); - } - } catch (...) { - try { - const auto ds = toml::find>>( - toml_data, - "grid", - "boundaries", - "match", - "ds"); - raise::ErrorIf(ds.size() != dim, - "invalid # in `grid.boundaries.match.ds`", - HERE); - for (auto d = 0u; d < dim; ++d) { - if (ds[d].size() == 1) { - ds_array.push_back({ ds[d][0], ds[d][0] }); - } else if (ds[d].size() == 2) { - ds_array.push_back({ ds[d][0], ds[d][1] }); - } else if (ds[d].size() == 0) { - ds_array.push_back({}); - } else { - raise::Error("invalid `grid.boundaries.match.ds`", HERE); - } - } - } catch (...) { - for (auto d = 0u; d < dim; ++d) { - ds_array.push_back({ default_ds, default_ds }); - } - } - } - set("grid.boundaries.match.ds", ds_array); - } else { - auto r_extent = extent_pairwise[0].second - extent_pairwise[0].first; - const auto ds = toml::find_or( - toml_data, - "grid", - "boundaries", - "match", - "ds", - r_extent * defaults::bc::match::ds_frac); - boundaries_t ds_array { - { ds, ds } - }; - set("grid.boundaries.match.ds", ds_array); - } - } - - if (isPromised("grid.boundaries.absorb.ds")) { - if (coord_enum == Coord::Cart) { - auto min_extent = std::numeric_limits::max(); - for (const auto& e : extent_pairwise) { - min_extent = std::min(min_extent, e.second - e.first); - } - set("grid.boundaries.absorb.ds", - toml::find_or(toml_data, - "grid", - "boundaries", - "absorb", - "ds", - min_extent * defaults::bc::absorb::ds_frac)); - } else { - auto r_extent = extent_pairwise[0].second - extent_pairwise[0].first; - set("grid.boundaries.absorb.ds", - toml::find_or(toml_data, - "grid", - "boundaries", - "absorb", - "ds", - r_extent * defaults::bc::absorb::ds_frac)); - } - } - - if (isPromised("grid.boundaries.atmosphere.temperature")) { - const auto atm_T = toml::find(toml_data, - "grid", - "boundaries", - "atmosphere", - "temperature"); - const auto atm_h = toml::find(toml_data, - "grid", - "boundaries", - "atmosphere", - "height"); - set("grid.boundaries.atmosphere.temperature", atm_T); - set("grid.boundaries.atmosphere.density", - toml::find(toml_data, "grid", "boundaries", "atmosphere", "density")); - set("grid.boundaries.atmosphere.ds", - toml::find_or(toml_data, "grid", "boundaries", "atmosphere", "ds", ZERO)); - set("grid.boundaries.atmosphere.height", atm_h); - set("grid.boundaries.atmosphere.g", atm_T / atm_h); - const auto atm_species = toml::find>( - toml_data, - "grid", - "boundaries", - "atmosphere", - "species"); - set("grid.boundaries.atmosphere.species", atm_species); - } - - // gca - if (isPromised("algorithms.gca.e_ovr_b_max")) { - set("algorithms.gca.e_ovr_b_max", - toml::find_or(toml_data, - "algorithms", - "gca", - "e_ovr_b_max", - defaults::gca::EovrB_max)); - set("algorithms.gca.larmor_max", - toml::find_or(toml_data, "algorithms", "gca", "larmor_max", ZERO)); - } - - // cooling - if (isPromised("algorithms.synchrotron.gamma_rad")) { - set("algorithms.synchrotron.gamma_rad", - toml::find_or(toml_data, - "algorithms", - "synchrotron", - "gamma_rad", - defaults::synchrotron::gamma_rad)); - } - if (isPromised("algorithms.compton.gamma_rad")) { - set("algorithms.compton.gamma_rad", - toml::find_or(toml_data, - "algorithms", - "compton", - "gamma_rad", - defaults::compton::gamma_rad)); - } - - // @TODO: disabling stats for non-Cartesian - if (coord_enum != Coord::Cart) { - set("output.stats.enable", false); - } - } - - void SimulationParams::setSetupParams(const toml::value& toml_data) { - /* [setup] -------------------------------------------------------------- */ - const auto setup = toml::find_or(toml_data, "setup", toml::table {}); - for (const auto& [key, val] : setup) { - if (val.is_boolean()) { - set("setup." + key, (bool)(val.as_boolean())); - } else if (val.is_integer()) { - set("setup." + key, (int)(val.as_integer())); - } else if (val.is_floating()) { - set("setup." + key, (real_t)(val.as_floating())); - } else if (val.is_string()) { - set("setup." + key, (std::string)(val.as_string())); - } else if (val.is_array()) { - const auto val_arr = val.as_array(); - if (val_arr.size() == 0) { - continue; - } else { - if (val_arr[0].is_integer()) { - std::vector vec; - for (const auto& v : val_arr) { - vec.push_back(v.as_integer()); - } - set("setup." + key, vec); - } else if (val_arr[0].is_floating()) { - std::vector vec; - for (const auto& v : val_arr) { - vec.push_back(v.as_floating()); - } - set("setup." + key, vec); - } else if (val_arr[0].is_boolean()) { - std::vector vec; - for (const auto& v : val_arr) { - vec.push_back(v.as_boolean()); - } - set("setup." + key, vec); - } else if (val_arr[0].is_string()) { - std::vector vec; - for (const auto& v : val_arr) { - vec.push_back(v.as_string()); - } - set("setup." + key, vec); - } else if (val_arr[0].is_array()) { - raise::Error("only 1D arrays allowed in [setup]", HERE); - } else { - raise::Error("invalid setup variable type", HERE); - } - } - } - } - } - - void SimulationParams::setCheckpointParams(bool is_resuming, - timestep_t start_step, - simtime_t start_time) { - set("checkpoint.is_resuming", is_resuming); - set("checkpoint.start_step", start_step); - set("checkpoint.start_time", start_time); - } - - void SimulationParams::checkPromises() const { - raise::ErrorIf(!promisesFulfilled(), - "Have not defined all the necessary variables", - HERE); - } - - void SimulationParams::saveTOML(const std::string& path, simtime_t time) const { - CallOnce([&]() { - std::ofstream metadata; - metadata.open(path); - metadata << "[metadata]\n" - << " time = " << time << "\n\n" - << data() << std::endl; - metadata.close(); - }); - } - -} // namespace ntt diff --git a/src/framework/parameters/algorithms.cpp b/src/framework/parameters/algorithms.cpp new file mode 100644 index 00000000..ce20136f --- /dev/null +++ b/src/framework/parameters/algorithms.cpp @@ -0,0 +1,157 @@ +#include "framework/parameters/algorithms.h" + +#include "defaults.h" +#include "global.h" + +#include "utils/numeric.h" +#include "utils/toml.h" + +#include "framework/parameters/parameters.h" + +namespace ntt { + namespace params { + + void Algorithms::read(real_t dx0, + const std::map& extra, + const toml::value& toml_data) { + CFL = toml::find_or(toml_data, "algorithms", "timestep", "CFL", defaults::cfl); + dt = CFL * dx0; + dt_correction_factor = toml::find_or(toml_data, + "algorithms", + "timestep", + "correction", + defaults::correction); + + number_of_current_filters = toml::find_or(toml_data, + "algorithms", + "current_filters", + defaults::current_filters); + + deposit_enable = toml::find_or(toml_data, "algorithms", "deposit", "enable", true); + deposit_order = static_cast(SHAPE_ORDER); + + fieldsolver_enable = toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "enable", + true); + + fieldsolver_stencil_coeffs["delta_x"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "delta_x", + defaults::fieldsolver::delta_x); + + fieldsolver_stencil_coeffs["delta_y"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "delta_y", + defaults::fieldsolver::delta_y); + + fieldsolver_stencil_coeffs["delta_z"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "delta_z", + defaults::fieldsolver::delta_z); + + fieldsolver_stencil_coeffs["beta_xy"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "beta_xy", + defaults::fieldsolver::beta_xy); + + fieldsolver_stencil_coeffs["beta_yx"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "beta_yx", + defaults::fieldsolver::beta_yx); + + fieldsolver_stencil_coeffs["beta_xz"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "beta_xz", + defaults::fieldsolver::beta_xz); + + fieldsolver_stencil_coeffs["beta_zx"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "beta_zx", + defaults::fieldsolver::beta_zx); + + fieldsolver_stencil_coeffs["beta_yz"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "beta_yz", + defaults::fieldsolver::beta_yz); + + fieldsolver_stencil_coeffs["beta_zy"] = toml::find_or( + toml_data, + "algorithms", + "fieldsolver", + "beta_zy", + defaults::fieldsolver::beta_zy); + + if (extra.at("gr")) { + gr_pusher_eps = toml::find_or(toml_data, + "algorithms", + "gr", + "pusher_eps", + defaults::gr::pusher_eps); + gr_pusher_niter = toml::find_or(toml_data, + "algorithms", + "gr", + "pusher_niter", + defaults::gr::pusher_niter); + } + + if (extra.at("gca")) { + gca_e_ovr_b_max = toml::find_or(toml_data, + "algorithms", + "gca", + "e_ovr_b_max", + defaults::gca::EovrB_max); + gca_larmor_max = toml::find_or(toml_data, + "algorithms", + "gca", + "larmor_max", + ZERO); + } + } + + void Algorithms::setParams(const std::map& extra, + SimulationParams* params) const { + params->set("algorithms.timestep.CFL", CFL); + params->set("algorithms.timestep.dt", dt); + params->set("algorithms.timestep.correction", dt_correction_factor); + + params->set("algorithms.current_filters", number_of_current_filters); + + params->set("algorithms.deposit.enable", deposit_enable); + params->set("algorithms.deposit.order", deposit_order); + + params->set("algorithms.fieldsolver.enable", fieldsolver_enable); + for (const auto& [key, value] : fieldsolver_stencil_coeffs) { + params->set("algorithms.fieldsolver." + key, value); + } + + if (extra.at("gr")) { + params->set("algorithms.gr.pusher_eps", gr_pusher_eps); + params->set("algorithms.gr.pusher_niter", gr_pusher_niter); + } + + if (extra.at("gca")) { + params->set("algorithms.gca.e_ovr_b_max", gca_e_ovr_b_max); + params->set("algorithms.gca.larmor_max", gca_larmor_max); + } + } + + } // namespace params +} // namespace ntt diff --git a/src/framework/parameters/algorithms.h b/src/framework/parameters/algorithms.h new file mode 100644 index 00000000..18bd2a66 --- /dev/null +++ b/src/framework/parameters/algorithms.h @@ -0,0 +1,56 @@ +/** + * @file framework/parameters/algorithms.h + * @brief Auxiliary functions for reading in algorithms parameters + * @implements + * - ntt::params::Algorithms + * @cpp: + * - algorithms.cpp + * @namespaces: + * - ntt::params:: + */ + +#ifndef FRAMEWORK_PARAMETERS_ALGORITHMS_H +#define FRAMEWORK_PARAMETERS_ALGORITHMS_H + +#include "global.h" + +#include "utils/toml.h" + +#include "framework/parameters/parameters.h" + +#include +#include + +namespace ntt { + namespace params { + + struct Algorithms { + real_t CFL; + real_t dt; + real_t dt_correction_factor; + + unsigned short number_of_current_filters; + + bool deposit_enable; + unsigned short deposit_order; + + bool fieldsolver_enable; + std::map fieldsolver_stencil_coeffs; + + real_t gr_pusher_eps; + unsigned short gr_pusher_niter; + + real_t gca_e_ovr_b_max; + real_t gca_larmor_max; + + real_t synchrotron_gamma_rad; + real_t compton_gamma_rad; + + void read(real_t, const std::map&, const toml::value&); + void setParams(const std::map&, SimulationParams*) const; + }; + + } // namespace params +} // namespace ntt + +#endif // FRAMEWORK_PARAMETERS_ALGORITHMS_H diff --git a/src/framework/parameters/extra.cpp b/src/framework/parameters/extra.cpp new file mode 100644 index 00000000..9a729d18 --- /dev/null +++ b/src/framework/parameters/extra.cpp @@ -0,0 +1,115 @@ +#include "framework/parameters/extra.h" + +#include "defaults.h" + +#include "utils/numeric.h" + +namespace ntt { + namespace params { + + void Extra::read(const std::map& extra, + const toml::value& toml_data) { + if (extra.at("synchrotron_drag")) { + synchrotron_gamma_rad = toml::find_or(toml_data, + "radiation", + "drag", + "synchrotron", + "gamma_rad", + defaults::synchrotron::gamma_rad); + } + + if (extra.at("compton_drag")) { + compton_gamma_rad = toml::find_or(toml_data, + "radiation", + "drag", + "compton", + "gamma_rad", + defaults::compton::gamma_rad); + } + + if (extra.at("synchrotron_emission")) { + synchrotron_gamma_qed = toml::find_or(toml_data, + "radiation", + "emission", + "synchrotron", + "gamma_qed", + defaults::synchrotron::gamma_qed); + synchrotron_energy_min = toml::find_or(toml_data, + "radiation", + "emission", + "synchrotron", + "photon_energy_min", + defaults::synchrotron::energy_min); + synchrotron_photon_weight = toml::find_or(toml_data, + "radiation", + "emission", + "synchrotron", + "photon_weight", + ONE); + synchrotron_photon_species = toml::find(toml_data, + "radiation", + "emission", + "synchrotron", + "photon_species"); + } + + if (extra.at("compton_emission")) { + compton_gamma_qed = toml::find_or(toml_data, + "radiation", + "emission", + "compton", + "gamma_qed", + defaults::compton::gamma_qed); + compton_energy_min = toml::find_or(toml_data, + "radiation", + "emission", + "compton", + "photon_energy_min", + defaults::compton::energy_min); + compton_photon_weight = toml::find_or(toml_data, + "radiation", + "emission", + "compton", + "photon_weight", + ONE); + compton_photon_species = toml::find(toml_data, + "radiation", + "emission", + "compton", + "photon_species"); + } + } + + void Extra::setParams(const std::map& extra, + SimulationParams* params) const { + if (extra.at("synchrotron_drag")) { + params->set("radiation.drag.synchrotron.gamma_rad", synchrotron_gamma_rad); + } + + if (extra.at("compton_drag")) { + params->set("radiation.drag.compton.gamma_rad", compton_gamma_rad); + } + + if (extra.at("synchrotron_emission")) { + params->set("radiation.emission.synchrotron.gamma_qed", + synchrotron_gamma_qed); + params->set("radiation.emission.synchrotron.photon_energy_min", + synchrotron_energy_min); + params->set("radiation.emission.synchrotron.photon_weight", + synchrotron_photon_weight); + params->set("radiation.emission.synchrotron.photon_species", + synchrotron_photon_species); + } + + if (extra.at("compton_emission")) { + params->set("radiation.emission.compton.gamma_qed", compton_gamma_qed); + params->set("radiation.emission.compton.photon_energy_min", + compton_energy_min); + params->set("radiation.emission.compton.photon_weight", + compton_photon_weight); + params->set("radiation.emission.compton.photon_species", + compton_photon_species); + } + } + } // namespace params +} // namespace ntt diff --git a/src/framework/parameters/extra.h b/src/framework/parameters/extra.h new file mode 100644 index 00000000..76fe9823 --- /dev/null +++ b/src/framework/parameters/extra.h @@ -0,0 +1,50 @@ +/** + * @file framework/parameters/algorithms.h + * @brief Auxiliary functions for reading in extra physics parameters + * @implements + * - ntt::params::Extra + * @cpp: + * - extra.cpp + * @namespaces: + * - ntt::params:: + */ + +#ifndef FRAMEWORK_PARAMETERS_EXTRA_H +#define FRAMEWORK_PARAMETERS_EXTRA_H + +#include "global.h" + +#include "utils/toml.h" + +#include "framework/parameters/parameters.h" + +#include +#include + +namespace ntt { + namespace params { + + struct Extra { + // radiative drag parameters + real_t synchrotron_gamma_rad; + real_t compton_gamma_rad; + + // emission parameters + real_t synchrotron_energy_min; + real_t synchrotron_gamma_qed; + real_t synchrotron_photon_weight; + spidx_t synchrotron_photon_species; + + real_t compton_energy_min; + real_t compton_gamma_qed; + real_t compton_photon_weight; + spidx_t compton_photon_species; + + void read(const std::map&, const toml::value&); + void setParams(const std::map&, SimulationParams*) const; + }; + + } // namespace params +} // namespace ntt + +#endif // FRAMEWORK_PARAMETERS_EXTRA_H diff --git a/src/framework/parameters/grid.cpp b/src/framework/parameters/grid.cpp new file mode 100644 index 00000000..705bb036 --- /dev/null +++ b/src/framework/parameters/grid.cpp @@ -0,0 +1,551 @@ +#include "framework/parameters/grid.h" + +#include "defaults.h" +#include "global.h" + +#include "utils/error.h" +#include "utils/formatting.h" +#include "utils/numeric.h" +#include "utils/toml.h" + +#include "metrics/kerr_schild.h" +#include "metrics/kerr_schild_0.h" +#include "metrics/minkowski.h" +#include "metrics/qkerr_schild.h" +#include "metrics/qspherical.h" +#include "metrics/spherical.h" + +#include "framework/parameters/parameters.h" + +#include +#include +#include +#include + +namespace ntt { + namespace params { + + template + auto get_dx0_V0( + const std::vector& resolution, + const boundaries_t& extent, + const std::map& params) -> std::pair { + const auto metric = M(resolution, extent, params); + const auto dx0 = metric.dxMin(); + coord_t x_corner { ZERO }; + for (auto d { 0u }; d < M::Dim; ++d) { + x_corner[d] = HALF; + } + const auto V0 = metric.sqrt_det_h(x_corner); + return { dx0, V0 }; + } + + auto GetBoundaryConditions(SimulationParams* params, + const SimEngine& engine_enum, + Dimension dim, + const Coord& coord_enum, + const toml::value& toml_data) + -> std::tuple, boundaries_t> { + auto flds_bc = toml::find>>( + toml_data, + "grid", + "boundaries", + "fields"); + { + raise::ErrorIf(flds_bc.size() < 1 || flds_bc.size() > 3, + "invalid `grid.boundaries.fields`", + HERE); + params->promiseToDefine("grid.boundaries.fields"); + auto atm_defined = false; + for (const auto& bcs : flds_bc) { + for (const auto& bc : bcs) { + if (fmt::toLower(bc) == "match") { + params->promiseToDefine("grid.boundaries.match.ds"); + } + if (fmt::toLower(bc) == "atmosphere") { + raise::ErrorIf(atm_defined, + "ATMOSPHERE is only allowed in one direction", + HERE); + atm_defined = true; + params->promiseToDefine("grid.boundaries.atmosphere.temperature"); + params->promiseToDefine("grid.boundaries.atmosphere.density"); + params->promiseToDefine("grid.boundaries.atmosphere.height"); + params->promiseToDefine("grid.boundaries.atmosphere.ds"); + params->promiseToDefine("grid.boundaries.atmosphere.species"); + params->promiseToDefine("grid.boundaries.atmosphere.g"); + } + } + } + } + + auto prtl_bc = toml::find>>( + toml_data, + "grid", + "boundaries", + "particles"); + { + raise::ErrorIf(prtl_bc.size() < 1 || prtl_bc.size() > 3, + "invalid `grid.boundaries.particles`", + HERE); + params->promiseToDefine("grid.boundaries.particles"); + auto atm_defined = false; + for (const auto& bcs : prtl_bc) { + for (const auto& bc : bcs) { + if (fmt::toLower(bc) == "absorb") { + params->promiseToDefine("grid.boundaries.absorb.ds"); + } + if (fmt::toLower(bc) == "atmosphere") { + raise::ErrorIf(atm_defined, + "ATMOSPHERE is only allowed in one direction", + HERE); + atm_defined = true; + params->promiseToDefine("grid.boundaries.atmosphere.temperature"); + params->promiseToDefine("grid.boundaries.atmosphere.density"); + params->promiseToDefine("grid.boundaries.atmosphere.height"); + params->promiseToDefine("grid.boundaries.atmosphere.ds"); + params->promiseToDefine("grid.boundaries.atmosphere.species"); + params->promiseToDefine("grid.boundaries.atmosphere.g"); + } + } + } + } + std::vector> flds_bc_enum; + std::vector> prtl_bc_enum; + if (coord_enum == Coord::Cart) { + raise::ErrorIf(flds_bc.size() != (std::size_t)dim, + "invalid `grid.boundaries.fields`", + HERE); + raise::ErrorIf(prtl_bc.size() != (std::size_t)dim, + "invalid `grid.boundaries.particles`", + HERE); + for (auto d { 0u }; d < (dim_t)dim; ++d) { + flds_bc_enum.push_back({}); + prtl_bc_enum.push_back({}); + const auto fbc = flds_bc[d]; + const auto pbc = prtl_bc[d]; + raise::ErrorIf(fbc.size() < 1 || fbc.size() > 2, + "invalid `grid.boundaries.fields`", + HERE); + raise::ErrorIf(pbc.size() < 1 || pbc.size() > 2, + "invalid `grid.boundaries.particles`", + HERE); + auto fbc_enum = FldsBC::pick(fmt::toLower(fbc[0]).c_str()); + auto pbc_enum = PrtlBC::pick(fmt::toLower(pbc[0]).c_str()); + if (fbc.size() == 1) { + raise::ErrorIf(fbc_enum != FldsBC::PERIODIC, + "invalid `grid.boundaries.fields`", + HERE); + flds_bc_enum.back().push_back(FldsBC(FldsBC::PERIODIC)); + flds_bc_enum.back().push_back(FldsBC(FldsBC::PERIODIC)); + } else { + raise::ErrorIf(fbc_enum == FldsBC::PERIODIC, + "invalid `grid.boundaries.fields`", + HERE); + flds_bc_enum.back().push_back(fbc_enum); + auto fbc_enum = FldsBC::pick(fmt::toLower(fbc[1]).c_str()); + raise::ErrorIf(fbc_enum == FldsBC::PERIODIC, + "invalid `grid.boundaries.fields`", + HERE); + flds_bc_enum.back().push_back(fbc_enum); + } + if (pbc.size() == 1) { + raise::ErrorIf(pbc_enum != PrtlBC::PERIODIC, + "invalid `grid.boundaries.particles`", + HERE); + prtl_bc_enum.back().push_back(PrtlBC(PrtlBC::PERIODIC)); + prtl_bc_enum.back().push_back(PrtlBC(PrtlBC::PERIODIC)); + } else { + raise::ErrorIf(pbc_enum == PrtlBC::PERIODIC, + "invalid `grid.boundaries.particles`", + HERE); + prtl_bc_enum.back().push_back(pbc_enum); + auto pbc_enum = PrtlBC::pick(fmt::toLower(pbc[1]).c_str()); + raise::ErrorIf(pbc_enum == PrtlBC::PERIODIC, + "invalid `grid.boundaries.particles`", + HERE); + prtl_bc_enum.back().push_back(pbc_enum); + } + } + } else { + raise::ErrorIf(flds_bc.size() > 1, "invalid `grid.boundaries.fields`", HERE); + raise::ErrorIf(prtl_bc.size() > 1, + "invalid `grid.boundaries.particles`", + HERE); + if (engine_enum == SimEngine::SRPIC) { + raise::ErrorIf(flds_bc[0].size() != 2, + "invalid `grid.boundaries.fields`", + HERE); + flds_bc_enum.push_back( + { FldsBC::pick(fmt::toLower(flds_bc[0][0]).c_str()), + FldsBC::pick(fmt::toLower(flds_bc[0][1]).c_str()) }); + flds_bc_enum.push_back({ FldsBC::AXIS, FldsBC::AXIS }); + if (dim == Dim::_3D) { + flds_bc_enum.push_back({ FldsBC::PERIODIC, FldsBC::PERIODIC }); + } + raise::ErrorIf(prtl_bc[0].size() != 2, + "invalid `grid.boundaries.particles`", + HERE); + prtl_bc_enum.push_back( + { PrtlBC::pick(fmt::toLower(prtl_bc[0][0]).c_str()), + PrtlBC::pick(fmt::toLower(prtl_bc[0][1]).c_str()) }); + prtl_bc_enum.push_back({ PrtlBC::AXIS, PrtlBC::AXIS }); + if (dim == Dim::_3D) { + prtl_bc_enum.push_back({ PrtlBC::PERIODIC, PrtlBC::PERIODIC }); + } + } else { + raise::ErrorIf(flds_bc[0].size() != 1, + "invalid `grid.boundaries.fields`", + HERE); + raise::ErrorIf(prtl_bc[0].size() != 1, + "invalid `grid.boundaries.particles`", + HERE); + flds_bc_enum.push_back( + { FldsBC::HORIZON, FldsBC::pick(fmt::toLower(flds_bc[0][0]).c_str()) }); + flds_bc_enum.push_back({ FldsBC::AXIS, FldsBC::AXIS }); + if (dim == Dim::_3D) { + flds_bc_enum.push_back({ FldsBC::PERIODIC, FldsBC::PERIODIC }); + } + prtl_bc_enum.push_back( + { PrtlBC::HORIZON, PrtlBC::pick(fmt::toLower(prtl_bc[0][0]).c_str()) }); + prtl_bc_enum.push_back({ PrtlBC::AXIS, PrtlBC::AXIS }); + if (dim == Dim::_3D) { + prtl_bc_enum.push_back({ PrtlBC::PERIODIC, PrtlBC::PERIODIC }); + } + } + } + + raise::ErrorIf(flds_bc_enum.size() != (std::size_t)dim, + "invalid inferred `grid.boundaries.fields`", + HERE); + raise::ErrorIf(prtl_bc_enum.size() != (std::size_t)dim, + "invalid inferred `grid.boundaries.particles`", + HERE); + boundaries_t flds_bc_pairwise; + boundaries_t prtl_bc_pairwise; + for (auto d { 0u }; d < (dim_t)dim; ++d) { + raise::ErrorIf( + flds_bc_enum[d].size() != 2, + fmt::format("invalid inferred `grid.boundaries.fields[%d]`", d), + HERE); + raise::ErrorIf( + prtl_bc_enum[d].size() != 2, + fmt::format("invalid inferred `grid.boundaries.particles[%d]`", d), + HERE); + flds_bc_pairwise.push_back({ flds_bc_enum[d][0], flds_bc_enum[d][1] }); + prtl_bc_pairwise.push_back({ prtl_bc_enum[d][0], prtl_bc_enum[d][1] }); + } + return { flds_bc_pairwise, prtl_bc_pairwise }; + } + + void Boundaries::read(Dimension dim, + const Coord& coord_enum, + const boundaries_t& extent_pairwise, + const toml::value& toml_data) { + if (needs_match_boundaries) { + if (coord_enum == Coord::Cart) { + auto min_extent = std::numeric_limits::max(); + for (const auto& e : extent_pairwise) { + min_extent = std::min(min_extent, e.second - e.first); + } + const auto default_ds = min_extent * defaults::bc::match::ds_frac; + try { + auto ds = toml::find(toml_data, "grid", "boundaries", "match", "ds"); + for (auto d = 0u; d < dim; ++d) { + match_ds_array.push_back({ ds, ds }); + } + } catch (...) { + try { + const auto ds = toml::find>>( + toml_data, + "grid", + "boundaries", + "match", + "ds"); + raise::ErrorIf(ds.size() != dim, + "invalid # in `grid.boundaries.match.ds`", + HERE); + for (auto d = 0u; d < dim; ++d) { + if (ds[d].size() == 1) { + match_ds_array.push_back({ ds[d][0], ds[d][0] }); + } else if (ds[d].size() == 2) { + match_ds_array.push_back({ ds[d][0], ds[d][1] }); + } else if (ds[d].size() == 0) { + match_ds_array.push_back({}); + } else { + raise::Error("invalid `grid.boundaries.match.ds`", HERE); + } + } + } catch (...) { + for (auto d = 0u; d < dim; ++d) { + match_ds_array.push_back({ default_ds, default_ds }); + } + } + } + } else { + auto r_extent = extent_pairwise[0].second - extent_pairwise[0].first; + const auto ds = toml::find_or( + toml_data, + "grid", + "boundaries", + "match", + "ds", + r_extent * defaults::bc::match::ds_frac); + match_ds_array.push_back({ ds, ds }); + } + } + + if (needs_absorb_boundaries) { + if (coord_enum == Coord::Cart) { + auto min_extent = std::numeric_limits::max(); + for (const auto& e : extent_pairwise) { + min_extent = std::min(min_extent, e.second - e.first); + } + absorb_ds = toml::find_or(toml_data, + "grid", + "boundaries", + "absorb", + "ds", + min_extent * defaults::bc::absorb::ds_frac); + } else { + auto r_extent = extent_pairwise[0].second - extent_pairwise[0].first; + absorb_ds = toml::find_or(toml_data, + "grid", + "boundaries", + "absorb", + "ds", + r_extent * defaults::bc::absorb::ds_frac); + } + } + + if (needs_atmosphere_boundaries) { + atmosphere_temperature = toml::find(toml_data, + "grid", + "boundaries", + "atmosphere", + "temperature"); + atmosphere_height = toml::find(toml_data, + "grid", + "boundaries", + "atmosphere", + "height"); + atmosphere_density = toml::find(toml_data, + "grid", + "boundaries", + "atmosphere", + "density"); + atmosphere_ds = + toml::find_or(toml_data, "grid", "boundaries", "atmosphere", "ds", ZERO); + atmosphere_g = atmosphere_temperature / atmosphere_height; + atmosphere_species = toml::find>( + toml_data, + "grid", + "boundaries", + "atmosphere", + "species"); + } + } + + void Boundaries::setParams(SimulationParams* params) const { + if (needs_match_boundaries) { + params->set("grid.boundaries.match.ds", match_ds_array); + } + if (needs_absorb_boundaries) { + params->set("grid.boundaries.absorb.ds", absorb_ds); + } + if (needs_atmosphere_boundaries) { + params->set("grid.boundaries.atmosphere.temperature", + atmosphere_temperature); + params->set("grid.boundaries.atmosphere.density", atmosphere_density); + params->set("grid.boundaries.atmosphere.height", atmosphere_height); + params->set("grid.boundaries.atmosphere.ds", atmosphere_ds); + params->set("grid.boundaries.atmosphere.g", atmosphere_g); + params->set("grid.boundaries.atmosphere.species", atmosphere_species); + } + } + + void Grid::read(const SimEngine& engine_enum, const toml::value& toml_data) { + /* domain decomposition ------------------------------------------------ */ + int default_ndomains = 1; +#if defined(MPI_ENABLED) + raise::ErrorIf(MPI_Comm_size(MPI_COMM_WORLD, &default_ndomains) != MPI_SUCCESS, + "MPI_Comm_size failed", + HERE); +#endif + number_of_domains = toml::find_or(toml_data, + "simulation", + "domain", + "number", + default_ndomains); + + domain_decomposition = toml::find_or>( + toml_data, + "simulation", + "domain", + "decomposition", + std::vector { -1, -1, -1 }); + + /* resolution and dimension ------------------------------------------- */ + resolution = toml::find>(toml_data, "grid", "resolution"); + raise::ErrorIf(resolution.size() < 1 || resolution.size() > 3, + "invalid `grid.resolution`", + HERE); + dim = static_cast(resolution.size()); + + if (domain_decomposition.size() > dim) { + domain_decomposition.erase(domain_decomposition.begin() + (std::size_t)(dim), + domain_decomposition.end()); + } + raise::ErrorIf(domain_decomposition.size() != dim, + "invalid `simulation.domain.decomposition`", + HERE); + + /* metric and coordinates -------------------------------------------- */ + metric_enum = Metric::pick( + fmt::toLower(toml::find(toml_data, "grid", "metric", "metric")) + .c_str()); + std::string coord; + if (metric_enum == Metric::Minkowski) { + raise::ErrorIf(engine_enum != SimEngine::SRPIC, + "minkowski metric is only supported for SRPIC", + HERE); + coord = "cart"; + } else if (metric_enum == Metric::QKerr_Schild or + metric_enum == Metric::QSpherical) { + // quasi-spherical geometry + raise::ErrorIf(dim == Dim::_1D, + "not enough dimensions for qspherical geometry", + HERE); + raise::ErrorIf(dim == Dim::_3D, + "3D not implemented for qspherical geometry", + HERE); + coord = "qsph"; + metric_params["qsph_r0"] = toml::find_or(toml_data, + "grid", + "metric", + "qsph_r0", + defaults::qsph::r0); + metric_params["qsph_h"] = toml::find_or(toml_data, + "grid", + "metric", + "qsph_h", + defaults::qsph::h); + } else { + // spherical geometry + raise::ErrorIf(dim == Dim::_1D, + "not enough dimensions for spherical geometry", + HERE); + raise::ErrorIf(dim == Dim::_3D, + "3D not implemented for spherical geometry", + HERE); + coord = "sph"; + } + if ((engine_enum == SimEngine::GRPIC) && + (metric_enum != Metric::Kerr_Schild_0)) { + const auto ks_a = toml::find_or(toml_data, + "grid", + "metric", + "ks_a", + defaults::ks::a); + metric_params["ks_a"] = ks_a; + metric_params["ks_rh"] = ONE + math::sqrt(ONE - SQR(ks_a)); + } + coord_enum = Coord::pick(coord.c_str()); + + /* extent ------------------------------------------------------------- */ + extent = toml::find>>(toml_data, + "grid", + "extent"); + + if (extent.size() > dim) { + extent.erase(extent.begin() + (std::size_t)(dim), extent.end()); + } + raise::ErrorIf(extent[0].size() != 2, "invalid `grid.extent[0]`", HERE); + if (coord_enum != Coord::Cart) { + raise::ErrorIf(extent.size() > 1, + "invalid `grid.extent` for non-cartesian geometry", + HERE); + extent.push_back({ ZERO, constant::PI }); + if (dim == Dim::_3D) { + extent.push_back({ ZERO, TWO * constant::PI }); + } + } + raise::ErrorIf(extent.size() != dim, "invalid inferred `grid.extent`", HERE); + for (auto d { 0u }; d < (dim_t)dim; ++d) { + raise::ErrorIf(extent[d].size() != 2, + fmt::format("invalid inferred `grid.extent[%d]`", d), + HERE); + extent_pairwise_.push_back({ extent[d][0], extent[d][1] }); + } + + /* metric parameters ------------------------------------------------------ */ + if (coord_enum == Coord::Qsph) { + metric_params_short_["r0"] = metric_params["qsph_r0"]; + metric_params_short_["h"] = metric_params["qsph_h"]; + } + if ((engine_enum == SimEngine::GRPIC) && + (metric_enum != Metric::Kerr_Schild_0)) { + metric_params_short_["a"] = metric_params["ks_a"]; + } + // set("grid.metric.params", params); + + std::pair dx0_V0; + if (metric_enum == Metric::Minkowski) { + if (dim == Dim::_1D) { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } else if (dim == Dim::_2D) { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } else { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } + } else if (metric_enum == Metric::Spherical) { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } else if (metric_enum == Metric::QSpherical) { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } else if (metric_enum == Metric::Kerr_Schild) { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } else if (metric_enum == Metric::Kerr_Schild_0) { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } else if (metric_enum == Metric::QKerr_Schild) { + dx0_V0 = get_dx0_V0>(resolution, + extent_pairwise_, + metric_params_short_); + } + auto [dx0, V0] = dx0_V0; + scale_dx0 = dx0; + scale_V0 = V0; + } + + void Grid::setParams(SimulationParams* params) const { + params->set("simulation.domain.number", number_of_domains); + params->set("simulation.domain.decomposition", domain_decomposition); + + params->set("grid.resolution", resolution); + params->set("grid.dim", dim); + params->set("grid.metric.metric", metric_enum); + params->set("grid.metric.coord", coord_enum); + for (const auto& [key, value] : metric_params) { + params->set("grid.metric." + key, value); + } + params->set("grid.metric.params", metric_params_short_); + params->set("grid.extent", extent_pairwise_); + + params->set("scales.dx0", scale_dx0); + params->set("scales.V0", scale_V0); + } + + } // namespace params +} // namespace ntt diff --git a/src/framework/parameters/grid.h b/src/framework/parameters/grid.h new file mode 100644 index 00000000..b436fd9d --- /dev/null +++ b/src/framework/parameters/grid.h @@ -0,0 +1,88 @@ +/** + * @file framework/parameters/grid.h + * @brief Auxiliary functions for reading in grid/box parameters + * @implements + * - ntt::params::Boundaries + * - ntt::params::GetBoundaryConditions -> (boundaries_t, boundaries_t) + * @cpp: + * - grid.cpp + * @namespaces: + * - ntt::params:: + */ +#ifndef FRAMEWORK_PARAMETERS_GRID_H +#define FRAMEWORK_PARAMETERS_GRID_H + +#include "enums.h" +#include "global.h" + +#include "utils/toml.h" + +#include "framework/parameters/parameters.h" + +#include +#include +#include + +namespace ntt { + namespace params { + + struct Boundaries { + const bool needs_match_boundaries; + boundaries_t match_ds_array; + + const bool needs_absorb_boundaries; + real_t absorb_ds; + + const bool needs_atmosphere_boundaries; + real_t atmosphere_temperature; + real_t atmosphere_height; + real_t atmosphere_density; + real_t atmosphere_g; + real_t atmosphere_ds; + std::pair atmosphere_species; + + Boundaries(bool needs_match, bool needs_absorb, bool needs_atmosphere) + : needs_match_boundaries { needs_match } + , needs_absorb_boundaries { needs_absorb } + , needs_atmosphere_boundaries { needs_atmosphere } {} + + void read(Dimension, + const Coord&, + const boundaries_t&, + const toml::value&); + void setParams(SimulationParams*) const; + }; + + struct Grid { + int number_of_domains; + std::vector domain_decomposition; + + std::vector resolution; + Dimension dim; + + std::vector> extent; + boundaries_t extent_pairwise_; + + Metric metric_enum = Metric::INVALID; + Coord coord_enum = Coord::INVALID; + std::map metric_params; + std::map metric_params_short_; + + real_t scale_dx0; + real_t scale_V0; + + void read(const SimEngine&, const toml::value&); + void setParams(SimulationParams*) const; + }; + + auto GetBoundaryConditions( + SimulationParams* params, + const SimEngine&, + Dimension, + const Coord&, + const toml::value&) -> std::tuple, boundaries_t>; + + } // namespace params +} // namespace ntt + +#endif // FRAMEWORK_PARAMETERS_GRID_H diff --git a/src/framework/parameters/output.cpp b/src/framework/parameters/output.cpp new file mode 100644 index 00000000..17fdf540 --- /dev/null +++ b/src/framework/parameters/output.cpp @@ -0,0 +1,203 @@ +#include "framework/parameters/output.h" + +#include "defaults.h" +#include "global.h" + +#include "utils/error.h" +#include "utils/log.h" +#include "utils/toml.h" + +namespace ntt { + namespace params { + + void Output::read(Dimension dim, std::size_t nspec, const toml::value& toml_data) { + format = toml::find_or(toml_data, "output", "format", defaults::output::format); + global_interval = toml::find_or(toml_data, + "output", + "interval", + defaults::output::interval); + global_interval_time = toml::find_or(toml_data, + "output", + "interval_time", + -1.0); + raise::ErrorIf( + not toml::find_or(toml_data, "output", "separate_files", true), + "separate_files=false is deprecated", + HERE); + + for (const auto& category : { "fields", "particles", "spectra", "stats" }) { + const auto q_int = toml::find_or(toml_data, + "output", + category, + "interval", + 0); + const auto q_int_time = toml::find_or(toml_data, + "output", + category, + "interval_time", + -1.0); + categories[category].enable = toml::find_or(toml_data, + "output", + category, + "enable", + true); + if ((q_int == 0) and (q_int_time == -1.0)) { + categories[category].interval = global_interval; + categories[category].interval_time = global_interval_time; + } else { + categories[category].interval = q_int; + categories[category].interval_time = q_int_time; + } + } + + /* Fields --------------------------------------------------------------- */ + const auto flds_out = toml::find_or(toml_data, + "output", + "fields", + "quantities", + std::vector {}); + const auto custom_flds_out = toml::find_or(toml_data, + "output", + "fields", + "custom", + std::vector {}); + if (flds_out.size() == 0) { + raise::Warning("No fields output specified", HERE); + } + fields_quantities = flds_out; + fields_custom_quantities = custom_flds_out; + fields_mom_smooth = toml::find_or(toml_data, + "output", + "fields", + "mom_smooth", + defaults::output::mom_smooth); + try { + auto field_dwn_ = toml::find>(toml_data, + "output", + "fields", + "downsampling"); + for (auto i = 0u; i < field_dwn_.size(); ++i) { + fields_downsampling.push_back(field_dwn_[i]); + } + } catch (...) { + try { + auto field_dwn_ = toml::find(toml_data, + "output", + "fields", + "downsampling"); + for (auto i = 0u; i < dim; ++i) { + fields_downsampling.push_back(field_dwn_); + } + } catch (...) { + for (auto i = 0u; i < dim; ++i) { + fields_downsampling.push_back(1u); + } + } + } + raise::ErrorIf(fields_downsampling.size() > 3, + "invalid `output.fields.downsampling`", + HERE); + if (fields_downsampling.size() > dim) { + fields_downsampling.erase(fields_downsampling.begin() + (std::size_t)(dim), + fields_downsampling.end()); + } + for (const auto& dwn : fields_downsampling) { + raise::ErrorIf(dwn == 0, "downsampling factor must be nonzero", HERE); + } + + /* Particles ------------------------------------------------------------ */ + auto all_specs = std::vector {}; + for (auto i = 0u; i < nspec; ++i) { + all_specs.push_back(static_cast(i + 1)); + } + particles_species = toml::find_or(toml_data, + "output", + "particles", + "species", + all_specs); + particles_stride = toml::find_or(toml_data, + "output", + "particles", + "stride", + defaults::output::prtl_stride); + + /* Spectra -------------------------------------------------------------- */ + spectra_e_min = toml::find_or(toml_data, + "output", + "spectra", + "e_min", + defaults::output::spec_emin); + spectra_e_max = toml::find_or(toml_data, + "output", + "spectra", + "e_max", + defaults::output::spec_emax); + spectra_log_bins = toml::find_or(toml_data, + "output", + "spectra", + "log_bins", + defaults::output::spec_log); + spectra_n_bins = toml::find_or(toml_data, + "output", + "spectra", + "n_bins", + defaults::output::spec_nbins); + + /* Stats ---------------------------------------------------------------- */ + stats_quantities = toml::find_or(toml_data, + "output", + "stats", + "quantities", + defaults::output::stats_quantities); + stats_custom_quantities = toml::find_or(toml_data, + "output", + "stats", + "custom", + std::vector {}); + + /* Debug ---------------------------------------------------------------- */ + debug_as_is = toml::find_or(toml_data, "output", "debug", "as_is", false); + debug_ghosts = toml::find_or(toml_data, "output", "debug", "ghosts", false); + if (debug_ghosts) { + for (const auto& dwn : fields_downsampling) { + raise::ErrorIf( + dwn != 1, + "full resolution required when outputting with ghost cells", + HERE); + } + } + } + + void Output::setParams(SimulationParams* params) const { + params->set("output.format", format); + params->set("output.interval", global_interval); + params->set("output.interval_time", global_interval_time); + for (const auto& [category, cat_params] : categories) { + params->set("output." + category + ".enable", cat_params.enable); + params->set("output." + category + ".interval", cat_params.interval); + params->set("output." + category + ".interval_time", + cat_params.interval_time); + } + + params->set("output.fields.quantities", fields_quantities); + params->set("output.fields.custom", fields_custom_quantities); + params->set("output.fields.mom_smooth", fields_mom_smooth); + params->set("output.fields.downsampling", fields_downsampling); + + params->set("output.particles.species", particles_species); + params->set("output.particles.stride", particles_stride); + + params->set("output.spectra.e_min", spectra_e_min); + params->set("output.spectra.e_max", spectra_e_max); + params->set("output.spectra.log_bins", spectra_log_bins); + params->set("output.spectra.n_bins", spectra_n_bins); + + params->set("output.stats.quantities", stats_quantities); + params->set("output.stats.custom", stats_custom_quantities); + + params->set("output.debug.as_is", debug_as_is); + params->set("output.debug.ghosts", debug_ghosts); + } + + } // namespace params +} // namespace ntt diff --git a/src/framework/parameters/output.h b/src/framework/parameters/output.h new file mode 100644 index 00000000..7a15baab --- /dev/null +++ b/src/framework/parameters/output.h @@ -0,0 +1,68 @@ +/** + * @file framework/parameters/output.h + * @brief Auxiliary functions for reading in output parameters + * @implements + * - ntt::params::Output + * - ntt::params::OutputCategory + * @cpp: + * - output.cpp + * @namespaces: + * - ntt::params:: + */ +#ifndef FRAMEWORK_PARAMETERS_OUTPUT_H +#define FRAMEWORK_PARAMETERS_OUTPUT_H + +#include "global.h" + +#include "utils/toml.h" + +#include "framework/parameters/parameters.h" + +#include +#include +#include + +namespace ntt { + namespace params { + + struct OutputCategory { + bool enable; + timestep_t interval; + simtime_t interval_time; + }; + + struct Output { + std::string format; + + timestep_t global_interval; + simtime_t global_interval_time; + + std::map categories; + + std::vector fields_quantities; + std::vector fields_custom_quantities; + unsigned short fields_mom_smooth; + std::vector fields_downsampling; + + std::vector particles_species; + npart_t particles_stride; + + real_t spectra_e_min; + real_t spectra_e_max; + bool spectra_log_bins; + std::size_t spectra_n_bins; + + std::vector stats_quantities; + std::vector stats_custom_quantities; + + bool debug_as_is; + bool debug_ghosts; + + void read(Dimension, std::size_t, const toml::value&); + void setParams(SimulationParams*) const; + }; + + } // namespace params +} // namespace ntt + +#endif // FRAMEWORK_PARAMETERS_OUTPUT_H diff --git a/src/framework/parameters/parameters.cpp b/src/framework/parameters/parameters.cpp new file mode 100644 index 00000000..9d4440d0 --- /dev/null +++ b/src/framework/parameters/parameters.cpp @@ -0,0 +1,302 @@ +#include "framework/parameters/parameters.h" + +#include "defaults.h" +#include "enums.h" +#include "global.h" + +#include "utils/error.h" +#include "utils/formatting.h" +#include "utils/numeric.h" +#include "utils/toml.h" + +#include "framework/containers/species.h" +#include "framework/parameters/algorithms.h" +#include "framework/parameters/extra.h" +#include "framework/parameters/grid.h" +#include "framework/parameters/output.h" +#include "framework/parameters/particles.h" + +#if defined(MPI_ENABLED) + #include +#endif + +#include +#include +#include +#include +#include + +namespace ntt { + + /* + * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + * Parameters that must not be changed during the checkpoint restart + * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + */ + void SimulationParams::setImmutableParams(const toml::value& toml_data) { + /* [simulation] --------------------------------------------------------- */ + const auto engine_enum = SimEngine::pick( + fmt::toLower(toml::find(toml_data, "simulation", "engine")).c_str()); + set("simulation.engine", engine_enum); + + /* grid and decomposition ------------------------------------------------ */ + params::Grid grid_params {}; + grid_params.read(engine_enum, toml_data); + grid_params.setParams(this); + + /* [scales] ------------------------------------------------------------- */ + const auto larmor0 = toml::find(toml_data, "scales", "larmor0"); + const auto skindepth0 = toml::find(toml_data, "scales", "skindepth0"); + raise::ErrorIf(larmor0 <= ZERO || skindepth0 <= ZERO, + "larmor0 and skindepth0 must be positive", + HERE); + set("scales.larmor0", larmor0); + set("scales.skindepth0", skindepth0); + set("scales.sigma0", SQR(skindepth0 / larmor0)); + set("scales.B0", ONE / larmor0); + set("scales.omegaB0", ONE / larmor0); + + /* [particles] ---------------------------------------------------------- */ + const auto ppc0 = toml::find(toml_data, "particles", "ppc0"); + set("particles.ppc0", ppc0); + raise::ErrorIf(ppc0 <= 0.0, "ppc0 must be positive", HERE); + set("particles.use_weights", + toml::find_or(toml_data, "particles", "use_weights", false)); + + set("scales.n0", ppc0 / get("scales.V0")); + set("scales.q0", get("scales.V0") / (ppc0 * SQR(skindepth0))); + + /* [particles.species] -------------------------------------------------- */ + std::vector species; + const auto species_tab = toml::find_or(toml_data, + "particles", + "species", + toml::array {}); + set("particles.nspec", species_tab.size()); + + spidx_t idx = 1; + for (const auto& sp : species_tab) { + species.emplace_back(params::GetParticleSpecies(this, engine_enum, idx, sp)); + idx += 1; + } + set("particles.species", species); + } + + /* + * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + * Parameters that may be changed during the checkpoint restart + * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + */ + void SimulationParams::setMutableParams(const toml::value& toml_data) { + const auto engine_enum = get("simulation.engine"); + const auto coord_enum = get("grid.metric.coord"); + const auto dim = get("grid.dim"); + const auto extent_pairwise = get>("grid.extent"); + + /* [simulation] --------------------------------------------------------- */ + set("simulation.name", + toml::find(toml_data, "simulation", "name")); + set("simulation.runtime", + toml::find(toml_data, "simulation", "runtime")); + + /* [grid.boundaraies] --------------------------------------------------- */ + const auto [flds_bc, prtl_bc] = params::GetBoundaryConditions(this, + engine_enum, + dim, + coord_enum, + toml_data); + set("grid.boundaries.fields", flds_bc); + set("grid.boundaries.particles", prtl_bc); + + /* [particles] ---------------------------------------------------------- */ + set("particles.clear_interval", + toml::find_or(toml_data, "particles", "clear_interval", defaults::clear_interval)); + const auto species_tab = toml::find_or(toml_data, + "particles", + "species", + toml::array {}); + std::vector species = get>( + "particles.species"); + raise::ErrorIf(species_tab.size() != species.size(), + "number of species changed after restart", + HERE); + + std::vector new_species; + + spidx_t idxM1 = 0; + for (const auto& sp : species_tab) { + const auto maxnpart_real = toml::find(sp, "maxnpart"); + const auto maxnpart = static_cast(maxnpart_real); + const auto particle_species = species[idxM1]; + new_species.emplace_back(particle_species.index(), + particle_species.label(), + particle_species.mass(), + particle_species.charge(), + maxnpart, + particle_species.pusher(), + particle_species.use_tracking(), + particle_species.radiative_drag_flags(), + particle_species.emission_policy_flag(), + particle_species.npld_r(), + particle_species.npld_i()); + idxM1++; + } + set("particles.species", new_species); + + /* [output] ------------------------------------------------------------- */ + params::Output output_params; + output_params.read(dim, get("particles.nspec"), toml_data); + output_params.setParams(this); + + /* [checkpoint] --------------------------------------------------------- */ + set("checkpoint.interval", + toml::find_or(toml_data, + "checkpoint", + "interval", + defaults::checkpoint::interval)); + set("checkpoint.interval_time", + toml::find_or(toml_data, "checkpoint", "interval_time", -1.0)); + set("checkpoint.keep", + toml::find_or(toml_data, "checkpoint", "keep", defaults::checkpoint::keep)); + auto walltime_str = toml::find_or(toml_data, + "checkpoint", + "walltime", + defaults::checkpoint::walltime); + if (walltime_str.empty()) { + walltime_str = defaults::checkpoint::walltime; + } + set("checkpoint.walltime", walltime_str); + + const auto checkpoint_write_path = toml::find_or( + toml_data, + "checkpoint", + "write_path", + fmt::format(defaults::checkpoint::write_path.c_str(), + get("simulation.name").c_str())); + set("checkpoint.write_path", checkpoint_write_path); + set("checkpoint.read_path", + toml::find_or(toml_data, "checkpoint", "read_path", checkpoint_write_path)); + + /* [diagnostics] -------------------------------------------------------- */ + set("diagnostics.interval", + toml::find_or(toml_data, "diagnostics", "interval", defaults::diag::interval)); + set("diagnostics.blocking_timers", + toml::find_or(toml_data, "diagnostics", "blocking_timers", false)); + set("diagnostics.colored_stdout", + toml::find_or(toml_data, "diagnostics", "colored_stdout", false)); + set("diagnostics.log_level", + toml::find_or(toml_data, "diagnostics", "log_level", defaults::diag::log_level)); + + /* inferred variables --------------------------------------------------- */ + params::Boundaries boundaries_params { + isPromised("grid.boundaries.match.ds"), + isPromised("grid.boundaries.absorb.ds"), + isPromised("grid.boundaries.atmosphere.temperature") + }; + boundaries_params.read(dim, coord_enum, extent_pairwise, toml_data); + boundaries_params.setParams(this); + + /* [algorithms] --------------------------------------------------------- */ + ntt::params::Algorithms alg_params {}; + std::map alg_extra_flags = { + { "gr", engine_enum == SimEngine::GRPIC }, + { "gca", isPromised("algorithms.gca.e_ovr_b_max") }, + }; + alg_params.read(get("scales.dx0"), alg_extra_flags, toml_data); + alg_params.setParams(alg_extra_flags, this); + + /* extra physics ------------------------------------------------------ */ + ntt::params::Extra extra_params {}; + std::map extra_extra_flags = { + { "synchrotron_drag",isPromised("radiation.drag.synchrotron.gamma_rad") }, + { "compton_drag", isPromised("radiation.drag.compton.gamma_rad") }, + { "synchrotron_emission", + isPromised("radiation.emission.synchrotron.photon_species") }, + { "compton_emission", isPromised("radiation.emission.compton.photon_species") } + }; + extra_params.read(extra_extra_flags, toml_data); + extra_params.setParams(extra_extra_flags, this); + + // @TODO: disabling stats for non-Cartesian + if (coord_enum != Coord::Cart) { + set("output.stats.enable", false); + } + } + + void SimulationParams::setSetupParams(const toml::value& toml_data) { + /* [setup] -------------------------------------------------------------- */ + const auto setup = toml::find_or(toml_data, "setup", toml::table {}); + for (const auto& [key, val] : setup) { + if (val.is_boolean()) { + set("setup." + key, (bool)(val.as_boolean())); + } else if (val.is_integer()) { + set("setup." + key, (int)(val.as_integer())); + } else if (val.is_floating()) { + set("setup." + key, (real_t)(val.as_floating())); + } else if (val.is_string()) { + set("setup." + key, (std::string)(val.as_string())); + } else if (val.is_array()) { + const auto val_arr = val.as_array(); + if (val_arr.size() == 0) { + continue; + } else { + if (val_arr[0].is_integer()) { + std::vector vec; + for (const auto& v : val_arr) { + vec.push_back(v.as_integer()); + } + set("setup." + key, vec); + } else if (val_arr[0].is_floating()) { + std::vector vec; + for (const auto& v : val_arr) { + vec.push_back(v.as_floating()); + } + set("setup." + key, vec); + } else if (val_arr[0].is_boolean()) { + std::vector vec; + for (const auto& v : val_arr) { + vec.push_back(v.as_boolean()); + } + set("setup." + key, vec); + } else if (val_arr[0].is_string()) { + std::vector vec; + for (const auto& v : val_arr) { + vec.push_back(v.as_string()); + } + set("setup." + key, vec); + } else if (val_arr[0].is_array()) { + raise::Error("only 1D arrays allowed in [setup]", HERE); + } else { + raise::Error("invalid setup variable type", HERE); + } + } + } + } + } + + void SimulationParams::setCheckpointParams(bool is_resuming, + timestep_t start_step, + simtime_t start_time) { + set("checkpoint.is_resuming", is_resuming); + set("checkpoint.start_step", start_step); + set("checkpoint.start_time", start_time); + } + + void SimulationParams::checkPromises() const { + raise::ErrorIf(!promisesFulfilled(), + "Have not defined all the necessary variables", + HERE); + } + + void SimulationParams::saveTOML(const std::string& path, simtime_t time) const { + CallOnce([&]() { + std::ofstream metadata; + metadata.open(path); + metadata << "[metadata]\n" + << " time = " << time << "\n\n" + << data() << std::endl; + metadata.close(); + }); + } + +} // namespace ntt diff --git a/src/framework/parameters.h b/src/framework/parameters/parameters.h similarity index 89% rename from src/framework/parameters.h rename to src/framework/parameters/parameters.h index 0af4fa40..e2a1d75e 100644 --- a/src/framework/parameters.h +++ b/src/framework/parameters/parameters.h @@ -1,5 +1,5 @@ /** - * @file framework/parameters.h + * @file framework/parameters/parameters.h * @brief Structure for defining and holding initial simulation parameters * @implements * - ntt::SimulationParams : ntt::Parameters @@ -14,8 +14,8 @@ * @note A proper metric is used to infer the minimum cell size/volume etc. */ -#ifndef FRAMEWORK_PARAMETERS_H -#define FRAMEWORK_PARAMETERS_H +#ifndef FRAMEWORK_PARAMETERS_PARAMETERS_H +#define FRAMEWORK_PARAMETERS_PARAMETERS_H #include "utils/param_container.h" #include "utils/toml.h" @@ -62,4 +62,4 @@ namespace ntt { } // namespace ntt -#endif // FRAMEWORK_PARAMETERS_H +#endif // FRAMEWORK_PARAMETERS_PARAMETERS_H diff --git a/src/framework/parameters/particles.cpp b/src/framework/parameters/particles.cpp new file mode 100644 index 00000000..ddc22135 --- /dev/null +++ b/src/framework/parameters/particles.cpp @@ -0,0 +1,219 @@ +#include "framework/parameters/particles.h" + +#include "defaults.h" +#include "enums.h" +#include "global.h" + +#include "utils/error.h" +#include "utils/formatting.h" +#include "utils/toml.h" + +#include "framework/containers/species.h" +#include "framework/parameters/parameters.h" + +#include + +namespace ntt { + + namespace params { + + /* + * Auxiliary functions + */ + auto getRadiativeDragFlags(const std::string& radiative_drag_str) + -> RadiativeDragFlags { + if (fmt::toLower(radiative_drag_str) == "none") { + return RadiativeDrag::NONE; + } else { + // separate comas + RadiativeDragFlags flags = RadiativeDrag::NONE; + std::string token; + std::istringstream tokenStream(radiative_drag_str); + while (std::getline(tokenStream, token, ',')) { + const auto token_lower = fmt::toLower(token); + if (token_lower == "synchrotron") { + flags |= RadiativeDrag::SYNCHROTRON; + } else if (token_lower == "compton") { + flags |= RadiativeDrag::COMPTON; + } else { + raise::Error(fmt::format("Invalid radiative_drag value: %s", + radiative_drag_str.c_str()), + HERE); + } + } + return flags; + } + } + + auto getPusherFlags(const std::string& particle_pusher_str) + -> ParticlePusherFlags { + if (fmt::toLower(particle_pusher_str) == "none") { + return ParticlePusher::NONE; + } else { + // separate comas + ParticlePusherFlags flags = ParticlePusher::NONE; + std::string token; + std::istringstream tokenStream(particle_pusher_str); + while (std::getline(tokenStream, token, ',')) { + const auto token_lower = fmt::toLower(token); + if (token_lower == "photon") { + flags |= ParticlePusher::PHOTON; + } else if (token_lower == "boris") { + flags |= ParticlePusher::BORIS; + } else if (token_lower == "vay") { + flags |= ParticlePusher::VAY; + } else if (token_lower == "gca") { + flags |= ParticlePusher::GCA; + } else { + raise::Error(fmt::format("Invalid pusher value: %s", + particle_pusher_str.c_str()), + HERE); + } + } + if (flags & ParticlePusher::PHOTON and flags & ParticlePusher::GCA) { + raise::Error("Photon pusher cannot be used with GCA", HERE); + } + return flags; + } + } + + auto getEmissionPolicyFlag(const std::string& emission_policy_str) + -> EmissionTypeFlag { + if (fmt::toLower(emission_policy_str) == "none") { + return EmissionType::NONE; + } else if (fmt::toLower(emission_policy_str) == "synchrotron") { + return EmissionType::SYNCHROTRON; + } else if (fmt::toLower(emission_policy_str) == "compton") { + return EmissionType::COMPTON; + } else if (fmt::toLower(emission_policy_str) == "strongfieldpp") { + return EmissionType::STRONGFIELDPP; + } else { + raise::Error(fmt::format("Invalid emission_policy value: %s", + emission_policy_str.c_str()), + HERE); + return EmissionType::NONE; + } + } + + auto GetParticleSpecies(SimulationParams* params, + const SimEngine& engine_enum, + spidx_t idx, + const toml::value& sp) -> ParticleSpecies { + const auto label = toml::find_or(sp, + "label", + "s" + std::to_string(idx)); + const auto mass = toml::find(sp, "mass"); + const auto charge = toml::find(sp, "charge"); + raise::ErrorIf((charge != 0.0f) && (mass == 0.0f), + "mass of the charged species must be non-zero", + HERE); + const auto is_massless = (mass == 0.0f) && (charge == 0.0f); + const auto def_pusher = (is_massless ? defaults::ph_pusher + : defaults::em_pusher); + const auto maxnpart_real = toml::find(sp, "maxnpart"); + const auto maxnpart = static_cast(maxnpart_real); + auto pusher_str = toml::find_or(sp, "pusher", std::string(def_pusher)); + const auto npayloads_real = toml::find_or(sp, + "n_payloads_real", + static_cast(0)); + const auto use_tracking = toml::find_or(sp, "tracking", false); + auto npayloads_int = toml::find_or(sp, + "n_payloads_int", + static_cast(0)); + if (use_tracking) { +#if !defined(MPI_ENABLED) + npayloads_int += 1; +#else + npayloads_int += 2; +#endif + } + auto radiative_drag_str = toml::find_or(sp, + "radiative_drag", + std::string("default")); + + const auto radiative_drag_defaulted = (fmt::toLower(radiative_drag_str) == + "default"); + if (radiative_drag_defaulted) { + radiative_drag_str = "none"; + } + + const auto emission_policy_str = toml::find_or(sp, + "emission", + std::string("none")); + raise::ErrorIf((fmt::toLower(radiative_drag_str) != "none") && is_massless, + "radiative drag is only applicable to massive particles", + HERE); + raise::ErrorIf((fmt::toLower(pusher_str) == "photon") && !is_massless, + "photon pusher is only applicable to massless particles", + HERE); + + auto particle_pusher_flags = getPusherFlags(pusher_str); + auto radiative_drag_flags = getRadiativeDragFlags(radiative_drag_str); + auto emission_policy_flag = getEmissionPolicyFlag(emission_policy_str); + + raise::ErrorIf((emission_policy_flag == EmissionType::STRONGFIELDPP) and + (not is_massless), + "Strong Field Pair Production emission policy is only " + "applicable to massless particles", + HERE); + raise::ErrorIf((emission_policy_flag == EmissionType::SYNCHROTRON or + emission_policy_flag == EmissionType::COMPTON) and + is_massless, + "Radiative emission policies are only applicable to " + "massive particles", + HERE); + + if (radiative_drag_defaulted) { + if (emission_policy_flag == EmissionType::SYNCHROTRON) { + radiative_drag_flags |= RadiativeDrag::SYNCHROTRON; + } else if (emission_policy_flag == EmissionType::COMPTON) { + radiative_drag_flags |= RadiativeDrag::COMPTON; + } + } + + if (radiative_drag_flags & RadiativeDrag::SYNCHROTRON) { + raise::ErrorIf(engine_enum != SimEngine::SRPIC, + "Synchrotron radiative drag is only supported for SRPIC", + HERE); + params->promiseToDefine("radiation.drag.synchrotron.gamma_rad"); + } + if (radiative_drag_flags & RadiativeDrag::COMPTON) { + raise::ErrorIf( + engine_enum != SimEngine::SRPIC, + "Inverse Compton radiative drag is only supported for SRPIC", + HERE); + params->promiseToDefine("radiation.drag.compton.gamma_rad"); + } + if (particle_pusher_flags & ParticlePusher::GCA) { + raise::ErrorIf(engine_enum != SimEngine::SRPIC, + "GCA pushers are only supported for SRPIC", + HERE); + params->promiseToDefine("algorithms.gca.e_ovr_b_max"); + params->promiseToDefine("algorithms.gca.larmor_max"); + } + + if (emission_policy_flag == EmissionType::SYNCHROTRON) { + params->promiseToDefine( + "radiation.emission.synchrotron.photon_species"); + params->promiseToDefine("radiation.drag.synchrotron.gamma_rad"); + } else if (emission_policy_flag == EmissionType::COMPTON) { + params->promiseToDefine("radiation.emission.compton.photon_species"); + params->promiseToDefine("radiation.drag.compton.gamma_rad"); + } + + return ParticleSpecies(idx, + label, + mass, + charge, + maxnpart, + particle_pusher_flags, + use_tracking, + radiative_drag_flags, + emission_policy_flag, + npayloads_real, + npayloads_int); + } + + } // namespace params + +} // namespace ntt diff --git a/src/framework/parameters/particles.h b/src/framework/parameters/particles.h new file mode 100644 index 00000000..6ceaebcd --- /dev/null +++ b/src/framework/parameters/particles.h @@ -0,0 +1,34 @@ +/** + * @file framework/parameters/particles.h + * @brief Auxiliary functions for reading in particle species parameters + * @implements + * - ntt::params::GetParticleSpecies -> ParticleSpecies + * @cpp: + * - particles.cpp + * @namespaces: + * - ntt::params:: + */ +#ifndef FRAMEWORK_PARAMETERS_PARTICLES_H +#define FRAMEWORK_PARAMETERS_PARTICLES_H + +#include "enums.h" + +#include "utils/toml.h" + +#include "framework/containers/species.h" +#include "framework/parameters/parameters.h" + +namespace ntt { + + namespace params { + + auto GetParticleSpecies(SimulationParams*, + const SimEngine&, + spidx_t, + const toml::value&) -> ParticleSpecies; + + } // namespace params + +} // namespace ntt + +#endif diff --git a/src/framework/simulation.h b/src/framework/simulation.h index d8696532..9f066440 100644 --- a/src/framework/simulation.h +++ b/src/framework/simulation.h @@ -16,11 +16,11 @@ #include "enums.h" -#include "arch/traits.h" #include "utils/error.h" #include "utils/toml.h" -#include "framework/parameters.h" +#include "engines/traits.h" +#include "framework/parameters/parameters.h" namespace ntt { diff --git a/src/framework/tests/parameters.cpp b/src/framework/tests/parameters.cpp index 78a3294b..41f854d3 100644 --- a/src/framework/tests/parameters.cpp +++ b/src/framework/tests/parameters.cpp @@ -1,4 +1,4 @@ -#include "framework/parameters.h" +#include "framework/parameters/parameters.h" #include "defaults.h" #include "enums.h" @@ -125,6 +125,21 @@ const auto sph_2d = u8R"( larmor0 = 0.01 skindepth0 = 0.01 +[radiation] + [radiation.drag] + [radiation.drag.synchrotron] + gamma_rad = 50.0 + + [radiation.drag.compton] + gamma_rad = 20.0 + + [radiation.emission] + [radiation.emission.synchrotron] + gamma_qed = 100.0 + photon_energy_min = 0.01 + photon_weight = 10.0 + photon_species = 3 + [algorithms] current_filters = 8 @@ -135,9 +150,6 @@ const auto sph_2d = u8R"( e_ovr_b_max = 0.95 larmor_max = 0.025 - [algorithms.synchrotron] - gamma_rad = 50.0 - [particles] ppc0 = 25.0 use_weights = true @@ -151,7 +163,8 @@ const auto sph_2d = u8R"( maxnpart = 1e2 pusher = "boris,gca" n_payloads_real = 3 - cooling = "synchrotron" + radiative_drag = "synchrotron,compton" + emission = "synchrotron" [[particles.species]] label = "e+" @@ -159,7 +172,7 @@ const auto sph_2d = u8R"( charge = 1.0 maxnpart = 1e2 pusher = "boris,gca" - cooling = "synchrotron" + radiative_drag = "synchrotron" n_payloads_int = 2 [[particles.species]] @@ -254,135 +267,135 @@ auto main(int argc, char* argv[]) -> int { try { using namespace ntt; - { - auto params_mink_1d = SimulationParams(); - params_mink_1d.setImmutableParams(mink_1d); - params_mink_1d.setMutableParams(mink_1d); - params_mink_1d.setSetupParams(mink_1d); - params_mink_1d.checkPromises(); - - assert_equal(params_mink_1d.get("grid.metric.metric"), - Metric::Minkowski, - "grid.metric.metric"); - // engine - assert_equal( - params_mink_1d.get("simulation.engine"), - SimEngine::SRPIC, - "simulation.engine"); - - assert_equal(params_mink_1d.get("scales.dx0"), - (real_t)0.0078125, - "scales.dx0"); - assert_equal(params_mink_1d.get("scales.V0"), - (real_t)0.0078125, - "scales.V0"); - boundaries_t fbc = { - { FldsBC::MATCH, FldsBC::MATCH } - }; - assert_equal( - params_mink_1d.get>("grid.boundaries.fields")[0].first, - fbc[0].first, - "grid.boundaries.fields[0].first"); - assert_equal( - params_mink_1d.get>("grid.boundaries.fields")[0].second, - fbc[0].second, - "grid.boundaries.fields[0].second"); - assert_equal( - params_mink_1d.get>("grid.boundaries.fields").size(), - fbc.size(), - "grid.boundaries.fields.size()"); - assert_equal( - params_mink_1d.get>("grid.boundaries.match.ds")[0].first, - (real_t)0.025, - "grid.boundaries.match.ds[0].first"); - assert_equal( - params_mink_1d.get>("grid.boundaries.match.ds")[0].second, - (real_t)0.1, - "grid.boundaries.match.ds[0].first"); - - const auto species = params_mink_1d.get>( - "particles.species"); - assert_equal(species[0].label(), "e-", "species[0].label"); - assert_equal(species[0].mass(), 1.0f, "species[0].mass"); - assert_equal(species[0].charge(), -1.0f, "species[0].charge"); - assert_equal(species[0].maxnpart(), 100, "species[0].maxnpart"); - assert_equal(species[0].pusher(), - PrtlPusher::BORIS, - "species[0].pusher"); - assert_equal(species[0].npld_r(), 3, "species[0].npld_r"); - assert_equal(species[0].npld_i(), 1, "species[0].npld_i"); - assert_equal(species[0].use_tracking(), true, "species[0].tracking"); - - assert_equal(species[1].label(), "p+", "species[1].label"); - assert_equal(species[1].mass(), 1.0f, "species[1].mass"); - assert_equal(species[1].charge(), 200.0f, "species[1].charge"); - assert_equal(species[1].maxnpart(), 100, "species[1].maxnpart"); - assert_equal(species[1].pusher(), - PrtlPusher::VAY, - "species[1].pusher"); - assert_equal(species[1].npld_r(), 0, "species[1].npld_r"); - - assert_equal(params_mink_1d.get("setup.myfloat"), - (real_t)(1e-2), - "setup.myfloat"); - assert_equal(params_mink_1d.get("setup.myint"), - (int)(123), - "setup.myint"); - assert_equal(params_mink_1d.get("setup.mybool"), - true, - "setup.mybool"); - const auto myarr = params_mink_1d.get>("setup.myarr"); - assert_equal(myarr[0], 1.0, "setup.myarr[0]"); - assert_equal(myarr[1], 2.0, "setup.myarr[1]"); - assert_equal(myarr[2], 3.0, "setup.myarr[2]"); - assert_equal(params_mink_1d.get("setup.mystr"), - "hi", - "setup.mystr"); - - const auto output_stride = params_mink_1d.get>( - "output.fields.downsampling"); - assert_equal(output_stride.size(), - 1, - "output.fields.downsampling.size()"); - assert_equal(output_stride[0], 4, "output.fields.downsampling[0]"); - - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.delta_x"), - (real_t)(1.0), - "algorithms.fieldsolver.delta_x"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.delta_y"), - (real_t)(2.0), - "algorithms.fieldsolver.delta_y"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.delta_z"), - (real_t)(3.0), - "algorithms.fieldsolver.delta_z"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.beta_xy"), - (real_t)(4.0), - "algorithms.fieldsolver.beta_xy"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.beta_yx"), - (real_t)(5.0), - "algorithms.fieldsolver.beta_yx"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.beta_xz"), - (real_t)(6.0), - "algorithms.fieldsolver.beta_xz"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.beta_zx"), - (real_t)(7.0), - "algorithms.fieldsolver.beta_zx"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.beta_yz"), - (real_t)(8.0), - "algorithms.fieldsolver.beta_yz"); - assert_equal( - params_mink_1d.get("algorithms.fieldsolver.beta_zy"), - (real_t)(9.0), - "algorithms.fieldsolver.beta_zy"); - } + // { + // auto params_mink_1d = SimulationParams(); + // params_mink_1d.setImmutableParams(mink_1d); + // params_mink_1d.setMutableParams(mink_1d); + // params_mink_1d.setSetupParams(mink_1d); + // params_mink_1d.checkPromises(); + // + // assert_equal(params_mink_1d.get("grid.metric.metric"), + // Metric::Minkowski, + // "grid.metric.metric"); + // // engine + // assert_equal( + // params_mink_1d.get("simulation.engine"), + // SimEngine::SRPIC, + // "simulation.engine"); + // + // assert_equal(params_mink_1d.get("scales.dx0"), + // (real_t)0.0078125, + // "scales.dx0"); + // assert_equal(params_mink_1d.get("scales.V0"), + // (real_t)0.0078125, + // "scales.V0"); + // boundaries_t fbc = { + // { FldsBC::MATCH, FldsBC::MATCH } + // }; + // assert_equal( + // params_mink_1d.get>("grid.boundaries.fields")[0].first, + // fbc[0].first, + // "grid.boundaries.fields[0].first"); + // assert_equal( + // params_mink_1d.get>("grid.boundaries.fields")[0].second, + // fbc[0].second, + // "grid.boundaries.fields[0].second"); + // assert_equal( + // params_mink_1d.get>("grid.boundaries.fields").size(), + // fbc.size(), + // "grid.boundaries.fields.size()"); + // assert_equal( + // params_mink_1d.get>("grid.boundaries.match.ds")[0].first, + // (real_t)0.025, + // "grid.boundaries.match.ds[0].first"); + // assert_equal( + // params_mink_1d.get>("grid.boundaries.match.ds")[0].second, + // (real_t)0.1, + // "grid.boundaries.match.ds[0].first"); + // + // const auto species = params_mink_1d.get>( + // "particles.species"); + // assert_equal(species[0].label(), "e-", "species[0].label"); + // assert_equal(species[0].mass(), 1.0f, "species[0].mass"); + // assert_equal(species[0].charge(), -1.0f, "species[0].charge"); + // assert_equal(species[0].maxnpart(), 100, "species[0].maxnpart"); + // assert_equal(species[0].pusher(), + // ParticlePusher::BORIS, + // "species[0].pusher"); + // assert_equal(species[0].npld_r(), 3, "species[0].npld_r"); + // assert_equal(species[0].npld_i(), 1, "species[0].npld_i"); + // assert_equal(species[0].use_tracking(), true, "species[0].tracking"); + // + // assert_equal(species[1].label(), "p+", "species[1].label"); + // assert_equal(species[1].mass(), 1.0f, "species[1].mass"); + // assert_equal(species[1].charge(), 200.0f, "species[1].charge"); + // assert_equal(species[1].maxnpart(), 100, "species[1].maxnpart"); + // assert_equal(species[1].pusher(), + // ParticlePusher::VAY, + // "species[1].pusher"); + // assert_equal(species[1].npld_r(), 0, "species[1].npld_r"); + // + // assert_equal(params_mink_1d.get("setup.myfloat"), + // (real_t)(1e-2), + // "setup.myfloat"); + // assert_equal(params_mink_1d.get("setup.myint"), + // (int)(123), + // "setup.myint"); + // assert_equal(params_mink_1d.get("setup.mybool"), + // true, + // "setup.mybool"); + // const auto myarr = params_mink_1d.get>("setup.myarr"); + // assert_equal(myarr[0], 1.0, "setup.myarr[0]"); + // assert_equal(myarr[1], 2.0, "setup.myarr[1]"); + // assert_equal(myarr[2], 3.0, "setup.myarr[2]"); + // assert_equal(params_mink_1d.get("setup.mystr"), + // "hi", + // "setup.mystr"); + // + // const auto output_stride = params_mink_1d.get>( + // "output.fields.downsampling"); + // assert_equal(output_stride.size(), + // 1, + // "output.fields.downsampling.size()"); + // assert_equal(output_stride[0], 4, "output.fields.downsampling[0]"); + // + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.delta_x"), + // (real_t)(1.0), + // "algorithms.fieldsolver.delta_x"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.delta_y"), + // (real_t)(2.0), + // "algorithms.fieldsolver.delta_y"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.delta_z"), + // (real_t)(3.0), + // "algorithms.fieldsolver.delta_z"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.beta_xy"), + // (real_t)(4.0), + // "algorithms.fieldsolver.beta_xy"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.beta_yx"), + // (real_t)(5.0), + // "algorithms.fieldsolver.beta_yx"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.beta_xz"), + // (real_t)(6.0), + // "algorithms.fieldsolver.beta_xz"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.beta_zx"), + // (real_t)(7.0), + // "algorithms.fieldsolver.beta_zx"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.beta_yz"), + // (real_t)(8.0), + // "algorithms.fieldsolver.beta_yz"); + // assert_equal( + // params_mink_1d.get("algorithms.fieldsolver.beta_zy"), + // (real_t)(9.0), + // "algorithms.fieldsolver.beta_zy"); + // } { auto params_sph_2d = SimulationParams(); @@ -454,9 +467,31 @@ auto main(int argc, char* argv[]) -> int { "algorithms.gca.larmor_max"); assert_equal( - params_sph_2d.get("algorithms.synchrotron.gamma_rad"), + params_sph_2d.get("radiation.drag.synchrotron.gamma_rad"), (real_t)50.0, - "algorithms.synchrotron.gamma_rad"); + "radiation.drag.synchrotron.gamma_rad"); + + assert_equal( + params_sph_2d.get("radiation.drag.compton.gamma_rad"), + (real_t)20.0, + "radiation.drag.compton.gamma_rad"); + + assert_equal(params_sph_2d.get( + "radiation.emission.synchrotron.photon_energy_min"), + (real_t)0.01, + "radiation.emission.synchrotron.photon_energy_min"); + assert_equal( + params_sph_2d.get("radiation.emission.synchrotron.gamma_qed"), + (real_t)100.0, + "radiation.emission.synchrotron.gamma_qed"); + assert_equal(params_sph_2d.get( + "radiation.emission.synchrotron.photon_weight"), + (real_t)10.0, + "radiation.emission.synchrotron.photon_weight"); + assert_equal(params_sph_2d.get( + "radiation.emission.synchrotron.photon_species"), + 3, + "radiation.emission.synchrotron.photon_species"); const auto species = params_sph_2d.get>( "particles.species"); @@ -464,27 +499,30 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[0].mass(), 1.0f, "species[0].mass"); assert_equal(species[0].charge(), -1.0f, "species[0].charge"); assert_equal(species[0].maxnpart(), 100, "species[0].maxnpart"); - assert_equal(species[0].pusher(), - PrtlPusher::BORIS, - "species[0].pusher"); - assert_equal(species[0].use_gca(), true, "species[0].use_gca"); + assert_equal(species[0].pusher(), + ParticlePusher::BORIS | ParticlePusher::GCA, + "species[0].pusher"); assert_equal(species[0].npld_r(), 3, "species[0].npld_r"); - assert_equal(species[0].cooling(), - Cooling::SYNCHROTRON, - "species[0].cooling"); + assert_equal(species[0].radiative_drag_flags(), + RadiativeDrag::SYNCHROTRON | + RadiativeDrag::COMPTON, + "species[0].radiative_drag_flags"); + + assert_equal(species[0].emission_policy_flag(), + EmissionType::SYNCHROTRON, + "species[0].emission_policy_flag"); assert_equal(species[1].label(), "e+", "species[1].label"); assert_equal(species[1].mass(), 1.0f, "species[1].mass"); assert_equal(species[1].charge(), 1.0f, "species[1].charge"); assert_equal(species[1].maxnpart(), 100, "species[1].maxnpart"); - assert_equal(species[1].pusher(), - PrtlPusher::BORIS, - "species[1].pusher"); - assert_equal(species[1].use_gca(), true, "species[1].use_gca"); + assert_equal(species[1].pusher(), + ParticlePusher::BORIS | ParticlePusher::GCA, + "species[1].pusher"); assert_equal(species[1].npld_r(), 0, "species[1].npld_r"); - assert_equal(species[1].cooling(), - Cooling::SYNCHROTRON, - "species[1].cooling"); + assert_equal(species[1].radiative_drag_flags(), + RadiativeDrag::SYNCHROTRON, + "species[1].radiative_drag_flags"); assert_equal(species[1].npld_i(), 2, "species[1].npld_i"); assert_equal(species[1].use_tracking(), false, "species[1].tracking"); @@ -492,9 +530,9 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[2].mass(), 0.0f, "species[2].mass"); assert_equal(species[2].charge(), 0.0f, "species[2].charge"); assert_equal(species[2].maxnpart(), 100, "species[2].maxnpart"); - assert_equal(species[2].pusher(), - PrtlPusher::PHOTON, - "species[2].pusher"); + assert_equal(species[2].pusher(), + ParticlePusher::PHOTON, + "species[2].pusher"); assert_equal(species[2].npld_r(), 0, "species[2].npld_r"); } @@ -601,18 +639,18 @@ auto main(int argc, char* argv[]) -> int { assert_equal(species[0].mass(), 1.0f, "species[0].mass"); assert_equal(species[0].charge(), -1.0f, "species[0].charge"); assert_equal(species[0].maxnpart(), 100, "species[0].maxnpart"); - assert_equal(species[0].pusher(), - PrtlPusher::BORIS, - "species[0].pusher"); + assert_equal(species[0].pusher(), + ParticlePusher::BORIS, + "species[0].pusher"); assert_equal(species[0].npld_r(), 0, "species[0].npld_r"); assert_equal(species[1].label(), "e+", "species[1].label"); assert_equal(species[1].mass(), 1.0f, "species[1].mass"); assert_equal(species[1].charge(), 1.0f, "species[1].charge"); assert_equal(species[1].maxnpart(), 100, "species[1].maxnpart"); - assert_equal(species[1].pusher(), - PrtlPusher::BORIS, - "species[1].pusher"); + assert_equal(species[1].pusher(), + ParticlePusher::BORIS, + "species[1].pusher"); assert_equal(species[1].npld_r(), 0, "species[1].npld_r"); } diff --git a/src/framework/tests/particles.cpp b/src/framework/tests/particles.cpp index ab3aa0e0..cce6af28 100644 --- a/src/framework/tests/particles.cpp +++ b/src/framework/tests/particles.cpp @@ -9,16 +9,18 @@ #include template -void testParticles(int index, - const std::string& label, - float m, - float ch, - std::size_t maxnpart, - const ntt::PrtlPusher& pusher, - bool use_tracking, - const ntt::Cooling& cooling, - unsigned short npld_r = 0, - unsigned short npld_i = 0) { +void testParticles( + int index, + const std::string& label, + float m, + float ch, + std::size_t maxnpart, + ntt::ParticlePusherFlags pusher, + bool use_tracking, + ntt::RadiativeDragFlags radiative_drag_flags, + ntt::EmissionTypeFlag emission_policy_flag = ntt::EmissionType::NONE, + unsigned short npld_r = 0, + unsigned short npld_i = 0) { using namespace ntt; auto p = Particles(index, label, @@ -27,8 +29,8 @@ void testParticles(int index, maxnpart, pusher, use_tracking, - false, - cooling, + radiative_drag_flags, + emission_policy_flag, npld_r, npld_i); raise::ErrorIf(p.index() != index, "Index mismatch", HERE); @@ -37,7 +39,12 @@ void testParticles(int index, raise::ErrorIf(p.charge() != ch, "Charge mismatch", HERE); raise::ErrorIf(p.maxnpart() != maxnpart, "Max number of particles mismatch", HERE); raise::ErrorIf(p.pusher() != pusher, "Pusher mismatch", HERE); - raise::ErrorIf(p.cooling() != cooling, "Cooling mismatch", HERE); + raise::ErrorIf(p.radiative_drag_flags() != radiative_drag_flags, + "Radiative drag mismatch", + HERE); + raise::ErrorIf(p.emission_policy_flag() != emission_policy_flag, + "Emission policy mismatch", + HERE); raise::ErrorIf(p.npart() != 0, "Number of particles mismatch", HERE); raise::ErrorIf(p.i1.extent(0) != maxnpart, "i1 incorrectly allocated", HERE); @@ -115,17 +122,19 @@ auto main(int argc, char** argv) -> int { 1.0, -1.0, 100, - PrtlPusher::BORIS, + ParticlePusher::BORIS, false, - Cooling::SYNCHROTRON); + RadiativeDrag::SYNCHROTRON); testParticles(2, "p+", 100.0, -1.0, 1000, - PrtlPusher::VAY, + ParticlePusher::VAY, true, - Cooling::SYNCHROTRON, + RadiativeDrag::SYNCHROTRON | + RadiativeDrag::COMPTON, + EmissionType::SYNCHROTRON, 2, 1); testParticles(3, @@ -133,18 +142,20 @@ auto main(int argc, char** argv) -> int { 0.0, 0.0, 100, - PrtlPusher::PHOTON, + ParticlePusher::PHOTON, false, - Cooling::NONE, + RadiativeDrag::NONE, + EmissionType::COMPTON, 5); testParticles(4, "e+", 1.0, 1.0, 100, - PrtlPusher::BORIS, + ParticlePusher::BORIS, true, - Cooling::NONE, + RadiativeDrag::NONE, + EmissionType::NONE, 2, 3); testParticles(5, @@ -152,9 +163,10 @@ auto main(int argc, char** argv) -> int { 1.0, 1.0, 100, - PrtlPusher::BORIS, + ParticlePusher::BORIS, false, - Cooling::NONE, + RadiativeDrag::NONE, + EmissionType::NONE, 1, 2); } catch (const std::exception& e) { diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h index e950af53..88a81a51 100644 --- a/src/global/arch/traits.h +++ b/src/global/arch/traits.h @@ -2,6 +2,9 @@ * @file arch/traits.h * @brief Defines a set of traits to check if a class satisfies certain conditions * @implements + * - traits::fieldsetter::ex1_t, ::ex2_t, ::ex3_t - checks for special methods in field setter classes + * - traits::fieldsetter::bx1_t, ::bx2_t, ::bx3_t - checks for special methods in field setter classes + * - traits::fieldsetter::dx1_t, ::dx2_t, ::dx3_t - checks for special methods in field setter classes * - traits::has_method<> * - traits::has_member<> * - traits::ex1_t, ::ex2_t, ::ex3_t @@ -22,78 +25,82 @@ #include "global.h" -#include "arch/kokkos_aliases.h" - +#include +#include #include #include namespace traits { - template