diff --git a/CHANGELOG.md b/CHANGELOG.md index b0732bb5..3778d514 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,20 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [2.2.02] 2025-10-18 +### Changed + +- Added a module to add explicit units (#338) +- fixed a bug that could lead to incorrect profiling information on non-blocking cuda loops (#341) +- fixed a bug that could lead to incorrect energy budget when shearing box and fargo were both enabled (#346) +- fixed a bug that led to incorrect BX2 reconstruction when axis is not used on both sides of the domain (#345) +- fixed a bug that led to incorrect reflective boundary conditions on B when DIMENSIONS < 3 (#345) +- fixed a bug that led to incorrect dust stopping time when the adiabatic equation of state is used with "size" drag law (#353) + +### Added + +- documentation for the continuous integration (#354) + ## [2.2.01] 2025-04-16 ### Changed diff --git a/CMakeLists.txt b/CMakeLists.txt index 6b62b754..364ee4ef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,9 +6,9 @@ set (CMAKE_CXX_STANDARD 17) set(Idefix_VERSION_MAJOR 2) set(Idefix_VERSION_MINOR 2) -set(Idefix_VERSION_PATCH 01) +set(Idefix_VERSION_PATCH 02) -project (idefix VERSION 2.2.01) +project (idefix VERSION 2.2.02) option(Idefix_MHD "enable MHD" OFF) option(Idefix_MPI "enable Message Passing Interface parallelisation" OFF) option(Idefix_HIGH_ORDER_FARGO "Force Fargo to use a PPM reconstruction scheme" OFF) diff --git a/doc/source/conf.py b/doc/source/conf.py index 6f0cffa4..595672c2 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -23,7 +23,7 @@ author = 'Geoffroy Lesur' # The full version, including alpha/beta/rc tags -release = '2.2.00' +release = '2.2.02' diff --git a/doc/source/faq.rst b/doc/source/faq.rst index 4c7bc056..4dd19831 100644 --- a/doc/source/faq.rst +++ b/doc/source/faq.rst @@ -61,9 +61,18 @@ How can I stop the code without loosing the current calculation? I'm doing performance measures. How do I disable all outputs in *Idefix*? Add ``-nowrite`` when you call *Idefix* executable. +VTK output appears corrupted when running with MPI (OpenMPI) + Some OpenMPI configurations (notably OpenMPI 4 with the ompio component) can produce corrupted VTK/VTU output when running with MPI enabled. This appears to be caused by bugs in OpenMPI's ompio I/O component. + Disable ompio so OpenMPI falls back to ROMIO (MPICH's MPI-IO), which is typically more stable: -Developement ------------- + .. code-block:: console + + mpirun --mca io ^ompio ... + + This has resolved intermittent corruption for several users. See issue #348 for discussion and reports. + +Development +----------- I have a serious bug (e.g. segmentation fault), in my setup, how do I proceed? Add ``-DIdefix_DEBUG=ON`` to ``cmake`` and recompile to find out exactly where the code crashes (see :ref:`debugging`). diff --git a/doc/source/index.rst b/doc/source/index.rst index 411d73d6..d36ab39b 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -122,6 +122,7 @@ The Idefix collaboration benefited from funding from the “Programme National d reference modules programmingguide + testing performances kokkos contributing diff --git a/doc/source/reference/idefix.ini.rst b/doc/source/reference/idefix.ini.rst index 1fa05cc6..8b558722 100644 --- a/doc/source/reference/idefix.ini.rst +++ b/doc/source/reference/idefix.ini.rst @@ -265,8 +265,10 @@ section as followed: +----------------+-------------------------+---------------------------------------------------------------------------------------------+ | Mcentral | real | | Mass of the central object when a central potential is enabled (see above). Default is 1. | +----------------+-------------------------+---------------------------------------------------------------------------------------------+ -| gravCst | real | | Set the value of the gravitational constant :math:`G_c` used by the central | +| gravCst | real or string | | Set the value of the gravitational constant :math:`G_c` used by the central | | | | | mass potential and self-gravitational potential (when enabled) ). Default is 1. | +| | | | If set to "units", Idefix will compute the gravitational constant from the user units | +| | | | defined in the units section (see :ref:`unitsSection`). | +----------------+-------------------------+---------------------------------------------------------------------------------------------+ | bodyForce | string | | Adds an acceleration vector to each cell of the domain. The only value allowed | | | | | is ``userdef``. The ``Gravity`` class then expects a user-defined bodyforce function to | @@ -310,7 +312,28 @@ of gravitational potential in the ``Gravity`` section (see :ref:`gravitySection` | | | | Default is 1 (i.e. self-gravity is computed at every cycle). | +----------------+-------------------------+---------------------------------------------------------------------------------------------+ +.. _unitsSection: +``Units`` section +------------------ + +This optional section controls how the code handles units. By default, Idefix work "unit free" and this section is not useful. +However, some physical processes like self-gravity or radiative transfer require the explicit use of units to define natural constants. +These entries define how one code units should be converted to CGS units. When used, the following values should be set: + ++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ +| Entry name | Parameter type | Comment | ++================+====================+===========================================================================================================+ +| length | float | The code unit length: 1 code unit = `length` cm. 1 by default. | ++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ +| velocity | float | The code unit velocity: 1 code unit = `velocity` cm/s. 1 by default. | ++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ +| density | float | The code unit density: 1 code unit = `density` g/cm^3. 1 by default. | ++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ + +.. note:: + The units module automatically reconstruct all of the units (time, temperature, magnetic field, etc.) from the above 3 fundamental constants. These can + all be obtained from the global object ``idfx::units``. Note however that the code outputs remain in code units, even if `[Units]` are defined. ``RKL`` section ------------------ diff --git a/doc/source/testing.rst b/doc/source/testing.rst new file mode 100644 index 00000000..888ac205 --- /dev/null +++ b/doc/source/testing.rst @@ -0,0 +1,117 @@ +Continuous Integration (CI) tests +================================ + +This document describes the GitHub Actions continuous-integration setup used to run the Idefix +test-suite. The CI is implemented by two workflows checked in .github/workflows: + +- .github/workflows/idefix-ci.yml +- .github/workflows/idefix-ci-jobs.yml + +Overview +-------- + +The CI is split in two layers: + +- A top-level workflow (.github/workflows/idefix-ci.yml) that: + + - runs a Linter job (pre-commit) on push / PR / manual dispatch, + - then calls a reusable workflow for different compiler/backends (intel, gcc, cuda) + providing two inputs: TESTME_OPTIONS and IDEFIX_COMPILER. + +- A reusable workflow (.github/workflows/idefix-ci-jobs.yml) that: + + - defines the actual test jobs grouped by physics domain (ShocksHydro, ParabolicHydro, + ShocksMHD, ParabolicMHD, Fargo, ShearingBox, SelfGravity, Planet, Dust, Braginskii, + Examples, Utils), + - runs test scripts on self-hosted runners, + - expects the repository to be checked out with submodules, + - invokes the repository-provided CI helper scripts to configure / build / run tests. + +Key configuration points +------------------------ + +- Inputs passed from the top-level workflow: + + - TESTME_OPTIONS (string): flags forwarded to the per-test runner (examples: -cuda, -Werror, + -intel, -all). + - IDEFIX_COMPILER (string): which compiler the tests should use (e.g. icc, gcc, nvcc). + +- Environment variables set by the reusable workflow: + + - IDEFIX_COMPILER, TESTME_OPTIONS, PYTHONPATH, IDEFIX_DIR + +- Linter job: + + - Runs only when repository is the main project (not arbitrary forks). + - Uses actions/setup-python and runs pre-commit (pre-commit/action@v3 and pre-commit-ci/lite). + - Prevents regressions in style and common mistakes before running heavy test jobs. + +- Test execution: + + - All test jobs call the repository script scripts/ci/run-tests with a test directory + and the TESTME_OPTIONS flags. Example invocation (from the workflows): + scripts/ci/run-tests $IDEFIX_DIR/test/HD/sod -all $TESTME_OPTIONS + + - The reusable workflow is written to execute many test directories in separate job steps, + so each physics group is kept logically separated in CI logs. + +Runners and prerequisites +------------------------- + +- The heavy numerical tests run on self-hosted runners (see runs-on: self-hosted). + The CI assumes appropriate hardware and dependencies are available on those runners + (compilers, MPI, GPUs when CUDA/HIP flags are used, required system libraries). + +- The workflows check out the repository and its submodules. Submodules must be available + on the CI machines. + +How tests are driven (testme scripts) +------------------------------------- + +Each test directory contains a small Python "testMe" driver that uses the helper Python +class documented in the repository: + +- See the test helper documentation: :doc:`idfxTest ` + +That helper (idfxTest) is responsible for: + +- parsing TESTME_OPTIONS-like flags (precision, MPI, CUDA, reconstruction, vector potential, etc.), +- calling configure / compile / run, +- performing standard python checks and non-regression (RMSE) comparisons against + reference dumps, +- optionally creating / updating reference dumps (init mode). + +Practical examples +------------------ + +- Example of a CI invocation (triggered by workflows): + + - Top-level workflow calls the reusable jobs workflow for each compiler/back-end, e.g. + TESTME_OPTIONS="-cuda -Werror" IDEFIX_COMPILER=nvcc + +- Running tests locally (developer machine) + - You can mimic what CI does by calling the repository helper script directly. Example: + scripts/ci/run-tests /path/to/idefix/test/HD/sod -all -mpi -dec 2 2 -reconstruction 3 -single + +Notes for maintainers +--------------------- + +- The reusable jobs workflow contains a commented concurrency block for optional cancellation + of in-flight runs — consider enabling it if you want to auto-cancel redundant CI runs. +- Because tests are run on self-hosted runners, ensure the pools have the required compilers, + MPI stacks and GPU drivers for the requested TESTME_OPTIONS. +- Keep TESTME_OPTIONS in sync with the options understood by the test helper documented in + :doc:`idfxTest `. + +Relevant files +-------------- + +- Workflow entry point: .github/workflows/idefix-ci.yml +- Reusable jobs: .github/workflows/idefix-ci-jobs.yml +- Test helper documentation: :doc:`idfxTest ` + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + testing/idfxTest.rst diff --git a/doc/source/testing/idfxTest.rst b/doc/source/testing/idfxTest.rst new file mode 100644 index 00000000..b425d3fe --- /dev/null +++ b/doc/source/testing/idfxTest.rst @@ -0,0 +1,172 @@ +========= +idfxTest +========= + +.. autoclass:: idfxTest + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The ``idfxTest`` class provides a high-level interface for automating the configuration, compilation, execution, and regression testing of Idefix simulations. It is designed to be used in test scripts (such as ``testme.py``) to streamline the testing workflow, including handling reference files and plotting differences. + +Constructor and Command-Line Options +------------------------------------ + +The constructor parses command-line arguments using ``argparse``. These options can be passed directly to the test script or via the command line. The following options are available: + +.. list-table:: + :header-rows: 1 + + * - Option + - Attribute + - Description + * - ``-noplot`` + - ``noplot`` + - Disable plotting in standard tests (default: True). + * - ``-ploterr`` + - ``ploterr`` + - Enable plotting of differences when regression tests fail. + * - ``-cmake OPT [OPT ...]`` + - ``cmake`` + - Extra CMake options (list of strings). + * - ``-definitions FILE`` + - ``definitions`` + - Specify a custom ``definitions.hpp`` file. + * - ``-dec NX NY NZ`` + - ``dec`` + - MPI domain decomposition (list of integers). + * - ``-check`` + - ``check`` + - Only perform regression tests without compilation. + * - ``-cuda`` + - ``cuda`` + - Enable CUDA backend for Nvidia GPUs. + * - ``-intel`` + - ``intel`` + - Use Intel OneAPI compilers. + * - ``-hip`` + - ``hip`` + - Enable HIP backend for AMD GPUs. + * - ``-single`` + - ``single`` + - Enable single precision. + * - ``-vectPot`` + - ``vectPot`` + - Enable vector potential formulation. + * - ``-reconstruction N`` + - ``reconstruction`` + - Set reconstruction scheme (2=PLM, 3=LimO3, 4=PPM). + * - ``-idefixDir PATH`` + - ``idefixDir`` + - Set the directory for Idefix source files (default: ``$IDEFIX_DIR``). + * - ``-mpi`` + - ``mpi`` + - Enable MPI parallelism. + * - ``-all`` + - ``all`` + - Run the full test suite with multiple configurations. + * - ``-init`` + - ``init`` + - Reinitialize reference files for non-regression tests. + * - ``-Werror`` + - ``Werror`` + - Treat compiler warnings as errors. + +Main Methods +------------ + +.. list-table:: + :header-rows: 1 + + * - Method + - Description + * - ``configure`` + - Runs CMake to configure the build system for Idefix, using options set by command-line flags (e.g., precision, MPI, CUDA, etc.). + * - ``compile`` + - Compiles the Idefix code using ``make`` with the specified number of parallel jobs. + * - ``run`` + - Executes the Idefix binary, optionally with MPI, using the provided input file and runtime options. + * - ``checkOnly`` + - Performs regression testing only, without compiling or running the code (useful for checking outputs after a manual run). + * - ``standardTest`` + - Runs any Python-based standard tests (e.g., ``testidefix.py``) present in the test directory for additional validation. + * - ``nonRegressionTest`` + - Compares the output dump file to a reference file using RMSE; fails if the error exceeds the tolerance. + * - ``compareDump`` + - Compares two arbitrary dump files using the same logic as ``nonRegressionTest``. + * - ``makeReference`` + - Copies the specified output file to the reference directory, updating the reference for future regression tests. + +Usage Example +------------- + +Below is an example inspired by ``testme.py`` from ``test/HD/sod/testme.py``. This demonstrates a typical workflow for running tests and performing regression checks. + +.. code-block:: python + + import pytools.idfx_test as tst + + name = "dump.0001.dmp" + + def testMe(test): + test.configure() + test.compile() + inifiles = ["idefix.ini", "idefix-hll.ini", "idefix-hllc.ini", "idefix-tvdlf.ini"] + if test.reconstruction == 4: + inifiles = ["idefix-rk3.ini", "idefix-hllc-rk3.ini"] + + # Loop over all ini files for this test + for ini in inifiles: + test.run(inputFile=ini) + if test.init: + test.makeReference(filename=name) + test.standardTest() + test.nonRegressionTest(filename=name) + + test = tst.idfxTest() + + if not test.all: + if test.check: + test.checkOnly(filename=name) + else: + testMe(test) + else: + test.noplot = True + for rec in range(2, 5): + test.vectPot = False + test.single = False + test.reconstruction = rec + test.mpi = False + testMe(test) + + # Test in single precision + test.reconstruction = 2 + test.single = True + testMe(test) + +How to Run +---------- + +You can run the test script from the command line, passing any of the supported options. For example: + +.. code-block:: bash + + python testme.py -mpi -dec 2 2 -reconstruction 3 -single -ploterr -idefixDir /path/to/idefix + +This will configure, compile, and run the test in MPI mode with a 2x2 domain decomposition, third-order reconstruction, single precision, and plotting enabled for regression errors. + +Reference File Management +------------------------ + +- Reference files are stored in ``$IDEFIX_DIR/reference/``. +- The filename is generated based on precision, reconstruction, input file, and vector potential settings. +- Use ``test.init`` to regenerate reference files (dangerous: overwrites existing references). + +Regression Testing +------------------ + +- The ``nonRegressionTest`` method compares output dumps to reference files using RMSE. +- If the error exceeds the tolerance, the test fails and (optionally) plots the difference. diff --git a/reference b/reference index c3289420..c4082b99 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit c328942083b618a85b84eb16ce9fd35e67c00597 +Subproject commit c4082b99a4c542def3177c96cb35b1c9d9002f18 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 727f2150..08dd93c6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -28,6 +28,8 @@ target_sources(idefix PUBLIC ${CMAKE_CURRENT_LIST_DIR}/reduce.hpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/setup.hpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/setup.cpp + PUBLIC ${CMAKE_CURRENT_LIST_DIR}/units.hpp + PUBLIC ${CMAKE_CURRENT_LIST_DIR}/units.cpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/timeIntegrator.hpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/timeIntegrator.cpp ) diff --git a/src/fluid/addSourceTerms.hpp b/src/fluid/addSourceTerms.hpp index a0cb5a71..e3f74586 100644 --- a/src/fluid/addSourceTerms.hpp +++ b/src/fluid/addSourceTerms.hpp @@ -83,9 +83,6 @@ struct Fluid_AddSourceTermsFunctor { Uc(MX1,k,j,i) += TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * Vc(VX2,k,j,i); Uc(MX2,k,j,i) += - TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * Vc(VX1,k,j,i); } - if(haveFargo) { - Uc(MX1,k,j,i) += TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * sbS * x1(i); - } #endif // fetch fargo velocity when required [[maybe_unused]] real fargoV = ZERO_F; diff --git a/src/fluid/boundary/axis.cpp b/src/fluid/boundary/axis.cpp index 3ad67f0a..e447aefc 100644 --- a/src/fluid/boundary/axis.cpp +++ b/src/fluid/boundary/axis.cpp @@ -358,65 +358,28 @@ void Axis::EnforceAxisBoundary(int side) { // Reconstruct Bx2s taking care of the sides where an axis is lying -void Axis::ReconstructBx2s() { - idfx::pushRegion("Axis::ReconstructBx2s"); +void Axis::RegularizeBX2s() { + idfx::pushRegion("Axis::RegularizeBX2s"); #if DIMENSIONS >= 2 && MHD == YES IdefixArray4D Vs = this->Vs; - IdefixArray3D Ax1=data->A[IDIR]; - IdefixArray3D Ax2=data->A[JDIR]; - IdefixArray3D Ax3=data->A[KDIR]; - int nstart = data->beg[JDIR]-1; - int nend = data->end[JDIR]; - int ntot = data->np_tot[JDIR]; - - bool haveleft = axisLeft; - bool haveright = axisRight; - - if(haveleft) { - // This loop is a copy of ReconstructNormalField, with the proper sign when we cross the axis - idefix_for("Axis::ReconstructBX2sLeft",0,data->np_tot[KDIR],0,data->np_tot[IDIR], - KOKKOS_LAMBDA (int k, int i) { - for(int j = nstart ; j>=0 ; j-- ) { - Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i) - +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i) , - , - - ( Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - - Ax3(k,j,i) * Vs(BX3s,k,j,i) )))); - } - } - ); - } - if(haveright) { - // This loop is a copy of ReconstructNormalField, with the proper sign when we cross the axis - idefix_for("Axis::ReconstructBX2sRight",0,data->np_tot[KDIR],0,data->np_tot[IDIR], - KOKKOS_LAMBDA (int k, int i) { - for(int j = nend ; jend[JDIR]; int jleft = data->beg[JDIR]; if(isTwoPi) { #ifdef AXIS_BX2S_USE_ATHENA_REGULARISATION - if(haveleft) FixBx2sAxisGhostAverage(left); - if(haveright) FixBx2sAxisGhostAverage(right); + if(axisLeft) FixBx2sAxisGhostAverage(left); + if(axisRight) FixBx2sAxisGhostAverage(right); #else - if(haveleft) FixBx2sAxis(left); - if(haveright) FixBx2sAxis(right); + if(axisLeft) FixBx2sAxis(left); + if(axisRight) FixBx2sAxis(right); #endif } else { + bool haveleft = axisLeft; + bool haveright = axisRight; idefix_for("Axis:BoundaryAvg",0,data->np_tot[KDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA (int k, int i) { if(haveleft) { diff --git a/src/fluid/boundary/axis.hpp b/src/fluid/boundary/axis.hpp index a1c3054b..b7469786 100644 --- a/src/fluid/boundary/axis.hpp +++ b/src/fluid/boundary/axis.hpp @@ -29,7 +29,7 @@ class Axis { void RegularizeEMFs(); // Regularize the EMF sitting on the axis void RegularizeCurrent(); // Regularize the currents along the axis void EnforceAxisBoundary(int side); // Enforce the boundary conditions (along X2) - void ReconstructBx2s(); // Reconstruct BX2s in the ghost zone using divB=0 + void RegularizeBX2s(); // Regularize BX2s on the axis void ShowConfig(); @@ -42,6 +42,8 @@ class Axis { void ExchangeMPI(int side); // Function has to be public for GPU, but its technically // a private function + bool haveLeftAxis() const { return axisLeft; } ///< Check if the axis is on the left + bool haveRightAxis() const { return axisRight; } ///< Check if the axis is on the right private: bool isTwoPi = false; @@ -147,7 +149,6 @@ Axis::Axis(Boundary *boundary) { Kokkos::deep_copy(symmetryVc, symmetryVcHost); if constexpr(Phys::mhd) { - idfx::cout << "Phys MHD" << std::endl; symmetryVs = IdefixArray1D("Axis:SymmetryVs",DIMENSIONS); IdefixArray1D::HostMirror symmetryVsHost = Kokkos::create_mirror_view(symmetryVs); // Init the array diff --git a/src/fluid/boundary/boundary.hpp b/src/fluid/boundary/boundary.hpp index 86909936..897c1ce1 100644 --- a/src/fluid/boundary/boundary.hpp +++ b/src/fluid/boundary/boundary.hpp @@ -112,6 +112,8 @@ class Boundary { IdefixArray4D Vs; ///< reference to face-centered array that we should sync std::unique_ptr axis; ///< Axis object, initialised if needed. bool haveAxis{false}; + bool haveLeftAxis{false}; ///< True if the left boundary is an axis + bool haveRightAxis{false}; ///< True if the right boundary is an axis private: friend class Axis; @@ -137,6 +139,8 @@ Boundary::Boundary(Fluid* fluid) { if(data->haveAxis) { this->axis = std::make_unique(this); this->haveAxis = true; + this->haveLeftAxis = axis->haveLeftAxis(); + this->haveRightAxis = axis->haveRightAxis(); } @@ -464,30 +468,34 @@ void Boundary::ReconstructNormalField(int dir) { if(dir==JDIR) { nstart = data->beg[JDIR]-1; nend = data->end[JDIR]; - if(!this->haveAxis) { - idefix_for("ReconstructBX2s",0,data->np_tot[KDIR],0,data->np_tot[IDIR], - KOKKOS_LAMBDA (int k, int i) { - if(reconstructLeft) { - for(int j = nstart ; j>=0 ; j-- ) { - Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i) - +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i) , - , - + Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) ))); - } + #if DIMENSIONS == 3 + const int signLeft = haveLeftAxis ? -1 : 1; // left axis is in the negative direction + const int signRight = haveRightAxis ? -1 : 1; // right axis is in the negative direction + #endif + + idefix_for("ReconstructBX2s",0,data->np_tot[KDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int i) { + if(reconstructLeft) { + for(int j = nstart ; j>=0 ; j-- ) { + Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i) + +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i) , + , + + signLeft*(Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) )))); } - if(reconstructRight) { - for(int j = nend ; jReconstructBx2s(); + axis->RegularizeBX2s(); } } #endif @@ -686,7 +694,7 @@ void Boundary::EnforceReflective(int dir, BoundarySide side ) { const int jref = (dir==JDIR) ? 2*(jghost + side*nxj) - j - 1 : j; const int kref = (dir==KDIR) ? 2*(kghost + side*nxk) - k - 1 : k; - const int sign = (n == VX1+dir) ? -1.0 : 1.0; + const int sign = (n == VX1+dir || (n >= BX1 && n != BX1+dir)) ? -1.0 : 1.0; Vc(n,k,j,i) = sign * Vc(n,kref,jref,iref); }); diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp index 38ef76d5..7777d1c2 100644 --- a/src/fluid/calcRightHandSide.hpp +++ b/src/fluid/calcRightHandSide.hpp @@ -472,12 +472,18 @@ struct Fluid_CalcRHSFunctor { // Body force if(needBodyForce) { - rhs[MX1+dir] += dt * Vc(RHO,k,j,i) * bodyForce(dir,k,j,i); + real bf = bodyForce(dir,k,j,i); + #if GEOMETRY == CARTESIAN + if(haveFargo && dir==IDIR) { + // Remove Body force that is already compensated by Fargo shear*Rotation + bf -= -2*Omega*sbS * x1(i); + } + #endif + rhs[MX1+dir] += dt * Vc(RHO,k,j,i) * bf; if constexpr(Phys::pressure) { // rho * v . f, where rhov is taken as a volume average of Flux(RHO) rhs[ENG] += HALF_F * dtdV * dl * - (Flux(RHO,k,j,i) + Flux(RHO, k+koffset, j+joffset, i+ioffset)) * - bodyForce(dir,k,j,i); + (Flux(RHO,k,j,i) + Flux(RHO, k+koffset, j+joffset, i+ioffset)) * bf; } // Pressure // Particular cases if we do not sweep all of the components diff --git a/src/fluid/drag.hpp b/src/fluid/drag.hpp index 29eca309..63b5e1bd 100644 --- a/src/fluid/drag.hpp +++ b/src/fluid/drag.hpp @@ -37,8 +37,8 @@ class GammaDrag { // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs // Get the sound speed #if HAVE_ENERGY == 1 - cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) - *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i))); + cs = std::sqrt( eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i)) + *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i)); #else cs = eos.GetWaveSpeed(k,j,i); #endif diff --git a/src/global.cpp b/src/global.cpp index 4d9499d9..1ed33ba1 100644 --- a/src/global.cpp +++ b/src/global.cpp @@ -11,6 +11,7 @@ #include "idefix.hpp" #include "global.hpp" #include "profiler.hpp" +#include "units.hpp" #ifdef WITH_MPI #include "mpi.hpp" @@ -29,6 +30,7 @@ IdefixOutStream cout; IdefixErrStream cerr; Profiler prof; LoopPattern defaultLoopPattern; +Units units; #ifdef DEBUG static int regionIndent = 0; @@ -64,6 +66,7 @@ int initialize() { void pushRegion(const std::string& kName) { Kokkos::Profiling::pushRegion(kName); if(prof.perfEnabled) { + Kokkos::fence(); prof.currentRegion = prof.currentRegion->GetChild(kName); prof.currentRegion->Start(); } diff --git a/src/global.hpp b/src/global.hpp index 265699b2..5a9ea70f 100644 --- a/src/global.hpp +++ b/src/global.hpp @@ -20,6 +20,7 @@ void safeExit(int ); // Exit the code class IdefixOutStream; class IdefixErrStream; class Profiler; +class Units; extern int prank; //< parallel rank extern int psize; @@ -29,6 +30,7 @@ extern Profiler prof; //< profiler (for memory & performance u extern double mpiCallsTimer; //< time significant MPI calls extern LoopPattern defaultLoopPattern; //< default loop patterns (for idefix_for loops) extern bool warningsAreErrors; //< whether warnings should be considered as errors +extern Units units; //< Units for the run void pushRegion(const std::string&); void popRegion(); diff --git a/src/gravity/gravity.cpp b/src/gravity/gravity.cpp index c0096c51..a39b0006 100644 --- a/src/gravity/gravity.cpp +++ b/src/gravity/gravity.cpp @@ -12,13 +12,27 @@ #include "dataBlock.hpp" #include "output.hpp" #include "input.hpp" +#include "units.hpp" Gravity::Gravity(Input &input, DataBlock *datain) { idfx::pushRegion("Gravity::Gravity"); this->data = datain; // Gravitational constant G - this->gravCst = input.GetOrSet("Gravity","gravCst",0, 1.0); + // should we compute it from the units? + if(input.CheckEntry("Gravity","gravCst")>=0) { + if(input.Get("Gravity","gravCst",0).compare("units") == 0) { + // User ask us to compute the gravitational constants from the units + + this->gravCst = idfx::units.GetDensity() * idfx::units.GetTime() * idfx::units.GetTime() + * idfx::units.G; + } else { + this->gravCst = input.Get("Gravity","gravCst",0); + } + } else { // default value to 1.0 + this->gravCst = 1.0; + } + // Gravitational potential int nPotential = input.CheckEntry("Gravity","potential"); if(nPotential >=0) { diff --git a/src/main.cpp b/src/main.cpp index d21793d1..164c134e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -30,6 +30,7 @@ #include "idefix.hpp" #include "profiler.hpp" #include "input.hpp" +#include "units.hpp" #include "grid.hpp" #include "gridHost.hpp" #include "fluid.hpp" @@ -79,6 +80,9 @@ int main( int argc, char* argv[] ) { input.PrintLogo(); idfx::cout << "Main: initialization stage." << std::endl; + // Init the units when needed + idfx::units.Init(input); + // Allocate the grid on device Grid grid(input); // Allocate the grid image on host @@ -112,6 +116,7 @@ int main( int argc, char* argv[] ) { << std::endl; } input.ShowConfig(); + idfx::units.ShowConfig(); grid.ShowConfig(); data.ShowConfig(); Tint.ShowConfig(); diff --git a/src/units.cpp b/src/units.cpp new file mode 100644 index 00000000..2baf97b1 --- /dev/null +++ b/src/units.cpp @@ -0,0 +1,55 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// Copyright(C) Geoffroy R. J. Lesur +// and other code contributors +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + +#include "units.hpp" +#include "idefix.hpp" +#include "input.hpp" + + +void idfx::Units::Init(Input &input) { + if(input.CheckBlock("Units")) { + this->_length = input.GetOrSet("Units","length",0,1.0); + this->_velocity = input.GetOrSet("Units","velocity",0,1.0); + this->_density = input.GetOrSet("Units","density",0,1.0); + this->ComputeUnits(); + } +} + +void idfx::Units::SetDensity(const real in) { + this->_density = in; + this->ComputeUnits(); +} + +void idfx::Units::SetLength(const real in) { + this->_length = in; + this->ComputeUnits(); +} + +void idfx::Units::SetVelocity(const real in) { + this->_velocity = in; + this->ComputeUnits(); +} + +void idfx::Units::ShowConfig() { + if(_isInitialized) { + idfx::cout << "Units: [Length] = " << this->_length << " cm" << std::endl; + idfx::cout << "Units: [Velocity] = " << this->_velocity << " cm/s" << std::endl; + idfx::cout << "Units: [Density] = " << this->_density << " g/cm3" << std::endl; + idfx::cout << "Units: [Energy] = " << this->_energy << " erg/cm3" << std::endl; + idfx::cout << "Units: [Time] = " << this->_time << " s" << std::endl; + idfx::cout << "Units: [Temperature] = " << this->_Kelvin << " K" << std::endl; + idfx::cout << "Units: [Mag. Field] = " << this->_magField << " G" << std::endl; + } +} + +void idfx::Units::ComputeUnits() { + this->_isInitialized = true; + this->_magField = std::sqrt(4*M_PI*_density*_velocity*_velocity); + this->_Kelvin = u * _velocity*_velocity/k_B; + this->_energy = _density*_velocity*_velocity; + this->_time = _length/_velocity; +} diff --git a/src/units.hpp b/src/units.hpp new file mode 100644 index 00000000..525e44bc --- /dev/null +++ b/src/units.hpp @@ -0,0 +1,87 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// Copyright(C) Geoffroy R. J. Lesur +// and other code contributors +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + +#ifndef UNITS_HPP_ +#define UNITS_HPP_ + +#include "idefix.hpp" + +// Units class implemented in CGS +// Constants taken from Astropy +class Input; +namespace idfx { +class Units { + public: + void Init(Input &input); + + + const real u{1.6605390666e-24}; // Atomic mass unit (g) + const real m_p{1.67262192369e-24}; // Proton mass unit (g) + const real m_n{1.67492749804e-24}; // neutron mass unit (g) + const real m_e{9.1093837015e-28}; // electron mass unit (g) + const real k_B{1.380649e-16}; // Boltzmann constant (erg/K) + const real sigma_sb{5.6703744191844314e-05}; // Stephan Boltzmann constant (g/(K^4 s^3)) + const real ar{7.5646e-15}; // Radiation constant + // = 4*sigma_sb/c (g/(K^4 s^2 cm)) + const real c{29979245800.0}; // Speed of light (cm/s) + const real M_sun{1.988409870698051e+33}; // Solar mass (g) + const real R_sun{69570000000.0}; // Solar radius (cm) + const real M_earth{5.972167867791379e+27}; // Earth mass (g) + const real R_earth{6.3781e8}; // Earth radius (cm) + const real G{6.674299999999999e-8}; // Gravitational constant (cm3 / (g s2)) + const real h{6.62607015e-27}; // Planck constant (erg.s) + const real pc{3.08568e+18}; // Parsec (cm) + const real au{1.49597892e13}; // Astronomical unit (cm) + const real e{4.80320425e-10}; // Elementary charge (cm^3/2 g^1/2 s^-1) + + + // User-defined units, non user-modifiable + // L (cm) = L (code) * Units::length + KOKKOS_INLINE_FUNCTION real GetLength() const {return _length;} + // V(cm/s) = V(code) * Units::velocity + KOKKOS_INLINE_FUNCTION real GetVelocity() const {return _velocity;} + // density (g/cm^3) = density(code) * Units::density + KOKKOS_INLINE_FUNCTION real GetDensity() const {return _density;} + + // Deduced units from user-defined units + // T(K) = P(code)/rho(code) * mu * Units::Kelvin + KOKKOS_INLINE_FUNCTION real GetKelvin() const {return _Kelvin;} + // B(G) = B(code) * Units::MagField + KOKKOS_INLINE_FUNCTION real GetMagField() const {return _magField;} + // energy(erg/cm3) = energy(code) * Units::energy + KOKKOS_INLINE_FUNCTION real GetEnergy() const {return _energy;} + // time(s) = time(code) * Units::time + KOKKOS_INLINE_FUNCTION real GetTime() const {return _time;} + + KOKKOS_INLINE_FUNCTION bool GetIsInitialized() const {return _isInitialized;} + + // code-style change of the units + void SetLength(const real); + void SetVelocity(const real); + void SetDensity(const real); + + // Show the configuration + void ShowConfig(); + + private: + bool _isInitialized{false}; + // read-write variables + real _length{1.0}; + real _velocity{1.0}; + real _density{1.0}; + + // Deduced units from user-defined units + real _Kelvin{1.0}; + real _magField{1.0}; + real _time{1.0}; + real _energy{1.0}; + + // Recompute deduced user-units + void ComputeUnits(); +}; +} // namespace idfx +#endif // UNITS_HPP_