diff --git a/.github/workflows/idefix-ci-jobs.yml b/.github/workflows/idefix-ci-jobs.yml index c0e622966..275ea1ce2 100644 --- a/.github/workflows/idefix-ci-jobs.yml +++ b/.github/workflows/idefix-ci-jobs.yml @@ -65,6 +65,8 @@ jobs: run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/OrszagTang -all $TESTME_OPTIONS - name: Orszag Tang 3D+restart tests run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/OrszagTang3D -all $TESTME_OPTIONS + - name: Axis Flux tube + run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/AxisFluxTube -all $TESTME_OPTIONS ParabolicMHD: runs-on: self-hosted @@ -177,6 +179,8 @@ jobs: run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/sphBragTDiffusion -all $TESTME_OPTIONS - name: Spherical anisotropic viscosity run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/sphBragViscosity -all $TESTME_OPTIONS + - name: Collisionless thermal conduction + run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/clessTDiffusion -all $TESTME_OPTIONS Examples: needs: [Fargo, Dust, Planet, ShearingBox, SelfGravity] @@ -201,3 +205,5 @@ jobs: run: scripts/ci/run-tests $IDEFIX_DIR/test/utils/lookupTable -all $TESTME_OPTIONS - name: Dump Image run: scripts/ci/run-tests $IDEFIX_DIR/test/utils/dumpImage -all $TESTME_OPTIONS + - name: Column density + run: scripts/ci/run-tests $IDEFIX_DIR/test/utils/columnDensity -all $TESTME_OPTIONS diff --git a/.gitignore b/.gitignore index 422647c3f..4de7fb76e 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,7 @@ doc/source/_static/* doc/source/_public/* doc/source/api/* doc/source/xml/* +doc/env/* # compiled files **/__pycache__ diff --git a/CHANGELOG.md b/CHANGELOG.md index 66f01f7be..b0732bb58 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,24 @@ 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.01] 2025-04-16 +### Changed + +- Fix a bug that led to instabilities in the RKL scheme with very small grid spacings (#323) +- Fix a bug that prevented Idefix from running with Sycl backend (required on Intel GPUs) (#331) +- Fix an error that led to incorrect electrical current regularisation around the polar axis in non-Ideal MHD (#333) +- Improve div(B) checks with a dimensionless implementation, avoiding too large divB errors in grids with large stretch factors (#334) + +### Added + +- Time-Implicit drag for multiple dust species, preventing small time steps for strongly coupled dust grains (#321) +- Collisionless heat flux added to the Braginskii module (#317) +- New global `idfx::DumpArray` debugging function to dump any Idefix array into a numpy array that can read from python (#318) +- Automatic benchmark plots in the documentation (#319) +- More CI tests of grid coarsening (#329) +- Dump outputs based on wallclock time (#335) + + ## [2.2.00] 2025-01-17 ### Changed diff --git a/CMakeLists.txt b/CMakeLists.txt index a8c7ff94e..6b62b7542 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 00) +set(Idefix_VERSION_PATCH 01) -project (idefix VERSION 2.2.00) +project (idefix VERSION 2.2.01) 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/python_requirements.txt b/doc/python_requirements.txt index 35f72dcf7..c75fa02b7 100644 --- a/doc/python_requirements.txt +++ b/doc/python_requirements.txt @@ -14,3 +14,4 @@ exhale==0.3.7 m2r2==0.3.2 sphinx-copybutton==0.5.2 #sphinxcontrib-applehelp==1.0.7 +matplotlib==3.10.0 diff --git a/doc/source/bench.json b/doc/source/bench.json new file mode 100644 index 000000000..08a3c8040 --- /dev/null +++ b/doc/source/bench.json @@ -0,0 +1,162 @@ +[ + { + "date": "2025-03-04_12:57:12", + "gpumodel": "v100", + "idefix_commit": "2bc09a0d218459f278e2b28506a09e4591b103ae", + "bench_commit": "37161676db15115c38fed3f35c94fa447cbac7bd", + "results": [ + { + "nbgpu": 1, + "cell_updates": 1.193720E+8 + }, + { + "nbgpu": 2, + "cell_updates": 1.178864E+8 + }, + { + "nbgpu": 4, + "cell_updates": 1.155336E+8 + }, + { + "nbgpu": 8, + "cell_updates": 1.014338E+8 + }, + { + "nbgpu": 16, + "cell_updates": 9.855007E+7 + }, + { + "nbgpu": 32, + "cell_updates": 9.012061E+7 + }, + { + "nbgpu": 64, + "cell_updates": 8.538461E+7 + }, + { + "nbgpu": 128, + "cell_updates": 8.531021E+7 + } + ] + }, + { + "date": "2025-03-04_13:07:10", + "gpumodel": "a100", + "idefix_commit": "2bc09a0d218459f278e2b28506a09e4591b103ae", + "bench_commit": "b536949200e50fac68d8a46d5db38fc8e3f02da5", + "results": [ + { + "nbgpu": 1, + "cell_updates": 2.044728E+8 + }, + { + "nbgpu": 2, + "cell_updates": 2.003563E+8 + }, + { + "nbgpu": 4, + "cell_updates": 1.963512E+8 + }, + { + "nbgpu": 8, + "cell_updates": 1.933039E+8 + }, + { + "nbgpu": 16, + "cell_updates": 9.759154E+7 + }, + { + "nbgpu": 32, + "cell_updates": 6.369645E+7 + }, + { + "nbgpu": 64, + "cell_updates": 4.629474E+7 + }, + { + "nbgpu": 128, + "cell_updates": 4.580281E+7 + } + ] + }, + { + "date": "2025-03-04_13:16:01", + "gpumodel": "h100", + "idefix_commit": "2bc09a0d218459f278e2b28506a09e4591b103ae", + "bench_commit": "b536949200e50fac68d8a46d5db38fc8e3f02da5", + "results": [ + { + "nbgpu": 1, + "cell_updates": 3.079643E+8 + }, + { + "nbgpu": 2, + "cell_updates": 3.012300E+8 + }, + { + "nbgpu": 4, + "cell_updates": 2.944091E+8 + }, + { + "nbgpu": 8, + "cell_updates": 2.837224E+8 + }, + { + "nbgpu": 16, + "cell_updates": 2.827778E+8 + }, + { + "nbgpu": 32, + "cell_updates": 2.822657E+8 + }, + { + "nbgpu": 64, + "cell_updates": 2.767820E+8 + }, + { + "nbgpu": 128, + "cell_updates": 2.767322E+8 + } + ] + }, + { + "date": "2025-03-06_11:21:56", + "gpumodel": "mi250x", + "idefix_commit": "2bc09a0d218459f278e2b28506a09e4591b103ae", + "bench_commit": "868be0a87c6fcda665cbb62db7020aeff70dc62d", + "results": [ + { + "nbgpu": 1, + "cell_updates": 1.436580E+8 + }, + { + "nbgpu": 2, + "cell_updates": 1.372499E+8 + }, + { + "nbgpu": 4, + "cell_updates": 1.344528E+8 + }, + { + "nbgpu": 8, + "cell_updates": 1.293602E+8 + }, + { + "nbgpu": 16, + "cell_updates": 1.260359E+8 + }, + { + "nbgpu": 32, + "cell_updates": 1.204980E+8 + }, + { + "nbgpu": 64, + "cell_updates": 1.163099E+8 + }, + { + "nbgpu": 128, + "cell_updates": 1.192343E+8 + } + ] + } +] diff --git a/doc/source/conf.py b/doc/source/conf.py index 24f9709bb..6f0cffa44 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -35,6 +35,7 @@ extensions = [ "sphinx_rtd_theme", 'sphinx_git', + 'matplotlib.sphinxext.plot_directive', "breathe", "exhale", "m2r2", diff --git a/doc/source/faq.rst b/doc/source/faq.rst index 75a07a68b..4c7bc056b 100644 --- a/doc/source/faq.rst +++ b/doc/source/faq.rst @@ -73,3 +73,6 @@ I want to test a modification to *Idefix* source code specific to my problem wit I want to use a lookup table from a CSV file in my idefix_loop. How could I proceed? Use the ``LookupTable`` class (see :ref:`LookupTableClass`) + +I want to compute a cumulative sum (e.g a column density) on the fly. How could I proceed? + Use the ``Column`` class (see :class:`::Column`) diff --git a/doc/source/index.rst b/doc/source/index.rst index 128d591ae..411d73d6d 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -84,8 +84,11 @@ Jonah Mauxion Clément Robert gitlab integration, linter -Jean Kempf & François Rincon - anisotropic diffusion +Jean Kempf, Victor Reville, François Rincon + anisotropic diffusion and collisionless thermal conduction + +Marc Coiffier + Continuous integration, automatic benchmarking ======================== About this documentation diff --git a/doc/source/modules/braginskii.rst b/doc/source/modules/braginskii.rst index e657cedd0..fd5f2b8fc 100644 --- a/doc/source/modules/braginskii.rst +++ b/doc/source/modules/braginskii.rst @@ -4,7 +4,7 @@ Braginskii module =================== Equations solved and methods ---------------------------- +---------------------------- The ``Braginskii`` module implements the anisotropic heat and momentum fluxes specific to weakly collisional, magnetised plasma like the intracluster medium @@ -37,7 +37,7 @@ though adapted to vector quantities. cell interface by a simple arithmetic average (Eq. (6)-(7) from Sharma & Hammett 2007). However in the same paper, the authors showed that this implementation can lead to - unphysical heat flux from high to low temperature regions. + unphysical heat flux from low to high temperature regions. So we implemented slope limiters for the computation of these transverse heat fluxes, as described in Eq. (17) from Sharma & Hammett (2007). Only the van Leer and the Monotonized Central (MC) limiters are available @@ -72,18 +72,68 @@ of the Braginskii heat flux and viscosity. .. _braginskiiParameterSection: + +Saturation with collisionless heat flux +--------------------------------------- + +The ``Braginskii`` module can include a collisionless saturation of the Braginskii heat flux, typically due to supra-thermal electrons. +The heat flux is then computed as follows: + +:math:`q = \alpha (q_B + q_\perp) + (1-\alpha)\beta*p*v`, + +where :math:`\alpha \in [0,1]` controls the transition between the Braginskii heat flux and the collisionless heat flux +and :math:`\beta` controls the amplitude of the collisionless heat flux (typically :math:`\beta \in [1,4]`, see Hollweg 1976). + +.. note:: + As a result, even with :math:`\kappa_\perp = 0`, the heat flux is no longer necessarilly strictly aligned with the magnetic field. +.. note:: + The collisionless heat flux is a hyperbolic term and the diffusion coefficient is set proportional to :math:`\alpha`. +.. note:: + If selected, slope limiters are also used in the collisionless flux, where an upwind scheme has been implemented for stability. +.. note:: + This saturation has been thought to be used mostly using the userdef function that takes four userdef arrays as input. + Main parameters of the module ----------------------------- The ``Braginskii`` module can be enabled adding one or two lines in the ``[Hydro]`` section -starting with the keyword -`bragTDiffusion` or/and *bragViscosity*. The following table summarises the different options +starting with the keyword `bragTDiffusion` or/and *bragViscosity*. The following tables summarise the different options associated to the activation of the Braginskii module: +--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ | Column | Entry name | Parameter type | Comment | +========+=======================+=========================+=======================================================================================+ -| 0 | bragModule | string | | Activates Braginskii diffusion. Can be ``bragTDiffusion`` or ``bragViscosity``. | +| 0 | bragTDiffusion | string | | Activates Braginskii thermal diffusion. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 1 | integration | string | | Specifies the type of scheme to be used to integrate the parabolic term. | +| | | | | Can be ``rkl`` or ``explicit``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 2 | slope limiter | string | | Choose the type of limiter to be used to compute anisotropic transverse flux terms. | +| | | | | Can be ``mc``, ``vanleer`` or ``nolimiter``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 3 | saturation mode | string | | Include or not collisionless saturation. Can be ``nosat`` or ``wcless``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 4 | diffusivity type | string | | Specifies the type of diffusivity wanted. Can be ``constant`` or ``userdef``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 5 | parallel diffusivity | real | | Mandatory if the diffusivity type is ``constant``. Not needed otherwise. | +| | | | | Value of the parallel diffusivity. Should be a real number. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 6 | normal diffusivity | real | | When bragModule ``bragTDiffusion`` and diffusivity type ``constant``, | +| | | | | value of the normal diffusivity. Should be a real number. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 7 | alpha collisionless | real | | If the diffusivity type is ``constant`` and saturation is ``wcless``. | +| | | | | Set to 1 if not provided. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 8 | beta collisionless | real | | If the diffusivity type is ``constant`` and saturation is ``wcless``. | +| | | | | Set to 0 if not provided. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ + +for the *bragViscosity*: + ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| Column | Entry name | Parameter type | Comment | ++========+=======================+=========================+=======================================================================================+ +| 0 | bragViscosity | string | | Activates Braginskii viscosity. | +--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ | 1 | integration | string | | Specifies the type of scheme to be used to integrate the parabolic term. | | | | | | Can be ``rkl`` or ``explicit``. | @@ -101,7 +151,7 @@ associated to the activation of the Braginskii module: +--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ Numerical checks ---------------- +---------------- In Cartesian geometry, the ``Braginskii`` module has been tested with many setups and in all configurations of magnetic polarisation: @@ -119,3 +169,5 @@ The same goes for the anisotropic heat flux in Cylindrical/Polar geometry while the anisotropic viscosity has *never* been tested in this geometry. In spherical geometry, both ``Braginskii`` operators have been partially validated (diffusion along the polar axis has not been directly tested). + +The collisionless saturation has been tested in 1D and 2D spherical geometry. diff --git a/doc/source/modules/dust.rst b/doc/source/modules/dust.rst index ca371dade..0cdac49ba 100644 --- a/doc/source/modules/dust.rst +++ b/doc/source/modules/dust.rst @@ -39,15 +39,45 @@ which guarantees total momentum conservation. -Drag CFL condition -------------------- -*Idefix* computes the drag terms with a time-explicit scheme. Hence, an addition CFL constraint arrises because of the drag: +Time Integration +---------------- + +*Idefix* can compute the drag with a 2nd order time-explicit (default) or 1st order time-implicit scheme. Each scheme has its pros and cons. The choice +is made by the user in the input file. + +Explicit scheme: drag CFL condition ++++++++++++++++++++++++++++++++++++ + +When using the time explicit scheme, an addition CFL constraint arrises because of the drag: .. math:: dt < \min(\frac{1}{\sum_i\gamma_i(\rho_i+\rho)}) -*Idefix* automatically adjusts the CFL to satisfy this inequality, in addition to the usual CFL condition. +*Idefix* automatically adjusts the CFL to satisfy this inequality, in addition to the usual CFL condition. This can severly limit the time step, +especially when the gas density is high. + +Implicit scheme ++++++++++++++++ + +When the drag is applied with a 1st order implicit scheme, the drag force is applied at the end of each step. In order to avoid +the complete inversion of the system, we follow a simplified inversion procedure, where we first update the gas momentum as: + +.. math:: + + v_g^{(n+1)}=\left(v_g^{(n)}+\sum_i\frac{\rho_i\gamma_i dt}{1+\rho \gamma_i dt}v_i^{(n)}\right)/\left(1+\sum_i\frac{\rho_i\gamma_i dt}{1+\rho\gamma_idt}\right) + +And then update the dust momentum as: + +.. math:: + + v_i^{(n+1)}=\left(v_i^{(n)}+\rho\gamma_i v_g^{(n+1)}dt\right)/(1+\rho\gamma_i dt) + +Note that the latter equation relies on the *updated* gas velocity. + +.. warning:: + While the implicit scheme is more stable than the explicit one, and it does not require any additional CFL condition, it is less accurate and + possibly lead to inacurrate dust velocities when :math:`dt\gg (\gamma_i\rho)^{-1}`. Use it at your own risk. Dust parameters --------------- @@ -66,6 +96,9 @@ The dust module can be enabled adding a block `[Dust]` in your input .ini file. +----------------+-------------------------+---------------------------------------------------------------------------------------------+ | drag_feedback | bool | | (optionnal) whether the gas feedback is enabled (default true). | +----------------+-------------------------+---------------------------------------------------------------------------------------------+ +| drag_implicit | bool | | (optionnal) whether the drag uses a 1st order implicit method. Otherwise use the | +| | | | 2nd order time-explicit scheme (default is false=time explicit) | ++----------------+-------------------------+---------------------------------------------------------------------------------------------+ The drag parameter :math:`\beta_i` above sets the functional form of :math:`\gamma_i(\rho, \rho_i, c_s)` depending on the drag type: diff --git a/doc/source/performances.rst b/doc/source/performances.rst index 84042685e..deae1194e 100644 --- a/doc/source/performances.rst +++ b/doc/source/performances.rst @@ -6,9 +6,8 @@ We report below the performances obtained on various architectures using Idefix. is the 3D MHD Orszag-Tang test problem with 2nd order reconstruction and uct_contact EMFS bundled in Idefix test suite, disabling passive tracers. The test is computed with a 128\ :sup:`3` resolution per MPI sub-domain on GPUs or 32\ :sup:`3` per MPI sub-domain on CPUs. All of the performances measures -have been obtained enabling MPI on *one full node*, but we report here the performance *per GPU* -(i.e. with 2 GCDs on AMD Mi250) or *per core* (on CPU), i.e. dividing the node performance by the number of GPU/core -to simplify the comparison with other clusters. +have been obtained enabling MPI and we reporte here the performance *per GPU*, *per GCD* (on Mi250) +or *per core* (on CPU). The complete scalability tests are available in Idefix `method paper `_. The performances mentionned below are updated for each major revision of Idefix, so they might slightly differ from the method paper. @@ -33,16 +32,14 @@ CPU performances | IDRIS/Jean Zay | Intel Cascade Lake | 0.62 | +---------------------+--------------------+----------------------------------------------------+ - GPU performances ================ -+----------------------+--------------------+----------------------------------------------------+ -| Cluster name | GPU | Performances (in 10\ :sup:`6` cell/s/GPU) | -+======================+====================+====================================================+ -| IDRIS/Jean Zay | NVIDIA V100 | 110 | -+----------------------+--------------------+----------------------------------------------------+ -| IDRIS/Jean Zay | NVIDIA A100 | 194 | -+----------------------+--------------------+----------------------------------------------------+ -| CINES/Adastra | AMD Mi250 | 250 | -+----------------------+--------------------+----------------------------------------------------+ +.. plot:: + + import plot_idefix_bench + plot_idefix_bench.do_plot('Performance on NVidia and AMD GPUs', 'bench.json', ['v100','a100','h100','mi250x']) + +.. note:: + + The inter-node communication on Jean Zay is not optimal on A100 nodes. A ticket is opened with IDRIS support to fix this issue. diff --git a/doc/source/plot_idefix_bench.py b/doc/source/plot_idefix_bench.py new file mode 100644 index 000000000..43482eb04 --- /dev/null +++ b/doc/source/plot_idefix_bench.py @@ -0,0 +1,26 @@ +import matplotlib.pyplot as plt +import json + +def do_plot(title, bench_file, gpumodels): + with open(bench_file, 'r') as f: + benches = json.load(f) + + plt.figure() + xmax=0 + ymax=0 + for gpumodel in gpumodels: + select = [bench for bench in benches if bench['gpumodel'] == gpumodel][-1] + + xs = [r['nbgpu'] for r in select['results']] + ys = [r['cell_updates'] for r in select['results']] + plt.plot(xs, ys,'o-',label=gpumodel) + xmax=max(xmax,max(xs)) + ymax=max(ymax,max(ys)) + + plt.xscale("log", base=2) + plt.ylim(0,ymax*1.1) + plt.xlim(1,xmax*1.1) + plt.legend() + plt.xlabel("Number of GPUs/GCDs") + plt.ylabel("Performance (cells / second / GPU)") + plt.title(title) diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index 18795b50a..4bd7c69f7 100644 --- a/doc/source/programmingguide.rst +++ b/doc/source/programmingguide.rst @@ -645,6 +645,27 @@ It may also be useful to implement debug-only safeguards with custom logic that fit `RUNTIME_CHECK_*` macros. This can be achieved by using the compiler directive `#ifdef RUNTIME_CHECKS` directly. +Dump an array to a file +----------------------- + +It is usually difficult to know what Idefix arrays effectively contains, especially when running on GPU. +To help with this difficulty, Idefix provides a global function ``DumpArray`` that can be used +to dump a ``IdefixArray`` to a numpy file (that can be read from python). This feature can be used +for debugging purpose as: + + +.. code-block:: c++ + + #include "idefix.hpp" + + IdefixArray3D myArray("debugMe",10,10,10); + + idfx::DumpArray("myFilename.npy",myArray); // Dump the array content to a numpy file named "myFilename.npy" + + +Note that the array is automatically transfered from the GPU, if needed. + + Minimal skeleton ================ diff --git a/doc/source/reference/idefix.ini.rst b/doc/source/reference/idefix.ini.rst index c9f9c5785..1fa05cc67 100644 --- a/doc/source/reference/idefix.ini.rst +++ b/doc/source/reference/idefix.ini.rst @@ -388,8 +388,10 @@ This section describes the outputs *Idefix* produces. For more details about eac +================+=========================+==================================================================================================+ | log | integer | | Time interval between log outputs, in code steps (default 100). | +----------------+-------------------------+--------------------------------------------------------------------------------------------------+ -| dmp | float | | Time interval between dump outputs, in code units. | -| | | | If negative, periodic dump outputs are disabled. | +| dmp | float, float+char | | 1st parameter: Code time interval between dump outputs, in code units. | +| | | | If negative, the first parameter is ignored. | +| | | | 2nd parameter (optional): Wallclock time interval between two dumps. The ending character | +| | | | can be "s" (seconds) "m" (minutes) "h" (hours) or "d" (days) | +----------------+-------------------------+--------------------------------------------------------------------------------------------------+ | dmp_dir | string | | directory for dump file outputs. Default to "./" | | | | | The directory is automatically created if it does not exist. | diff --git a/doc/source/reference/makefile.rst b/doc/source/reference/makefile.rst index aadc205c3..23052d92c 100644 --- a/doc/source/reference/makefile.rst +++ b/doc/source/reference/makefile.rst @@ -111,7 +111,8 @@ We recommend the following modules and environement variables on AdAstra: module load cpe/24.07 module load craype-accel-amd-gfx90a craype-x86-trento module load PrgEnv-cray - module load amd-mixed + module load amd-mixed/6.1.2 + module load rocm/6.1.2 module load cray-python/3.11.7 Finally, *Idefix* can be configured to run on Mi250 by enabling HIP and the desired architecture with the following options to ccmake: diff --git a/reference b/reference index b675bceaa..c32894208 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit b675bceaa6aabc01dded346e2d631857f698dc76 +Subproject commit c328942083b618a85b84eb16ce9fd35e67c00597 diff --git a/src/dataBlock/coarsen.cpp b/src/dataBlock/coarsen.cpp index 2da3049cb..4c6b60ccd 100644 --- a/src/dataBlock/coarsen.cpp +++ b/src/dataBlock/coarsen.cpp @@ -66,27 +66,33 @@ void DataBlock::CheckCoarseningLevels() { idfx::pushRegion("DataBlock::CheckCoarseningLevels()"); // Check that the coarsening levels we have are valid // NB: this is a costly procedure, we can't repeat it at each loop! - DataBlockHost d(*this); - d.SyncFromDevice(); for(int dir = 0 ; dir < DIMENSIONS ; dir++) { if(mygrid->coarseningDirection[dir]) { - IdefixHostArray2D arr = d.coarseningLevel[dir]; - const int Xt = (dir == IDIR ? JDIR : IDIR); - const int Xb = (dir == KDIR ? JDIR : KDIR); - for(int i = beg[Xt] ; i < end[Xt] ; i++) { - for(int j = beg[Xb] ; j < end[Xb] ; j++) { + IdefixHostArray2D arr = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), + coarseningLevel[dir]); + for(int j = 0 ; j < arr.extent(0) ; j++) { + for(int i = 0 ; i < arr.extent(1) ; i++) { + if(std::isnan(arr(j,i))) { + std::stringstream str; + str << "Nan in grid coarsening levels" << std::endl; + str << "at (i,j)=("<< i << "," << j << "): Coarsening level is NaN!" << std::endl; + IDEFIX_ERROR(str); + } if(arr(j,i) < 1) { std::stringstream str; - str << "Incorrect grid coarsening levels" << std::endl; - str << "at (i,j)=("<< i << "," << j << "): Coarsening level < 1!" << std::endl; + str << "Coarsening level < 1!" << std::endl; + str << "at (i,j)=("<< i << "," << j << "): "; + str << "coarsening level= " << arr(j,i) << std::endl; IDEFIX_ERROR(str); } const int factor = 1 << (arr(j,i) - 1); if(np_int[dir] % factor != 0) { std::stringstream str; - str << "local grid size not divisible by coarsening level" << std::endl; - str << "at (i,j)=("<< i << "," << j << "): Coarsening level: "; - str << arr(j,i) << std::endl; + str << "Local grid size not divisible by coarsening level." << std::endl; + str << "at (i,j)=("<< i << "," << j << "): "; + str << "coarsening level= " << arr(j,i) << std::endl; + str << np_int[dir] << " cannot be divided by 2^" << arr(j,i)-1; + str << " = " << factor << std::endl; IDEFIX_ERROR(str); } } diff --git a/src/dataBlock/evolveStage.cpp b/src/dataBlock/evolveStage.cpp index 981c060ae..c97c9f713 100644 --- a/src/dataBlock/evolveStage.cpp +++ b/src/dataBlock/evolveStage.cpp @@ -19,6 +19,16 @@ void DataBlock::EvolveStage() { for(int i = 0 ; i < dust.size() ; i++) { dust[i]->EvolveStage(this->t,this->dt); } + // Add implicit term for dust drag + if(dust[0]->drag->IsImplicit()) { + for(int i = 0 ; i < dust.size() ; i++) { + dust[i]->drag->AddImplicitBackReaction(this->dt,dust[0]->drag->implicitFactor); + } + dust[0]->drag->NormalizeImplicitBackReaction(this->dt); + for(int i = 0 ; i < dust.size() ; i++) { + dust[i]->drag->AddImplicitFluidMomentum(this->dt); + } + } } idfx::popRegion(); diff --git a/src/fluid/boundary/axis.cpp b/src/fluid/boundary/axis.cpp index b77fe3b5f..3ad67f0a0 100644 --- a/src/fluid/boundary/axis.cpp +++ b/src/fluid/boundary/axis.cpp @@ -128,12 +128,12 @@ void Axis::RegularizeCurrentSide(int side) { real deltaPhi = data->mygrid->xend[KDIR] - data->mygrid->xbeg[KDIR]; // Use the circulation around the pole of Bphi to determine Jr on the pole: - // Delta phi r^2(1-cos theta) Jr = int r sin(theta) Bphi dphi + // 1/2*Delta phi r^2 sin theta^2 Jr = int r sin(theta) Bphi dphi idefix_for("fixJ",0,data->np_tot[KDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA(int k,int i) { real th = x2(jc); - real fact = sign*sin(th)/(deltaPhi*x1(i)*(1-cos(th))); + real fact = 2*sign/(deltaPhi*x1(i)*sin(th)); J(IDIR, k,js,i) = BAvg(i)*fact; }); diff --git a/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index 82d5dac94..f251d0b0c 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.cpp +++ b/src/fluid/braginskii/bragThermalDiffusion.cpp @@ -19,7 +19,6 @@ #include "eos.hpp" - void BragThermalDiffusion::ShowConfig() { if(status.status==Constant) { idfx::cout << "Braginskii Thermal Diffusion: ENABLED with constant diffusivity kpar=" @@ -27,11 +26,11 @@ void BragThermalDiffusion::ShowConfig() { } else if (status.status==UserDefFunction) { idfx::cout << "Braginskii Thermal Diffusion: ENABLED with user-defined diffusivity function." << std::endl; - if(!diffusivityFunc) { + if(!bragDiffusivityFunc) { IDEFIX_ERROR("No braginskii thermal diffusion function has been enrolled"); } } else { - IDEFIX_ERROR("Unknown braginskii thermal diffusion mode"); + IDEFIX_ERROR("Unknown Braginskii thermal diffusion mode"); } if(status.isExplicit) { idfx::cout << "Braginskii Thermal Diffusion: uses an explicit time integration." << std::endl; @@ -42,16 +41,27 @@ void BragThermalDiffusion::ShowConfig() { IDEFIX_ERROR("Unknown time integrator for braginskii thermal diffusion."); } if(haveSlopeLimiter) { - idfx::cout << "Braginskii Thermal Diffusion: uses a slope limiter." << std::endl; + if(haveMonotonizedCentral) { + idfx::cout << "Braginskii Thermal Diffusion: " + "uses the monotonized central slope limiter." << std::endl; + } else if(haveVanLeer) { + idfx::cout << "Braginskii Thermal Diffusion: uses the van Leer slope limiter." << std::endl; + } else { + IDEFIX_ERROR("Unknown slope limiter for braginskii thermal diffusion."); + } + } + if(includeCollisionlessTD) { + idfx::cout << "Braginskii Thermal Diffusion: saturation" + " with collisionless flux is enabled." << std::endl; } } void BragThermalDiffusion::EnrollBragThermalDiffusivity(BragDiffusivityFunc myFunc) { if(this->status.status != UserDefFunction) { IDEFIX_WARNING("Braginskii thermal diffusivity enrollment requires Hydro/BragThermalDiffusion " - "to be set to userdef in .ini file"); + "to be set to userdef in .ini file"); } - this->diffusivityFunc = myFunc; + this->bragDiffusivityFunc = myFunc; } void BragThermalDiffusion::AddBragDiffusiveFlux(int dir, const real t, diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 1fe774a41..116ff9a03 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -9,6 +9,7 @@ #define FLUID_BRAGINSKII_BRAGTHERMALDIFFUSION_HPP_ #include +#include #include "idefix.hpp" #include "input.hpp" @@ -40,6 +41,8 @@ class BragThermalDiffusion { IdefixArray3D heatSrc; // Source terms of the thermal operator IdefixArray3D knorArr; IdefixArray3D kparArr; + IdefixArray3D clessAlphaArr; // Transition from Brag to collisionless tc + IdefixArray3D clessBetaArr; // Collisionless tc flux coefficient (before pv) // pre-computed geometrical factors in non-cartesian geometry IdefixArray1D one_dmu; @@ -50,9 +53,14 @@ class BragThermalDiffusion { // status of the module ParabolicModuleStatus &status; - BragDiffusivityFunc diffusivityFunc; + BragDiffusivityFunc bragDiffusivityFunc; + + bool includeCollisionlessTD{false}; bool haveSlopeLimiter{false}; + bool haveMonotonizedCentral{false}; + bool haveVanLeer{false}; + bool haveMinmod{false}; // helper array IdefixArray4D &Vc; @@ -61,6 +69,7 @@ class BragThermalDiffusion { // constant diffusion coefficient (when needed) real knor, kpar; + real clessAlpha, clessBeta; // equation of state (required to get the heat capacity) EquationOfState *eos; @@ -84,9 +93,11 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid if(input.CheckEntry("Hydro","bragTDiffusion")>=0) { if(input.Get("Hydro","bragTDiffusion",1).compare("mc") == 0) { this->haveSlopeLimiter = true; + this->haveMonotonizedCentral = true; limiter = PLMLimiter::McLim; } else if(input.Get("Hydro","bragTDiffusion",1).compare("vanleer") == 0) { this->haveSlopeLimiter = true; + this->haveVanLeer = true; limiter = PLMLimiter::VanLeer; } else if(input.Get("Hydro","bragTDiffusion",1).compare("minmod") == 0) { IDEFIX_ERROR("The minmod slope limiter is not available because it has been " @@ -98,21 +109,55 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid IDEFIX_ERROR("Unknown braginskii thermal diffusion limiter in idefix.ini. " "Can only be vanleer, mc or nolimiter."); } - if(input.Get("Hydro","bragTDiffusion",2).compare("constant") == 0) { - this->kpar = input.Get("Hydro","bragTDiffusion",3); - this->knor = input.GetOrSet("Hydro","bragTDiffusion",4,0.); - this->status.status = Constant; - } else if(input.Get("Hydro","bragTDiffusion",2).compare("userdef") == 0) { - this->status.status = UserDefFunction; - this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); + if(input.Get("Hydro","bragTDiffusion",2).compare("nosat") == 0) { + if(input.Get("Hydro","bragTDiffusion",3).compare("constant") == 0) { + this->kpar = input.Get("Hydro","bragTDiffusion",4); + this->knor = input.GetOrSet("Hydro","bragTDiffusion",5,0.); + this->status.status = Constant; + } else if(input.Get("Hydro","bragTDiffusion",3).compare("userdef") == 0) { + this->status.status = UserDefFunction; + this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } else { + IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " + "Can only be constant or userdef."); + } + } else if(input.Get("Hydro","bragTDiffusion",2).compare("wcless") == 0.0) { + if(input.Get("Hydro","bragTDiffusion",3).compare("constant") == 0) { + this->includeCollisionlessTD = true; + this->kpar = input.Get("Hydro","bragTDiffusion",4); + this->knor = input.GetOrSet("Hydro","bragTDiffusion",5,0.); + this->clessAlpha = input.GetOrSet("Hydro","bragTDiffusion",6,1.); + this->clessBeta = input.GetOrSet("Hydro","bragTDiffusion",7,0.); + this->status.status = Constant; + } else if(input.Get("Hydro","bragTDiffusion",3).compare("userdef") == 0) { + this->includeCollisionlessTD = true; + this->status.status = UserDefFunction; + this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->clessAlphaArr = IdefixArray3D("ClessThermalDiffusionAlphaArray", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->clessBetaArr = IdefixArray3D("ClessThermalDiffusionBetaArray", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } else { + IDEFIX_ERROR("Unknown braginskii/collisionless thermal diffusion definition in idefix.ini. " + "Can only be constant or userdef."); + } } else { - IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " - "Can only be constant or userdef."); + IDEFIX_ERROR("Unknown braginskii thermal diffusion saturation in idefix.ini. " + "Can only be nosat or wcless."); } } else { IDEFIX_ERROR("I cannot create a BragThermalDiffusion object without bragTDiffusion defined" @@ -135,7 +180,7 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid // and therefore not available as such in the DataBlock. // It is rather defined as PRS/RHO. // Special spatial derivative macros are therefore needed and defined here -// directly at the right cell interface according to the direciton of the flux. +// directly at the right cell interface according to the direction of the flux. #define D_DX_I_T(q) (q(PRS,k,j,i)/q(RHO,k,j,i) - q(PRS,k,j,i - 1)/q(RHO,k,j,i - 1)) #define D_DY_J_T(q) (q(PRS,k,j,i)/q(RHO,k,j,i) - q(PRS,k,j - 1,i)/q(RHO,k,j - 1,i)) #define D_DZ_K_T(q) (q(PRS,k,j,i)/q(RHO,k,j,i) - q(PRS,k - 1,j,i)/q(RHO,k - 1,j,i)) @@ -179,7 +224,7 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid //We now define spatial average macros for the magnetic field. // The magnetic field appears in the expression of the Braginskii heat flux. -// It is therefore needed at the right cell interface according to the direction of the flux. +// It is therefore needed at the left cell interface according to the direction of the flux. #define BX_I Vs(BX1s,k,j,i) #define BY_J Vs(BX2s,k,j,i) #define BZ_K Vs(BX3s,k,j,i) @@ -210,6 +255,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, HydroModuleStatus haveThermalDiffusion = this->status.status; bool haveSlopeLimiter = this->haveSlopeLimiter; + bool includeCollisionlessTD = this->includeCollisionlessTD; using SL = SlopeLimiter; int ibeg, iend, jbeg, jend, kbeg, kend; @@ -250,16 +296,33 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, #endif real knorConstant = this->knor; real kparConstant = this->kpar; + real clessAlphaConst = this->clessAlpha; + real clessBetaConst = this->clessBeta; IdefixArray3D knorArr = this->knorArr; IdefixArray3D kparArr = this->kparArr; + IdefixArray3D clessAlphaArr = this->clessAlphaArr; + IdefixArray3D clessBetaArr = this->clessBetaArr; - if(haveThermalDiffusion == UserDefFunction && dir == IDIR) { - if(diffusivityFunc) { + if(includeCollisionlessTD == false && haveThermalDiffusion == UserDefFunction && dir == IDIR) { + if(bragDiffusivityFunc) { idfx::pushRegion("UserDef::BragThermalDiffusivityFunction"); - diffusivityFunc(*this->data, t, kparArr, knorArr); + std::vector> userdefArr = {kparArr, knorArr}; + bragDiffusivityFunc(*this->data, t, userdefArr); + idfx::popRegion(); + } + else { + IDEFIX_ERROR("No user-defined Braginskii thermal diffusion function has been enrolled"); + } + } else if(includeCollisionlessTD == true && haveThermalDiffusion == UserDefFunction + && dir == IDIR) { + if (bragDiffusivityFunc) { + idfx::pushRegion("UserDef::ClessThermalDiffusivityFunction"); + std::vector> userdefArr = {kparArr, knorArr, clessAlphaArr, clessBetaArr}; + bragDiffusivityFunc(*this->data, t, userdefArr); idfx::popRegion(); } else { - IDEFIX_ERROR("No user-defined thermal diffusion function has been enrolled"); + IDEFIX_ERROR("No user-defined Braginskii/collisionless " + "thermal diffusion function has been enrolled"); } } @@ -273,6 +336,8 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, [[maybe_unused]] real dTi, dTj, dTk, dTn; dTi = dTj = dTk = dTn = ZERO_F; + /* For the collisionless / saturated flux */ + [[maybe_unused]] real clessAlpha, clessBeta, dvp, dvm, dpp, dpm, Pn, Vn; real locdmax = 0; /////////////////////////////////////////// @@ -287,14 +352,22 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } else { kpar = HALF_F*(kparArr(k,j,i-1)+kparArr(k,j,i)); } + if(includeCollisionlessTD) { + clessAlpha=HALF_F*(clessAlphaArr(k,j,i-1)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k,j,i-1)+clessBetaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } - EXPAND( Bi = BX_I; , - Bj = BY_I; , - Bk = BZ_I; ) + D_EXPAND( Bi = BX_I; , + Bj = BY_I; , + Bk = BZ_I; ) Bn = BX_I; #if GEOMETRY == CARTESIAN @@ -370,7 +443,9 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } dTi = D_DX_I_T(Vc)/dx1(i); + #if DIMENSIONS >= 2 + if (haveSlopeLimiter) { dTj = 1./x1(i-1)*SL::PLMLim(SL_DY_T(Vc,k,j,i-1)/dx2(j), SL_DY_T(Vc,k,j+1,i-1)/dx2(j+1)); @@ -399,6 +474,30 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, locdmax = FMAX(kpar,knor)/(0.5*(Vc(RHO,k,j,i) + Vc(RHO,k,j,i - 1)))*gamma_m1; dTn = dTi; + + /* For collisionless / saturated heat flux */ + if (includeCollisionlessTD) { + if(haveSlopeLimiter) { + dvp = Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i); + dvm = Vc(VX1,k,j,i) - Vc(VX1,k,j,i-1); + + dpp = Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); + + /* Upwind scheme */ + if (Vc(VX1,k,j,i) > 0.0) { + Vn = Vc(VX1,k,j,i-1)+HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i-1)+HALF_F*SL::PLMLim(dpm, dpp); + } + else { + Vn = Vc(VX1,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); + } + } else { + Vn = HALF_F*(Vc(VX1,k,j,i) + Vc(VX1,k,j,i-1)); + Pn = HALF_F*(Vc(PRS,k,j,i) + Vc(PRS,k,j,i-1)); + } + } } else if(dir == JDIR) { ////////////// // JDIR @@ -411,9 +510,17 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } else { kpar = HALF_F*(kparArr(k,j-1,i)+kparArr(k,j,i)); } + if(includeCollisionlessTD) { + clessAlpha=HALF_F*(clessAlphaArr(k,j-1,i)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k,j-1,i)+clessBetaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } EXPAND( Bi = BX_J; , @@ -533,6 +640,29 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, locdmax = FMAX(kpar,knor)/(0.5*(Vc(RHO,k,j,i) + Vc(RHO,k,j - 1,i)))*gamma_m1; dTn = dTj; + + /* For the collisionless / saturated flux */ + if(includeCollisionlessTD) { + if(haveSlopeLimiter) { + dvp = Vc(VX2,k,j+1,i) - Vc(VX2,k,j,i); + dvm = Vc(VX2,k,j,i) - Vc(VX2,k,j-1,i); + + dpp = Vc(PRS,k,j+1,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); + + /* Upwind scheme */ + if (Vc(VX2,k,j,i) > 0.0) { + Vn = Vc(VX2,k,j-1,i)+HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j-1,i)+HALF_F*SL::PLMLim(dpm, dpp); + } else { + Vn = Vc(VX2,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); + } + } else { + Vn = HALF_F*(Vc(VX2,k,j-1,i) + Vc(VX2,k,j,i)); + Pn = HALF_F*(Vc(PRS,k,j-1,i) + Vc(PRS,k,j,i)); + } + } } else if(dir == KDIR) { ////////////// // KDIR @@ -544,9 +674,17 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } else { kpar = HALF_F*(kparArr(k-1,j,i)+kparArr(k,j,i)); } + if(includeCollisionlessTD) { + clessAlpha=HALF_F*(clessAlphaArr(k-1,j,i)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k-1,j,i)+clessBetaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } Bi = BX_K; @@ -625,11 +763,36 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, locdmax = FMAX(kpar,knor)/(0.5*(Vc(RHO,k,j,i) + Vc(RHO,k-1,j,i)))*gamma_m1; dTn = dTk; + + /* For the collisionless / saturated flux */ + if(includeCollisionlessTD) { + if(haveSlopeLimiter) { + dvp = Vc(VX3,k+1,j,i) - Vc(VX3,k,j,i); + dvm = Vc(VX3,k,j,i) - Vc(VX3,k-1,j,i); + + dpp = Vc(PRS,k+1,j,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); + + /* Upwind scheme */ + if (Vc(VX3,k,j,i) > 0.0) { + Vn = Vc(VX3,k-1,j,i)+HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k-1,j,i)+HALF_F*SL::PLMLim(dpp, dpm); + } else { + Vn = Vc(VX3,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); + } + } else { + Vn = HALF_F*(Vc(VX3,k-1,j,i) + Vc(VX3,k,j,i)); + Pn = HALF_F*(Vc(PRS,k-1,j,i) + Vc(PRS,k,j,i)); + } + } } // From here, gradients and normal have been computed, so we just need to get the fluxes - bgradT = EXPAND( Bi*dTi , + Bj*dTj, +Bk*dTk); - Bmag = EXPAND( Bi*Bi , + Bj*Bj, + Bk*Bk); + bgradT = D_EXPAND( Bi*dTi , + Bj*dTj, +Bk*dTk); + Bmag = D_EXPAND( Bi*Bi , + Bj*Bj, + Bk*Bk); + // EXPAND can yield unexpected behaviour when DIMENSIONS < COMPONENTS + //printf("%f , %f\n", Bmag, EXPAND(Bi*Bi, + Bj*Bj, + Bk*Bk)); Bmag = sqrt(Bmag); Bmag = FMAX(1e-6*SMALL_NUMBER,Bmag); @@ -637,10 +800,15 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, bn = Bn/Bmag; /* -- unit vector component -- */ q = kpar*bgradT*bn + knor*(dTn - bn*bgradT); - + if(includeCollisionlessTD) { + q = clessAlpha*q + (1-clessAlpha)*clessBeta*Pn*Vn; + } Flux(ENG, k, j, i) -= q; - dMax(k,j,i) = FMAX(dMax(k,j,i),locdmax); + + if(includeCollisionlessTD) { + dMax(k,j,i) *= clessAlpha; + } }); idfx::popRegion(); } diff --git a/src/fluid/braginskii/bragViscosity.cpp b/src/fluid/braginskii/bragViscosity.cpp index 6d5007b50..e53ab2c65 100644 --- a/src/fluid/braginskii/bragViscosity.cpp +++ b/src/fluid/braginskii/bragViscosity.cpp @@ -59,7 +59,14 @@ void BragViscosity::ShowConfig() { IDEFIX_ERROR("Unknown time integrator for braginskii viscosity."); } if(haveSlopeLimiter) { - idfx::cout << "Braginskii Viscosity: uses a slope limiter." << std::endl; + if(haveMonotizedCentral) { + idfx::cout << "Braginskii Viscosity: uses the " + "monotonized central slope limiter." << std::endl; + } else if(haveVanLeer) { + idfx::cout << "Braginskii Viscosity: uses the van Leer slope limiter." << std::endl; + } else { + IDEFIX_ERROR("Unknown slope limiter for braginskii viscosity."); + } } #if GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR diff --git a/src/fluid/braginskii/bragViscosity.hpp b/src/fluid/braginskii/bragViscosity.hpp index 7a8d2b756..5340110a8 100644 --- a/src/fluid/braginskii/bragViscosity.hpp +++ b/src/fluid/braginskii/bragViscosity.hpp @@ -52,6 +52,8 @@ class BragViscosity { DiffusivityFunc bragViscousDiffusivityFunc; bool haveSlopeLimiter{false}; + bool haveMonotizedCentral{false}; + bool haveVanLeer{false}; IdefixArray4D &Vc; IdefixArray4D &Vs; @@ -78,9 +80,11 @@ BragViscosity::BragViscosity(Input &input, Grid &grid, Fluid *hydroin): if(input.CheckEntry("Hydro","bragViscosity")>=0) { if(input.Get("Hydro","bragViscosity",1).compare("vanleer") == 0) { this->haveSlopeLimiter = true; + this->haveVanLeer = true; limiter = PLMLimiter::VanLeer; } else if(input.Get("Hydro","bragViscosity",1).compare("mc") == 0) { this->haveSlopeLimiter = true; + this->haveMonotizedCentral = true; limiter = PLMLimiter::McLim; } else if (input.Get("Hydro","bragViscosity",1).compare("nolimiter") == 0) { this->haveSlopeLimiter = false; diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp index c8738b4b2..38ef76d5d 100644 --- a/src/fluid/calcRightHandSide.hpp +++ b/src/fluid/calcRightHandSide.hpp @@ -149,47 +149,49 @@ struct Fluid_CorrectFluxFunctor { Flux(MX1+meanDir,k,j,i) += meanV * Flux(RHO,k,j,i); } // Fargo & Rotation corrections - real Ax = A(k,j,i); + ////////////////////////////////////////////// + // Define correction factor for the fluxes + ////////////////////////////////////////////// - for(int nv = 0 ; nv < Phys::nvar ; nv++) { - Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax; - } + real Ax[Phys::nvar]; // corrected Area + for(int nv = 0 ; nv < Phys::nvar ; nv++) Ax[nv] = A(k,j,i); // Curvature terms -#if (GEOMETRY == POLAR && COMPONENTS >= 2) \ - || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) - if constexpr (dir==IDIR) { - // Conserve angular momentum, hence flux is R*Bphi - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - if constexpr(Phys::mhd) { - if(Ax= 2) \ + || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) + if constexpr (dir==IDIR) { + // Conserve angular momentum, hence flux is R*Bphi + Ax[iMPHI] *= FABS(x1m(i)); + if constexpr(Phys::mhd) { + Ax[iBPHI] = 1.0; // corrected Flux is simply Bphi + } } - } -#endif // GEOMETRY==POLAR OR CYLINDRICAL + #endif // GEOMETRY==POLAR OR CYLINDRICAL -#if GEOMETRY == SPHERICAL - if constexpr(dir==IDIR) { - #if COMPONENTS == 3 - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - #endif // COMPONENTS == 3 - if constexpr(Phys::mhd) { - if(Ax::CheckDivB() { data->beg[IDIR], data->end[IDIR], KOKKOS_LAMBDA (int k, int j, int i, real &divBmax) { [[maybe_unused]] real dB1,dB2,dB3; + [[maybe_unused]] real d1, d2, d3; + [[maybe_unused]] real B1,B2,B3; dB1=dB2=dB3=ZERO_F; + d1=d2=d3=ZERO_F; + B1=B2=B3=ZERO_F; + D_EXPAND( dB1=(Ax1(k,j,i+1)*Vs(BX1s,k,j,i+1)-Ax1(k,j,i)*Vs(BX1s,k,j,i)); , dB2=(Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i)-Ax2(k,j,i)*Vs(BX2s,k,j,i)); , dB3=(Ax3(k+1,j,i)*Vs(BX3s,k+1,j,i)-Ax3(k,j,i)*Vs(BX3s,k,j,i)); ) - divBmax=FMAX(FABS(D_EXPAND(dB1, +dB2, +dB3))/dV(k,j,i),divBmax); + D_EXPAND( d1=0.5*(Ax1(k,j,i+1) + Ax1(k,j,i)); , + d2=0.5*(Ax2(k,j+1,i) + Ax2(k,j,i)); , + d3=0.5*(Ax3(k+1,j,i) + Ax3(k,j,i)); ) + + D_EXPAND( B1=0.5*(Vs(BX1s,k,j,i+1) + Vs(BX1s,k,j,i)); , + B2=0.5*(Vs(BX2s,k,j+1,i) + Vs(BX2s,k,j,i)); , + B3=0.5*(Vs(BX3s,k+1,j,i) + Vs(BX3s,k,j,i)); ) + + real amplitude = 1e-40; + amplitude += D_EXPAND( std::fabs(B1)*d1, + std::fabs(B2)*d2, + std::fabs(B3)*d3 ); + + divBmax=FMAX(FABS(D_EXPAND(dB1, +dB2, +dB3))/amplitude,divBmax); }, Kokkos::Max(divB) // reduction ); diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index d24e0bea7..4cf879cdd 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -5,7 +5,9 @@ // Licensed under CeCILL 2.1 License, see COPYING for more information // *********************************************************************************** #include "drag.hpp" +#include #include "physics.hpp" + void Drag::AddDragForce(const real dt) { idfx::pushRegion("Drag::AddDragForce"); @@ -15,83 +17,182 @@ void Drag::AddDragForce(const real dt) { auto VcDust = this->VcDust; auto InvDt = this->InvDt; - const Type type = this->type; - real dragCoeff = this->dragCoeff; bool feedback = this->feedback; - EquationOfState eos = *(this->eos); - - auto userGammai = this->gammai; - if(type == Type::Userdef) { - if(userDrag != NULL) { - idfx::pushRegion("Drag::UserDrag"); - userDrag(data, dragCoeff, userGammai); - idfx::popRegion(); - } else { - IDEFIX_ERROR("No User-defined drag function has been enrolled"); - } + if(implicit) { + IDEFIX_ERROR("Add DragForce should not be called when drag is implicit"); } + + auto gammaDrag = this->gammaDrag; + gammaDrag.RefreshUserDrag(data); + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) // Where gamma is computed according to the choice of drag type idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA (int k, int j, int i) { - real gamma; // The drag coefficient - if(type == Type::Gamma) { - gamma = dragCoeff; - - } else if(type == Type::Tau) { - // In this case, the coefficient is the stopping time (assumed constant) - gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); - } else if(type == Type::Size) { - real cs; - // 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))); - #else - cs = eos.GetWaveSpeed(k,j,i); - #endif - gamma = cs/dragCoeff; - } else if(type == Type::Userdef) { - gamma = userGammai(k,j,i); - } + real gamma = gammaDrag.GetGamma(k,j,i); // The drag coefficient real dp = dt * gamma * VcDust(RHO,k,j,i) * VcGas(RHO,k,j,i); for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { real dv = VcDust(n,k,j,i) - VcGas(n,k,j,i); UcDust(n,k,j,i) -= dp*dv; - if(feedback) UcGas(n,k,j,i) += dp*dv; + if(feedback) { + UcGas(n,k,j,i) += dp*dv; + #if HAVE_ENERGY == 1 + // We add back the energy dissipated for the dust which is not accounted for + // (since there is no energy equation for dust grains) + + // TODO(GL): this should be disabled in the case of a true multifluid system where + // both fluids have a proper energy equation + UcGas(ENG,k,j,i) += dp*dv*VcDust(n,k,j,i); + #endif + } // feedback + } + // Cfl constraint + real idt = gamma*VcGas(RHO,k,j,i); + if(feedback) idt += gamma*VcDust(RHO,k,j,i); + InvDt(k,j,i) += idt; + }); + idfx::popRegion(); +} +// +// $ p_g^{(n+1)}=p_g^{(n)}+\sum_i\frac{\rho_g\gamma_i dt}{1+\rho_g \gamma_i dt}p_i^{(n)} $ +// We accumulate in the array "prefactor" +// $ \sum_i\frac{\rho_i\gamma_i dt}{1+\rho_g\gamma_i dt} +// + +void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) { + if(!feedback) { + // no feedback, no need for anything + return; + } + idfx::pushRegion("AddImplicitFluidMomentum"); + + auto UcGas = this->UcGas; + auto VcGas = this->VcGas; + auto VcDust = this->VcDust; + auto UcDust = this->UcDust; + + bool isFirst = this->instanceNumber == 0; + + + if(!implicit) { + IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); + } + + auto gammaDrag = this->gammaDrag; + gammaDrag.RefreshUserDrag(data); + + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) + // Where gamma is computed according to the choice of drag type + idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + real gamma = gammaDrag.GetGamma(k,j,i); // The drag coefficient + + + const real factor = UcDust(RHO,k,j,i)*gamma*dt/(1+UcGas(RHO,k,j,i)*gamma*dt); + if(isFirst) { + preFactor(k,j,i) = factor; + } else { + preFactor(k,j,i) += factor; + } + + for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { + UcGas(n,k,j,i) += dt * gamma * UcGas(RHO,k,j,i) * UcDust(n,k,j,i) / + (1 + UcGas(RHO,k,j,i)*dt*gamma); + } + }); + idfx::popRegion(); +} + +void Drag::NormalizeImplicitBackReaction(const real dt) { + if(!feedback) { + // no feedback, no need for anything + return; + } + idfx::pushRegion("AddImplicitFluidMomentum"); + + auto UcGas = this->UcGas; + auto preFactor = this->implicitFactor; + + if(!implicit) { + IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); + } + + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) + // Where gamma is computed according to the choice of drag type + idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + const real factor = 1+preFactor(k,j,i); + for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { + UcGas(n,k,j,i) /= factor; + } + }); + idfx::popRegion(); +} + +void Drag::AddImplicitFluidMomentum(const real dt) { + idfx::pushRegion("AddImplicitFluidMomentum"); + + auto UcGas = this->UcGas; + auto VcGas = this->VcGas; + auto UcDust = this->UcDust; + auto VcDust = this->VcDust; + auto InvDt = this->InvDt; + + bool feedback = this->feedback; + + if(!implicit) { + IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); + } + + auto gammaDrag = this->gammaDrag; + if(!feedback) gammaDrag.RefreshUserDrag(data); + + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) + // Where gamma is computed according to the choice of drag type + idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + real gamma = gammaDrag.GetGamma(k,j,i); // The drag coefficient + + for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { + real oldUc = UcDust(n,k,j,i); + UcDust(n,k,j,i) = (oldUc + dt * gamma * UcDust(RHO,k,j,i) * UcGas(n,k,j,i)) / + (1 + UcGas(RHO,k,j,i)*dt*gamma); + #if HAVE_ENERGY == 1 + real dp = UcDust(n,k,j,i) - oldUc; + // We add back the energy dissipated for the dust which is not accounted for // (since there is no energy equation for dust grains) // TODO(GL): this should be disabled in the case of a true multifluid system where // both fluids have a proper energy equation - UcGas(ENG,k,j,i) += dp*dv*VcDust(n,k,j,i); + if(feedback) UcGas(ENG,k,j,i) -= dp*VcDust(n,k,j,i); #endif } - // Cfl constraint - real idt = gamma*VcGas(RHO,k,j,i); - if(feedback) idt += gamma*VcDust(RHO,k,j,i); - InvDt(k,j,i) += idt; }); idfx::popRegion(); } void Drag::ShowConfig() { idfx::cout << "Drag: Using "; - switch(type) { - case Type::Gamma: + if(implicit) { + idfx::cout << " IMPLICIT "; + } else { + idfx::cout << " EXPLICIT "; + } + switch(gammaDrag.type) { + case GammaDrag::Type::Gamma: idfx::cout << "constant gamma"; break; - case Type::Tau: + case GammaDrag::Type::Tau: idfx::cout << "constant stopping time"; break; - case Type::Size: + case GammaDrag::Type::Size: idfx::cout << "constant dust size"; break; - case Type::Userdef: + case GammaDrag::Type::Userdef: idfx::cout << "user-defined"; break; } @@ -105,8 +206,60 @@ void Drag::ShowConfig() { } void Drag::EnrollUserDrag(UserDefDragFunc func) { + gammaDrag.EnrollUserDrag(func); +} + +//////////////////////////////////////////// +// GammaDrag function definitions +//////////////////////////////////////////// + +void GammaDrag::RefreshUserDrag(DataBlock *data) { + if(type == Type::Userdef) { + if(userDrag != NULL) { + idfx::pushRegion("GammaDrag::UserDrag"); + userDrag(data, instanceNumber, dragCoeff, gammai); + idfx::popRegion(); + } else { + IDEFIX_ERROR("No User-defined drag function has been enrolled"); + } + } +} + +void GammaDrag::EnrollUserDrag(UserDefDragFunc func) { if(type != Type::Userdef) { IDEFIX_ERROR("User-defined drag function requires drag entry to be set to \"userdef\""); } this->userDrag = func; } + +GammaDrag::GammaDrag(Input &input, std::string BlockName, int instanceNumber, DataBlock *data) { + if(input.CheckEntry(BlockName,"drag")>=0) { + std::string dragType = input.Get(BlockName,"drag",0); + if(dragType.compare("gamma") == 0) { + this->type = Type::Gamma; + } else if(dragType.compare("tau") == 0) { + this->type = Type::Tau; + this->VcGas = data->hydro->Vc; + } else if(dragType.compare("size") == 0) { + this->type = Type::Size; + this->eos = *(data->hydro->eos.get()); + this->VcGas = data->hydro->Vc; + } else if(dragType.compare("userdef") == 0) { + this->type = Type::Userdef; + this->gammai = IdefixArray3D("UserDrag", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } else { + std::stringstream msg; + msg << "Unknown drag type \"" << dragType + << "\" in your input file." << std::endl + << "Allowed values are: gamma, tau, epstein, stokes, userdef." << std::endl; + + IDEFIX_ERROR(msg); + } + } + // Fetch the drag coefficient for the current specie. + this->dragCoeff = input.Get(BlockName,"drag",instanceNumber+1); + this->instanceNumber = instanceNumber; +} diff --git a/src/fluid/drag.hpp b/src/fluid/drag.hpp index d5f7a5a73..29eca3097 100644 --- a/src/fluid/drag.hpp +++ b/src/fluid/drag.hpp @@ -13,38 +13,88 @@ #include "fluid_defs.hpp" #include "eos.hpp" -using UserDefDragFunc = void (*) (DataBlock *, real beta, IdefixArray3D &gammai); +using UserDefDragFunc = void (*) (DataBlock *, int n, real beta, IdefixArray3D &gammai); + +/* A class that holds the kind of drag law we intend to use */ +class GammaDrag { + public: + enum class Type{Gamma, Tau, Size, Userdef}; + GammaDrag() = default; + GammaDrag(Input &, std::string BlockName, int instanceNumber, DataBlock *data); + void EnrollUserDrag(UserDefDragFunc); + void RefreshUserDrag(DataBlock *); + + KOKKOS_INLINE_FUNCTION real GetGamma(const int k, const int j, const int i) const { + real gamma; // The drag coefficient + if(type == Type::Gamma) { + gamma = dragCoeff; + + } else if(type == Type::Tau) { + // In this case, the coefficient is the stopping time (assumed constant) + gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); + } else if(type == Type::Size) { + real cs; + // 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))); + #else + cs = eos.GetWaveSpeed(k,j,i); + #endif + gamma = cs/dragCoeff; + } else if(type == Type::Userdef) { + gamma = gammai(k,j,i); + } + return gamma; + } + + Type type; + real dragCoeff; + EquationOfState eos; + IdefixArray3D gammai; + IdefixArray4D VcGas; + + int instanceNumber; + UserDefDragFunc userDrag{NULL}; +}; class Drag { public: - enum class Type{Gamma, Tau, Size, Userdef}; // Different types of implementation for the drag force. template Drag(Input &, Fluid *); void ShowConfig(); // print configuration void AddDragForce(const real); + + ////////////////////////// + // Implicit functions + void AddImplicitBackReaction(const real, IdefixArray3D); // Add the back reaction + void NormalizeImplicitBackReaction(const real); // Normalize the implicit back reaction + void AddImplicitFluidMomentum(const real); // Add the implicit drag force on dust grains + ///////////////////////// + void EnrollUserDrag(UserDefDragFunc); // User defined drag function enrollment + bool IsImplicit() const { return implicit; } // Check if the drag is implicit IdefixArray4D UcDust; // Dust conservative quantities IdefixArray4D UcGas; // Gas conservative quantities IdefixArray4D VcDust; // Gas primitive quantities IdefixArray4D VcGas; // Gas primitive quantities IdefixArray3D InvDt; // The InvDt of current dust specie - IdefixArray3D gammai; // the drag coefficient (only used for user-defined dust grains) - Type type; + IdefixArray3D implicitFactor; // The prefactor used by the implicit timestepping + + GammaDrag gammaDrag; // The drag law private: DataBlock* data; - real dragCoeff; bool feedback{false}; - - UserDefDragFunc userDrag{NULL}; - - // Sound speed computation - EquationOfState *eos; + bool implicit{false}; + int instanceNumber; }; + #include "fluid.hpp" template @@ -58,45 +108,30 @@ Drag::Drag(Input &input, Fluid *hydroin): // Save the parent hydro object this->data = hydroin->data; - // We use the EOS of the basic fluid instance, as there is no EOS for dust! - this->eos = hydroin->data->hydro->eos.get(); // Check in which block we should fetch our information - std::string BlockName; + std::string blockName; if(Phys::dust) { - BlockName = "Dust"; + blockName = "Dust"; } else { IDEFIX_ERROR("Drag is currently implemented only for dusty fuids"); } - if(input.CheckEntry(BlockName,"drag")>=0) { - std::string dragType = input.Get(BlockName,"drag",0); - if(dragType.compare("gamma") == 0) { - this->type = Type::Gamma; - } else if(dragType.compare("tau") == 0) { - this->type = Type::Tau; - } else if(dragType.compare("size") == 0) { - this->type = Type::Size; - } else if(dragType.compare("userdef") == 0) { - this->type = Type::Userdef; - this->gammai = IdefixArray3D("UserDrag", + if(input.CheckEntry(blockName,"drag")>=0) { + this->instanceNumber = hydroin->instanceNumber; + this->gammaDrag = GammaDrag(input,blockName,instanceNumber,data); + + // Feedback is true by default, but can be switched off. + this->feedback = input.GetOrSet(blockName,"drag_feedback",0,true); + this->implicit = input.GetOrSet(blockName,"drag_implicit",0,false); + + if(implicit && instanceNumber == 0) { + this->implicitFactor = IdefixArray3D("ImplicitFactor", data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); - } else { - std::stringstream msg; - msg << "Unknown drag type \"" << dragType - << "\" in your input file." << std::endl - << "Allowed values are: gamma, tau, epstein, stokes, userdef." << std::endl; - - IDEFIX_ERROR(msg); } - // Fetch the drag coefficient for the current specie. - const int n = hydroin->instanceNumber; - this->dragCoeff = input.Get(BlockName,"drag",n+1); - // Feedback is true by default, but can be switched off. - this->feedback = input.GetOrSet(BlockName,"drag_feedback",0,true); } else { IDEFIX_ERROR("A [Drag] block is required in your input file to define the drag force."); } diff --git a/src/fluid/evolveStage.hpp b/src/fluid/evolveStage.hpp index 06754d0b7..31dae5b7d 100644 --- a/src/fluid/evolveStage.hpp +++ b/src/fluid/evolveStage.hpp @@ -61,7 +61,11 @@ void Fluid::EvolveStage(const real t, const real dt) { if(haveSourceTerms) AddSourceTerms(t, dt); // Step 5: add drag when needed - if(haveDrag) drag->AddDragForce(dt); + if(haveDrag) { + if(!drag->IsImplicit()) { + drag->AddDragForce(dt); + } + } if constexpr(Phys::mhd) { #if DIMENSIONS >= 2 diff --git a/src/fluid/fluid.hpp b/src/fluid/fluid.hpp index ae809b489..44cace6bf 100644 --- a/src/fluid/fluid.hpp +++ b/src/fluid/fluid.hpp @@ -585,52 +585,55 @@ Fluid::Fluid(Grid &grid, Input &input, DataBlock *datain, int n) { } for(int i = 0 ; i < Phys::nvar+nTracer ; i++) { - switch(i) { - case RHO: - VcName.push_back("RHO"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-RHO", RHO); - break; - case VX1: - VcName.push_back("VX1"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX1", VX1, IDIR); - break; - case VX2: - VcName.push_back("VX2"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX2", VX2, JDIR); - break; - case VX3: - VcName.push_back("VX3"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX3", VX3, KDIR); - break; - case BX1: - VcName.push_back("BX1"); - // never save cell-centered BX1 in dumps - break; - case BX2: - VcName.push_back("BX2"); - #if DIMENSIONS < 2 - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX2", BX2, JDIR); - #endif - break; - case BX3: - VcName.push_back("BX3"); - #if DIMENSIONS < 3 - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX3", BX3, KDIR); - #endif - break; - case PRS: - VcName.push_back("PRS"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-PRS", PRS); - break; - default: - if(i>=Phys::nvar) { - std::string tracerLabel = std::string("TR")+std::to_string(i-Phys::nvar); // ="TRn" - VcName.push_back(tracerLabel); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-"+tracerLabel, i); - } else { + if(i < Phys::nvar) { + // These are standard variable names + switch(i) { + case RHO: + VcName.push_back("RHO"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-RHO", RHO); + break; + case VX1: + VcName.push_back("VX1"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX1", VX1, IDIR); + break; + case VX2: + VcName.push_back("VX2"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX2", VX2, JDIR); + break; + case VX3: + VcName.push_back("VX3"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX3", VX3, KDIR); + break; + case BX1: + VcName.push_back("BX1"); + // never save cell-centered BX1 in dumps + break; + case BX2: + VcName.push_back("BX2"); + #if DIMENSIONS < 2 + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX2", BX2, JDIR); + #endif + break; + case BX3: + VcName.push_back("BX3"); + #if DIMENSIONS < 3 + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX3", BX3, KDIR); + #endif + break; + case PRS: + VcName.push_back("PRS"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-PRS", PRS); + break; + default: + // unknown variable name (this should never happen, but who knows...) VcName.push_back("Vc-"+std::to_string(i)); data->vtk->RegisterVariable(Vc, outputPrefix+"Vc-"+std::to_string(i), i); - } + break; + } + } else { //if(i>=Phys::nvar) + std::string tracerLabel = std::string("TR")+std::to_string(i-Phys::nvar); // ="TRn" + VcName.push_back(tracerLabel); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-"+tracerLabel, i); } data->vtk->RegisterVariable(Vc, outputPrefix+VcName[i], i); #ifdef WITH_HDF5 diff --git a/src/fluid/fluid_defs.hpp b/src/fluid/fluid_defs.hpp index 2269fc803..3f68c9b49 100644 --- a/src/fluid/fluid_defs.hpp +++ b/src/fluid/fluid_defs.hpp @@ -8,6 +8,7 @@ #ifndef FLUID_FLUID_DEFS_HPP_ #define FLUID_FLUID_DEFS_HPP_ +#include #include "../idefix.hpp" @@ -25,7 +26,7 @@ class Fluid; // Parabolic terms can have different status -enum HydroModuleStatus {Disabled, Constant, UserDefFunction }; +enum HydroModuleStatus {Disabled, Constant, UserDefFunction}; // Structure to describe the status of parabolic modules struct ParabolicModuleStatus { @@ -49,8 +50,9 @@ using SrcTermFunc = void (*) (Fluid *, const real t, const real dt); using EmfBoundaryFunc = void (*) (DataBlock &, const real t); using DiffusivityFunc = void (*) (DataBlock &, const real t, IdefixArray3D &); -using BragDiffusivityFunc = void (*) (DataBlock &, const real t, - IdefixArray3D &, IdefixArray3D &); + +using BragDiffusivityFunc = void (*) (DataBlock &, const real, + std::vector> &); // Deprecated signatures using SrcTermFuncOld = void (*) (DataBlock &, const real t, const real dt); diff --git a/src/global.hpp b/src/global.hpp index 8651f9748..265699b29 100644 --- a/src/global.hpp +++ b/src/global.hpp @@ -11,6 +11,7 @@ #include #include #include "arrays.hpp" +#include "npy.hpp" namespace idfx { int initialize(); // Initialisation routine for idefix @@ -44,6 +45,20 @@ IdefixArray1D ConvertVectorToIdefixArray(std::vector &inputVector) { return(outArr); } +///< dump Idefix array to a numpy array on disk +template +void DumpArray(std::string filename, ArrayType array) { + auto hArray = Kokkos::create_mirror(array); + Kokkos::deep_copy(hArray, array); + + std::array shape; + bool fortran_order{false}; + for (size_t i = 0; i < ArrayType::rank; ++i) { + shape[i] = array.extent(i); + } + npy::SaveArrayAsNumpy(filename, fortran_order, ArrayType::rank, shape.data(), hArray.data()); +} + } // namespace idfx class idfx::IdefixOutStream { diff --git a/src/loop.hpp b/src/loop.hpp index 8ad345dea..42192fa84 100644 --- a/src/loop.hpp +++ b/src/loop.hpp @@ -34,6 +34,12 @@ typedef Kokkos::TeamPolicy<>::member_type member_type; // Check if the user requested a specific loop unrolling strategy #if defined(LOOP_PATTERN_SIMD) + // Check that Idefix Arrays can be assigned from SIMD loop + static_assert( + Kokkos::SpaceAccessibility::accessible, + "Idefix Arrays cannot be accessed from SIMD loop. You should try another loop pattern." + ); constexpr LoopPattern defaultLoop = LoopPattern::SIMDFOR; #elif defined(LOOP_PATTERN_1DRANGE) constexpr LoopPattern defaultLoop = LoopPattern::RANGE; @@ -51,7 +57,16 @@ typedef Kokkos::TeamPolicy<>::member_type member_type; constexpr LoopPattern defaultLoop = LoopPattern::RANGE; #elif defined(KOKKOS_ENABLE_HIP) constexpr LoopPattern defaultLoop = LoopPattern::RANGE; + #elif defined(KOKKOS_ENABLE_SYCL) + constexpr LoopPattern defaultLoop = LoopPattern::RANGE; #elif defined(KOKKOS_ENABLE_SERIAL) + // Check that Idefix Arrays can be assigned from SIMD loop + static_assert( + Kokkos::SpaceAccessibility::accessible, + "Idefix Arrays cannot be accessed from Host, but Device is unknown/untested. " + "Ask the developers to add support." + ); constexpr LoopPattern defaultLoop = LoopPattern::SIMDFOR; #else #warning "Unknown target architeture: default to MDrange" diff --git a/src/main.cpp b/src/main.cpp index 6b65c4f1c..d21793d1d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -9,7 +9,7 @@ //@HEADER // ************************************************************************ // -// IDEFIX v 2.2.00 +// IDEFIX v 2.2.01 // // ************************************************************************ //@HEADER diff --git a/src/output/output.cpp b/src/output/output.cpp index a67cb1666..9df188926 100644 --- a/src/output/output.cpp +++ b/src/output/output.cpp @@ -55,8 +55,13 @@ Output::Output(Input &input, DataBlock &data) dumpPeriod = input.Get("Output","dmp",0); dumpLast = data.t - dumpPeriod; // dump something in the next CheckForWrite() if(input.CheckEntry("Output","dmp")>1) { - dumpTimePeriod = input.Get("Output","dmp",1); - std::string dumpTimeExtension = input.Get("Output","dmp",2); + std::string dumpString = input.Get("Output","dmp",1); + try { + dumpTimePeriod = std::stod(dumpString.substr(0, dumpString.size()-1), NULL); + } catch(const std::exception& e) { + IDEFIX_ERROR("The dump time period should be a number followed by a unit (s, m, h or d)"); + } + std::string dumpTimeExtension = dumpString.substr(dumpString.size()-1,1); if(dumpTimeExtension.compare("s")==0) { dumpTimePeriod *= 1.0; // Dump time period is in seconds by default } else if (dumpTimeExtension.compare("m")==0) { diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp index a2f6734f9..c8dd8cb96 100644 --- a/src/rkl/rkl.hpp +++ b/src/rkl/rkl.hpp @@ -762,32 +762,42 @@ void RKLegendre::CalcParabolicRHS(real t) { data->beg[IDIR],data->end[IDIR]+ioffset, KOKKOS_LAMBDA (int n, int k, int j, int i) { real Ax = A(k,j,i); - -#if GEOMETRY != CARTESIAN - if(Ax= 2) \ - || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) - if(dir==IDIR && nv==iMPHI) { - // Conserve angular momentum, hence flux is R*Vphi - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - } -#endif // GEOMETRY==POLAR OR CYLINDRICAL + #if (GEOMETRY == POLAR && COMPONENTS >= 2) \ + || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) + if(dir==IDIR) { + if(nv==iMPHI) { + Ax *= FABS(x1m(i)); + } + if constexpr(Phys::mhd) { + if(nv==iBPHI) { + // No area for this one + Ax = 1; + } + } + } + #endif // GEOMETRY==POLAR OR CYLINDRICAL + + #if GEOMETRY == SPHERICAL + if(dir==IDIR && nv==VX3) { + Ax *= FABS(x1m(i)); + } else if(dir==JDIR && nv==VX3) { + Ax *= FABS(sm(j)); + } + + if constexpr(Phys::mhd) { + if(dir == IDIR && (nv==BX3 || nv == BX2)) { + Ax = x1m(i); + } + if(dir==JDIR && nv==BX3) { + Ax = 1.0; + } + } + #endif // GEOMETRY == SPHERICAL -#if GEOMETRY == SPHERICAL && COMPONENTS == 3 - if(dir==IDIR && nv==iMPHI) { - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - } else if(dir==JDIR && nv==iMPHI) { - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(sm(j)); - } -#endif // GEOMETRY == SPHERICAL && COMPONENTS == 3 + Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax; } ); @@ -818,17 +828,42 @@ void RKLegendre::CalcParabolicRHS(real t) { } #if GEOMETRY != CARTESIAN - #ifdef iMPHI - if((dir==IDIR) && (nv == iMPHI)) { + if(dir==IDIR) { + #ifdef iMPHI + if((nv == iMPHI)) { rhs /= x1(i); } + #endif // iMPHI + real dx_ = dx(i); + real x1_ = x1(i); + + if constexpr(Phys::mhd) { + #if (GEOMETRY == POLAR || GEOMETRY == CYLINDRICAL) && (defined iBPHI) + if(nv==iBPHI) rhs = - 1 / dx_ * (Flux(iBPHI, k, j, i+1) - Flux(iBPHI, k, j, i) ); + + #elif (GEOMETRY == SPHERICAL) + real q = 1 / (x1_*dx_); + if(nv == BX2 || nv == BX3) { + rhs = -q * ((Flux(nv, k, j, i+1) - Flux(nv, k, j, i) )); + } + #endif // GEOMETRY + } // MHD + } // dir==IDIR + if(dir==JDIR) { #if (GEOMETRY == SPHERICAL) && (COMPONENTS == 3) - if((dir==JDIR) && (nv == iMPHI)) { + if(nv == iMPHI) { rhs /= FABS(s(j)); } + real dx_ = dx(j); + real rt_ = rt(i); + if constexpr(Phys::mhd) { + if(nv == iBPHI) { + rhs = - 1 / (rt_*dx_) * (Flux(nv, k, j+1, i) - Flux(nv, k, j, i)); + } + } #endif // GEOMETRY - // Nothing for KDIR - #endif // iMPHI + } // dir==JDIR + // Nothing for KDIR #endif // GEOMETRY != CARTESIAN // store the field components diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index b36f8a08e..afb4f57d4 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -5,4 +5,6 @@ target_sources(idefix PUBLIC ${CMAKE_CURRENT_LIST_DIR}/dumpImage.cpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/dumpImage.hpp PUBLIC ${CMAKE_CURRENT_LIST_DIR}/lookupTable.hpp + PUBLIC ${CMAKE_CURRENT_LIST_DIR}/column.cpp + PUBLIC ${CMAKE_CURRENT_LIST_DIR}/column.hpp ) diff --git a/src/utils/column.cpp b/src/utils/column.cpp new file mode 100644 index 000000000..34d55f56b --- /dev/null +++ b/src/utils/column.cpp @@ -0,0 +1,240 @@ +// *********************************************************************************** +// 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 +#include + +#include "column.hpp" +#include "idefix.hpp" +#include "input.hpp" +#include "output.hpp" +#include "grid.hpp" +#include "dataBlock.hpp" +#include "dataBlockHost.hpp" + +Column::Column(int dir, int sign, DataBlock *data) + : direction(dir), sign(sign) { + idfx::pushRegion("Column::Column"); + this->np_tot = data->np_tot; + this->np_int = data->np_int; + this->beg = data->beg; + + if(dir>= DIMENSIONS || dir < IDIR) IDEFIX_ERROR("Unknown direction for Column constructor"); + + // Allocate the array on which we do the average + this->ColumnArray = IdefixArray3D("ColumnArray",np_tot[KDIR], np_tot[JDIR], np_tot[IDIR]); + + // Elementary volumes and area + this->Volume = data->dV; + this->Area = data->A[dir]; + + // allocate helper array + if(dir == IDIR) { + localSum = IdefixArray2D("localSum",np_tot[KDIR], np_tot[JDIR]); + } + if(dir == JDIR) { + localSum = IdefixArray2D("localSum",np_tot[KDIR], np_tot[IDIR]); + } + if(dir == KDIR) { + localSum = IdefixArray2D("localSum",np_tot[JDIR], np_tot[IDIR]); + } + #ifdef WITH_MPI + // Create sub-MPI communicator dedicated to scan + int remainDims[3] = {false, false, false}; + remainDims[dir] = true; + MPI_Cart_sub(data->mygrid->CartComm, remainDims, &this->ColumnComm); + MPI_Comm_rank(this->ColumnComm, &this->MPIrank); + MPI_Comm_size(this->ColumnComm, &this->MPIsize); + + // create MPI class for boundary Xchanges + int ntarget = 0; + std::vector mapVars; + mapVars.push_back(ntarget); + + this->mpi.Init(data->mygrid, mapVars, data->nghost.data(), data->np_int.data()); + this->nproc = data->mygrid->nproc; + #endif + idfx::popRegion(); +} + +void Column::ComputeColumn(IdefixArray4D in, const int var) { + idfx::pushRegion("Column::ComputeColumn"); + const int nk = np_int[KDIR]; + const int nj = np_int[JDIR]; + const int ni = np_int[IDIR]; + const int kb = beg[KDIR]; + const int jb = beg[JDIR]; + const int ib = beg[IDIR]; + const int ke = kb+nk; + const int je = jb+nj; + const int ie = ib+ni; + + const int direction = this->direction; + auto column = this->ColumnArray; + auto dV = this->Volume; + auto A = this->Area; + auto localSum = this->localSum; + + if(direction==IDIR) { + // Inspired from loop.hpp + Kokkos::parallel_for("ColumnX1", team_policy (nk*nj, Kokkos::AUTO), + KOKKOS_LAMBDA (member_type team_member) { + int k = team_member.league_rank() / nj; + int j = team_member.league_rank() - k*nj + jb; + k += kb; + Kokkos::parallel_scan(Kokkos::TeamThreadRange<>(team_member,ib,ie), + [=] (int i, real &partial_sum, bool is_final) { + partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k,j,i+1))); + if(is_final) column(k,j,i) = partial_sum; + //partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k,j,i+1))); + }); + }); + } + if(direction==JDIR) { + // Inspired from loop.hpp + Kokkos::parallel_for("ColumnX2", team_policy (nk*ni, Kokkos::AUTO), + KOKKOS_LAMBDA (member_type team_member) { + int k = team_member.league_rank() / ni; + int i = team_member.league_rank() - k*ni + ib; + k += kb; + Kokkos::parallel_scan(Kokkos::TeamThreadRange<>(team_member,jb,je), + [=] (int j, real &partial_sum, bool is_final) { + partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k,j+1,i))); + if(is_final) column(k,j,i) = partial_sum; + }); + }); + } + if(direction==KDIR) { + // Inspired from loop.hpp + Kokkos::parallel_for("ColumnX2", team_policy (nj*ni, Kokkos::AUTO), + KOKKOS_LAMBDA (member_type team_member) { + int j = team_member.league_rank() / ni; + int i = team_member.league_rank() - j*ni + ib; + j += jb; + Kokkos::parallel_scan(Kokkos::TeamThreadRange<>(team_member,kb,ke), + [=] (int k, real &partial_sum, bool is_final) { + partial_sum += in(var,k,j,i)*dV(k,j,i) / (0.5*(A(k,j,i)+A(k+1,j,i))); + if(is_final) column(k,j,i) = partial_sum; + }); + }); + } + + #ifdef WITH_MPI + // Load the current sum + int dst,src; + MPI_Cart_shift(this->ColumnComm,0, 1, &src, &dst); + int size = localSum.extent(0)*localSum.extent(1); + if(MPIrank>0) { + MPI_Status status; + // Get the cumulative sum from previous processes + Kokkos::fence(); + MPI_Recv(localSum.data(), size, realMPI, src, 20, ColumnComm, &status); + // Add this to our cumulative sum + idefix_for("Addsum",kb,ke,jb,je,ib,ie, + KOKKOS_LAMBDA(int k, int j, int i) { + if(direction == IDIR) column(k,j,i) += localSum(k,j); + if(direction == JDIR) column(k,j,i) += localSum(k,i); + if(direction == KDIR) column(k,j,i) += localSum(j,i); + }); + } // rank =0 + // Get the final sum + if(MPIrank < MPIsize-1) { + if(direction==IDIR) { + idefix_for("Loadsum",kb,ke,jb,je, + KOKKOS_LAMBDA(int k, int j) { + localSum(k,j) = column(k,j,ie-1); + }); + } + if(direction==JDIR) { + idefix_for("Loadsum",kb,ke,ib,ie, + KOKKOS_LAMBDA(int k, int i) { + localSum(k,i) = column(k,je-1,i); + }); + } + if(direction==KDIR) { + idefix_for("Loadsum",jb,je,ib,ie, + KOKKOS_LAMBDA(int j, int i) { + localSum(j,i) = column(ke-1,j,i); + }); + } + // And send it + Kokkos::fence(); + MPI_Send(localSum.data(),size, realMPI, dst, 20, ColumnComm); + } // MPIrank small enough + #endif + // If we need it backwards + if(sign<0) { + // Compute total cumulative sum in the last subdomain + if(direction == IDIR) { + idefix_for("Loadsum",kb,ke,jb,je, + KOKKOS_LAMBDA(int k, int j) { + localSum(k,j) = column(k,j,ie-1) + in(var,k,j,ie-1) * dV(k,j,ie-1) + / (0.5*(A(k,j,ie-1)+A(k,j,ie))); + }); + } + if(direction == JDIR) { + idefix_for("Loadsum",kb,ke,ib,ie, + KOKKOS_LAMBDA(int k, int i) { + localSum(k,i) = column(k,je-1,i) + in(var,k,je-1,i) * dV(k,je-1,i) + / (0.5*(A(k,je-1,i)+A(k,je,i))); + }); + } + if(direction == KDIR) { + idefix_for("Loadsum",jb,je,ib,ie, + KOKKOS_LAMBDA(int j, int i) { + localSum(j,i) = column(ke-1,j,i) + in(var,ke-1,j,i) * dV(ke-1,j,i) + / (0.5*(A(ke-1,j,i)+A(ke,j,i))); + }); + } + + #ifdef WITH_MPI + Kokkos::fence(); + MPI_Bcast(localSum.data(),size, realMPI, MPIsize-1, ColumnComm); + #endif + // All substract the local column from the full column + + idefix_for("Subsum",kb,ke,jb,je,ib,ie, + KOKKOS_LAMBDA(int k, int j, int i) { + if(direction==IDIR) column(k,j,i) = localSum(k,j)-column(k,j,i); + if(direction==JDIR) column(k,j,i) = localSum(k,i)-column(k,j,i); + if(direction==KDIR) column(k,j,i) = localSum(j,i)-column(k,j,i); + }); + } // sign <0 + // Xchange boundary elements when using MPI to ensure that column + // density in the ghost zones are coherent + #ifdef WITH_MPI + // Create a 4D array that contains our column data + IdefixArray4D arr4D(column.data(), 1, this->np_tot[KDIR], + this->np_tot[JDIR], + this->np_tot[IDIR]); + + for(int dir = 0 ; dir < DIMENSIONS ; dir++) { + // MPI Exchange data when needed + if(nproc[dir]>1) { + switch(dir) { + case 0: + this->mpi.ExchangeX1(arr4D); + break; + case 1: + this->mpi.ExchangeX2(arr4D); + break; + case 2: + this->mpi.ExchangeX3(arr4D); + break; + } + } + } + #endif + idfx::popRegion(); +} + +void Column::ComputeColumn(IdefixArray3D in) { + // 4D alias + IdefixArray4D arr4D(in.data(), 1, in.extent(0), in.extent(1), in.extent(2)); + return this->ComputeColumn(arr4D,0); +} diff --git a/src/utils/column.hpp b/src/utils/column.hpp new file mode 100644 index 000000000..43383cbe6 --- /dev/null +++ b/src/utils/column.hpp @@ -0,0 +1,74 @@ +// *********************************************************************************** +// 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 UTILS_COLUMN_HPP_ +#define UTILS_COLUMN_HPP_ + +#include +#include + +#include "idefix.hpp" +#include "dataBlock.hpp" + + + +// A class to implement a parralel cumulative sum +class Column { + public: + //////////////////////////////////////////////////////////////////////////////////// + /// @brief Constructor of a cumulative sum(=column density) from setup's arguments + /// @param dir direction along which the integration is performed + /// @param sign: +1 for an integration from left to right, -1 for an integration from right to + /// left i.e (backwards) + /////////////////////////////////////////////////////////////////////////////////// + Column(int dir, int sign, DataBlock *); + + /////////////////////////////////////////////////////////////////////////////////// + /// @brief Effectively compute integral from the input array in argument + /// @param in: 4D input array + /// @param variable: index of the variable along which we do the integral (since the + /// intput array is 4D) + /////////////////////////////////////////////////////////////////////////////////// + void ComputeColumn(IdefixArray4D in, int variable); + + /////////////////////////////////////////////////////////////////////////////////// + /// @brief Effectively compute integral from the input array in argument + /// @param in: 3D input array + /////////////////////////////////////////////////////////////////////////////////// + void ComputeColumn(IdefixArray3D in); + + /////////////////////////////////////////////////////////////////////////////////// + /// @brief Get a reference to the computed column density array + /////////////////////////////////////////////////////////////////////////////////// + IdefixArray3D GetColumn() { + return (this->ColumnArray); + } + + private: + IdefixArray3D ColumnArray; + int direction; // The direction along which the column is computed + int sign; // whether we integrate from the left or from the right + std::array np_tot; + std::array np_int; + std::array beg; + + IdefixArray3D Area; + IdefixArray3D Volume; + + IdefixArray2D localSum; + #ifdef WITH_MPI + Mpi mpi; // Mpi object when WITH_MPI is set + MPI_Comm ColumnComm; + int MPIrank; + int MPIsize; + + std::array nproc; // 3D size of the MPI cartesian geometry + + #endif +}; + +#endif // UTILS_COLUMN_HPP_ diff --git a/test/Dust/DustEnergy/idefix-implicit.ini b/test/Dust/DustEnergy/idefix-implicit.ini new file mode 100644 index 000000000..9ba3c781a --- /dev/null +++ b/test/Dust/DustEnergy/idefix-implicit.ini @@ -0,0 +1,36 @@ +# This test checks that the total energy (thermal+dust kinetic+gas kinetic) +# is effectively conserved when drag is present + +[Grid] +X1-grid 1 0.0 500 u 1.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 1.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +gamma 1.4 + +[Dust] +nSpecies 1 +drag tau 1.0 +drag_feedback yes +drag_implicit yes + +[Boundary] +X1-beg periodic +X1-end periodic +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +dmp 1.0 +analysis 0.01 +log 1000 diff --git a/test/Dust/DustEnergy/testme.py b/test/Dust/DustEnergy/testme.py index f3b413e83..7187d299b 100755 --- a/test/Dust/DustEnergy/testme.py +++ b/test/Dust/DustEnergy/testme.py @@ -15,7 +15,7 @@ def testMe(test): test.configure() test.compile() - inifiles=["idefix.ini"] + inifiles=["idefix.ini","idefix-implicit.ini"] # loop on all the ini files for this test for ini in inifiles: diff --git a/test/Dust/DustyShock/definitions.hpp b/test/Dust/DustyShock/definitions.hpp new file mode 100644 index 000000000..cdbe23d6e --- /dev/null +++ b/test/Dust/DustyShock/definitions.hpp @@ -0,0 +1,5 @@ +#define COMPONENTS 1 +#define DIMENSIONS 1 +#define ISOTHERMAL + +#define GEOMETRY CARTESIAN diff --git a/test/Dust/DustyShock/idefix-implicit.ini b/test/Dust/DustyShock/idefix-implicit.ini new file mode 100644 index 000000000..2d92a0a52 --- /dev/null +++ b/test/Dust/DustyShock/idefix-implicit.ini @@ -0,0 +1,36 @@ +# This test checks the behaviour of a dust sound shock +# following the 4 fluids test of Benitez-Llambay+ 2019 + +[Grid] +X1-grid 1 0.0 400 u 40.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 500.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +csiso constant 1.0 + +[Dust] +nSpecies 3 +drag userdef 1.0 3.0 5.0 +drag_feedback yes +drag_implicit yes + +[Boundary] +X1-beg userdef +X1-end userdef +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +dmp 500.0 +vtk 500.0 +log 1000 diff --git a/test/Dust/DustyShock/idefix.ini b/test/Dust/DustyShock/idefix.ini new file mode 100644 index 000000000..a98ebd097 --- /dev/null +++ b/test/Dust/DustyShock/idefix.ini @@ -0,0 +1,35 @@ +# This test checks the behaviour of a dust sound shock +# following the 4 fluids test of Benitez-Llambay+ 2019 + +[Grid] +X1-grid 1 0.0 400 u 40.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 500.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +csiso constant 1.0 + +[Dust] +nSpecies 3 +drag userdef 1.0 3.0 5.0 +drag_feedback yes + +[Boundary] +X1-beg userdef +X1-end userdef +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +dmp 500.0 +vtk 500.0 +log 1000 diff --git a/test/Dust/DustyShock/python/testidefix.py b/test/Dust/DustyShock/python/testidefix.py new file mode 100755 index 000000000..8fca15a3f --- /dev/null +++ b/test/Dust/DustyShock/python/testidefix.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Mar 5 11:29:41 2020 + +@author: glesur +""" + +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) +from pytools.vtk_io import readVTK +import argparse +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp + +parser = argparse.ArgumentParser() +parser.add_argument("-noplot", + default=False, + help="disable plotting", + action="store_true") + + +args, unknown=parser.parse_known_args() + +V=readVTK('../data.0001.vtk', geometry='cartesian') + +# Find where the shock starts +istart = np.argwhere(V.data['RHO'][:,0,0]>2.0)[0][0]-1 + +# Compute analytical solution +M=2 # Mach number +K=np.asarray([1,3,5])/M # Drag coefficients + +# (56) from Benitez-Llambay 2018 +def get_wg(Mach, wi): + coeff=[1,np.sum(wi-1)-1/Mach**2-1,1/Mach**2] + solution=np.roots(coeff) + return(np.min(solution)) + +# (55) from Benitez-Llambay 2018 +def rhs(t,wi,Mach,K): + wg = get_wg(Mach,wi) + return K*(wg-wi) + +wi0=np.asarray([1.0,1.0,1.0]) +x=V.x[istart:] +arguments=M,K + +sol = solve_ivp(rhs, t_span=[x[0],x[-1]], y0=wi0, t_eval=x,args=arguments) + +# compute gas velocity +wg=np.zeros(len(x)) +for i in range(0,len(wg)): + wg[i]=get_wg(M,sol.y[:,i]) + + +if(not args.noplot): + plt.figure(1) + plt.plot(V.x,V.data['RHO'][:,0,0],'o',markersize=4,color='tab:purple') + plt.plot(x,1/wg,color='tab:purple') + plt.plot(V.x,V.data['Dust0_RHO'][:,0,0],'o',markersize=4,color='tab:orange') + plt.plot(x,1/sol.y[0,:],color='tab:orange') + plt.plot(V.x,V.data['Dust1_RHO'][:,0,0],'o',markersize=4,color='tab:blue') + plt.plot(x,1/sol.y[1,:],color='tab:blue') + plt.plot(V.x,V.data['Dust2_RHO'][:,0,0],'o',markersize=4,color='tab:green') + plt.plot(x,1/sol.y[2,:],color='tab:green') + plt.xlim([0,10]) + plt.title('Density') + + plt.figure(2) + plt.plot(V.x,V.data['VX1'][:,0,0],'o',markersize=4,color='tab:purple') + plt.plot(x,M*wg,color='tab:purple') + plt.plot(V.x[::4],V.data['Dust0_VX1'][::4,0,0],'o',markersize=4,color='tab:orange') + plt.plot(x,M*sol.y[0,:],color='tab:orange') + plt.plot(V.x[::4],V.data['Dust1_VX1'][::4,0,0],'o',markersize=4,color='tab:blue') + plt.plot(x,M*sol.y[1,:],color='tab:blue') + plt.plot(V.x[::4],V.data['Dust2_VX1'][::4,0,0],'o',markersize=4,color='tab:green') + plt.plot(x,M*sol.y[2,:],color='tab:green') + + plt.xlim([0,10]) + plt.title('Velocity') + + plt.ioff() + plt.show() + +# Compute the error on the dust densities + +error=np.mean(np.fabs(V.data['Dust0_RHO'][istart:,0,0]-1/sol.y[0,:]) + np.fabs(V.data['Dust1_RHO'][istart:,0,0]-1/sol.y[1,:]) + np.fabs(V.data['Dust2_RHO'][istart:,0,0]-1/sol.y[2,:])) / 3 +print("Error=%e"%error) +if error<3.6e-2: + print("SUCCESS!") + sys.exit(0) +else: + print("FAILURE!") + sys.exit(1) diff --git a/test/Dust/DustyShock/setup.cpp b/test/Dust/DustyShock/setup.cpp new file mode 100644 index 000000000..f8632248a --- /dev/null +++ b/test/Dust/DustyShock/setup.cpp @@ -0,0 +1,99 @@ +#include "idefix.hpp" +#include "setup.hpp" + +#define FILENAME "timevol.dat" + +void MyDrag(DataBlock *data, int n, real beta, IdefixArray3D &gamma) { + // Compute the drag coefficient gamma from the input beta + auto VcGas = data->hydro->Vc; + auto VcDust = data->dust[n]->Vc; + + idefix_for("MyDrag",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + gamma(k,j,i) = beta/VcGas(RHO,k,j,i)/VcDust(RHO,k,j,i); + }); +} + + +void ApplyBoundary(DataBlock *data, IdefixArray4D Vc, int dir, BoundarySide side) { + if(dir==IDIR) { + int iref,ibeg,iend; + if(side == left) { + iref = data->beg[IDIR]; + ibeg = 0; + iend = data->beg[IDIR]; + } else { + iref =data->end[IDIR] - 1; + ibeg=data->end[IDIR]; + iend=data->np_tot[IDIR]; + } + int nvar = Vc.extent(0); + idefix_for("UserDefBoundary",0,data->np_tot[KDIR],0,data->np_tot[JDIR],ibeg,iend, + KOKKOS_LAMBDA (int k, int j, int i) { + for(int n = 0 ; n < nvar ; n++) { + Vc(n,k,j,i) = Vc(n,k,j,iref ); + } + }); + } +} + +void UserdefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) { + ApplyBoundary(hydro->data, hydro->Vc, dir, side); +} + +void UserdefBoundaryDust(Fluid *dust, int dir, BoundarySide side, real t) { + ApplyBoundary(dust->data, dust->Vc, dir, side); +} + + +// Initialisation routine. Can be used to allocate +// Arrays or variables which are used later on +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + + data.hydro->EnrollUserDefBoundary(&UserdefBoundary); + if(data.haveDust) { + int nSpecies = data.dust.size(); + for(int n = 0 ; n < nSpecies ; n++) { + data.dust[n]->EnrollUserDefBoundary(&UserdefBoundaryDust); + data.dust[n]->drag->EnrollUserDrag(&MyDrag); + } + } + +} + +// This routine initialize the flow +// Note that data is on the device. +// One can therefore define locally +// a datahost and sync it, if needed +void Setup::InitFlow(DataBlock &data) { + // Create a host copy + DataBlockHost d(data); + + const real x0 = 4.0; + + for(int k = 0; k < d.np_tot[KDIR] ; k++) { + for(int j = 0; j < d.np_tot[JDIR] ; j++) { + for(int i = 0; i < d.np_tot[IDIR] ; i++) { + real x = d.x[IDIR](i); + + d.Vc(RHO,k,j,i) = (x < x0) ? 1.0 : 16.0; + d.Vc(VX1,k,j,i) = (x < x0) ? 2.0 : 0.125; + + for(int n = 0 ; n < data.dust.size(); n++) { + d.dustVc[n](RHO,k,j,i) = (x < x0) ? 1.0 : 16.0; + d.dustVc[n](VX1,k,j,i) = (x < x0) ? 2.0 : 0.125; + } + + + } + } + } + + // Send it all, if needed + d.SyncToDevice(); +} + +// Analyse data to produce an output +void MakeAnalysis(DataBlock & data) { + +} diff --git a/test/Dust/DustyShock/testme.py b/test/Dust/DustyShock/testme.py new file mode 100755 index 000000000..2b03d1636 --- /dev/null +++ b/test/Dust/DustyShock/testme.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) + +import pytools.idfx_test as tst + +name="dump.0001.dmp" + +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini","idefix-implicit.ini"] + + # loop on all the 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,tolerance=1e-14) + + +test=tst.idfxTest() + +if not test.all: + if(test.check): + test.checkOnly(filename=name) + else: + testMe(test) +else: + test.noplot = True + test.reconstruction=2 + testMe(test) diff --git a/test/Dust/DustyWave/idefix-implicit.ini b/test/Dust/DustyWave/idefix-implicit.ini new file mode 100644 index 000000000..726585e0b --- /dev/null +++ b/test/Dust/DustyWave/idefix-implicit.ini @@ -0,0 +1,36 @@ +# This test checks the dissipation of a sound wave by a dust grains +# partially coupled to the gas (Riols & Lesur 2018, appendix A) + +[Grid] +X1-grid 1 0.0 500 u 1.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 10.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +csiso constant 1.0 + +[Dust] +nSpecies 1 +drag tau 1.0 +drag_feedback yes +drag_implicit yes + +[Boundary] +X1-beg periodic +X1-end periodic +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +dmp 10.0 +analysis 0.01 +log 1000 diff --git a/test/Dust/DustyWave/python/testidefix.py b/test/Dust/DustyWave/python/testidefix.py index c98bc8213..ad52d3306 100755 --- a/test/Dust/DustyWave/python/testidefix.py +++ b/test/Dust/DustyWave/python/testidefix.py @@ -35,9 +35,6 @@ # Get the minimal decay rate (this should be the one that pops up) tau=np.amax(np.imag(sol)) - - - # load the dat file produced by the setup raw=np.loadtxt('../timevol.dat',skiprows=1) t=raw[:,0] @@ -58,10 +55,10 @@ # Compute decay rate: tau_measured=t[-1]/(2*np.log(etot[-1]/etot[0])) # error on the decay rate: -error=(tau_measured-tau)/tau +error=np.abs((tau_measured-tau)/tau) print("error=%f"%error) -if(error<0.065): +if(error<0.07): print("Success!") else: print("Failure!") diff --git a/test/Dust/DustyWave/testme.py b/test/Dust/DustyWave/testme.py index b6a75ce3c..2b03d1636 100755 --- a/test/Dust/DustyWave/testme.py +++ b/test/Dust/DustyWave/testme.py @@ -15,7 +15,7 @@ def testMe(test): test.configure() test.compile() - inifiles=["idefix.ini"] + inifiles=["idefix.ini","idefix-implicit.ini"] # loop on all the ini files for this test for ini in inifiles: diff --git a/test/Dust/FargoPlanet/setup.cpp b/test/Dust/FargoPlanet/setup.cpp index af611911d..8447efe0b 100644 --- a/test/Dust/FargoPlanet/setup.cpp +++ b/test/Dust/FargoPlanet/setup.cpp @@ -42,7 +42,7 @@ void MySoundSpeed(DataBlock &data, const real t, IdefixArray3D &cs) { // a = 20 cm * beta * (Sigma_phys/100 g.cm^-2) / (rho_s / 2 g.cm^-3) // NB: checked against A. Johansen+ (2007) -void MyDrag(DataBlock *data, real beta, IdefixArray3D &gamma) { +void MyDrag(DataBlock *data, int nSpecie, real beta, IdefixArray3D &gamma) { // Compute the drag coefficient gamma from the input beta auto x1 = data->x[IDIR]; real sigma0 = sigma0Glob; diff --git a/test/MHD/AxisFluxTube/setup.cpp b/test/MHD/AxisFluxTube/setup.cpp index 34009da85..4c30506d3 100644 --- a/test/MHD/AxisFluxTube/setup.cpp +++ b/test/MHD/AxisFluxTube/setup.cpp @@ -112,7 +112,7 @@ void CoarsenFunction(DataBlock &data) { Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { output.EnrollUserDefVariables(&ComputeUserVars); data.hydro->EnrollUserDefBoundary(&UserdefBoundary); - output.EnrollAnalysis(&Analysis); + //output.EnrollAnalysis(&Analysis); Rtorus = input.Get("Setup","Rtorus",0); Ztorus = input.Get("Setup","Ztorus",0); Rin = input.Get("Setup","Rin",0); diff --git a/test/MHD/MTI/idefix-rkl.ini b/test/MHD/MTI/idefix-rkl.ini index 1b9b11eb2..09e9eb72d 100644 --- a/test/MHD/MTI/idefix-rkl.ini +++ b/test/MHD/MTI/idefix-rkl.ini @@ -12,7 +12,7 @@ nstages 3 [Hydro] solver hlld gamma 1.6666666666666666 -bragTDiffusion rkl nolimiter userdef +bragTDiffusion rkl nolimiter nosat userdef bragViscosity rkl nolimiter userdef [Gravity] diff --git a/test/MHD/MTI/idefix-sl.ini b/test/MHD/MTI/idefix-sl.ini index f2c9deb2b..aa707914f 100644 --- a/test/MHD/MTI/idefix-sl.ini +++ b/test/MHD/MTI/idefix-sl.ini @@ -12,7 +12,7 @@ nstages 3 [Hydro] solver hlld gamma 1.6666666666666666 -bragTDiffusion rkl mc userdef +bragTDiffusion rkl mc nosat userdef bragViscosity rkl vanleer userdef [Gravity] diff --git a/test/MHD/MTI/idefix.ini b/test/MHD/MTI/idefix.ini index 7a09c406c..a159c4f70 100644 --- a/test/MHD/MTI/idefix.ini +++ b/test/MHD/MTI/idefix.ini @@ -12,7 +12,7 @@ nstages 3 [Hydro] solver hlld gamma 1.6666666666666666 -bragTDiffusion explicit nolimiter userdef +bragTDiffusion explicit nolimiter nosat userdef bragViscosity explicit nolimiter userdef [Gravity] diff --git a/test/MHD/MTI/setup.cpp b/test/MHD/MTI/setup.cpp index e631f7b2b..552f76465 100644 --- a/test/MHD/MTI/setup.cpp +++ b/test/MHD/MTI/setup.cpp @@ -39,9 +39,13 @@ void Potential(DataBlock& data, const real t, IdefixArray1D& x1, IdefixArr } -void MyBragThermalConductivity(DataBlock &data, const real t, IdefixArray3D &kparArr, IdefixArray3D &knorArr) { +void MyBragThermalConductivity(DataBlock &data, const real t, std::vector> &userdefArr) { IdefixArray4D Vc = data.hydro->Vc; IdefixArray1D x2 = data.x[JDIR]; + + IdefixArray3D kparArr = userdefArr.at(0); + IdefixArray3D knorArr = userdefArr.at(1); + real ksi = ksiGlob; idefix_for("MyThConductivity",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], KOKKOS_LAMBDA (int k, int j, int i) { diff --git a/test/MHD/clessTDiffusion/CMakeLists.txt b/test/MHD/clessTDiffusion/CMakeLists.txt new file mode 100644 index 000000000..629aed2d1 --- /dev/null +++ b/test/MHD/clessTDiffusion/CMakeLists.txt @@ -0,0 +1 @@ +enable_idefix_property(Idefix_MHD) diff --git a/test/MHD/clessTDiffusion/definitions.hpp b/test/MHD/clessTDiffusion/definitions.hpp new file mode 100644 index 000000000..21b1f8733 --- /dev/null +++ b/test/MHD/clessTDiffusion/definitions.hpp @@ -0,0 +1,5 @@ +#define COMPONENTS 3 +#define DIMENSIONS 1 +#define GEOMETRY SPHERICAL + +#define SMALL_PRESSURE_TEMPERATURE (0.1) diff --git a/test/MHD/clessTDiffusion/idefix.ini b/test/MHD/clessTDiffusion/idefix.ini new file mode 100644 index 000000000..9c444b928 --- /dev/null +++ b/test/MHD/clessTDiffusion/idefix.ini @@ -0,0 +1,41 @@ +[Grid] +X1-grid 1 1.0 256 l 40.0 +X2-grid 1 1.4708 1 u 1.6708 +X3-grid 1 0.9892 1 u 1.0108 + +[TimeIntegrator] +CFL 0.4 +CFL_max_var 1.1 # not used +tstop 300.0 +first_dt 1.e-6 +nstages 2 + +[Hydro] +solver hlld +gamma 1.6666666666666 +bragTDiffusion rkl mc wcless userdef + +[Gravity] +potential central +Mcentral 1.0 + +[Boundary] +X1-beg userdef +X1-end outflow +X2-beg periodic +X2-end periodic +X3-beg periodic +X3-end periodic + +[Setup] +UNIT_DENSITY 1.6726e-16 +UNIT_LENGTH 6.9570e10 +UNIT_VELOCITY 4.3670e7 +cs_vesc 0.26 +va_vesc 0.3 +k0 9e-7 + +[Output] +vtk 100.0 +dmp 300.0 +log 10 diff --git a/test/MHD/clessTDiffusion/python/testidefix.py b/test/MHD/clessTDiffusion/python/testidefix.py new file mode 100644 index 000000000..9427aca82 --- /dev/null +++ b/test/MHD/clessTDiffusion/python/testidefix.py @@ -0,0 +1,40 @@ +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) +from pytools.vtk_io import readVTK +import numpy as np +import inifix + +conf = inifix.load("../idefix.ini") +gamma = conf["Hydro"]["gamma"] +kappa = conf["Hydro"]["bragTDiffusion"][-1] + +current_VTK = readVTK('../data.0003.vtk', geometry='spherical') + +rho = current_VTK.data['RHO'].squeeze() +prs = current_VTK.data['PRS'].squeeze() +r = current_VTK.r + +idx1=np.where(r > 20)[0][0] +idx2=np.where(r > 30)[0][0] + +def gamma_prime(gamma, beta): + return (gamma+beta*(gamma-1))/(1+beta*(gamma-1)) + +success = True +eps = 3e-4 + +gamma_eff = np.gradient(prs, r)/np.gradient(rho, r)*rho/prs + +Error=np.abs(gamma_eff[idx1:idx2]-gamma_prime(gamma, 1.5)).mean() +if Error > eps: + success=False + +if success: + print("SUCCESS") + print("Error: {0}".format(Error)) + sys.exit(0) +else: + print("Failed") + print("Error: {0}".format(Error)) + sys.exit(1) diff --git a/test/MHD/clessTDiffusion/setup.cpp b/test/MHD/clessTDiffusion/setup.cpp new file mode 100644 index 000000000..c540c0def --- /dev/null +++ b/test/MHD/clessTDiffusion/setup.cpp @@ -0,0 +1,158 @@ +#include "idefix.hpp" +#include "setup.hpp" + +real udenGlob; +real ulenGlob; +real uvelGlob; +real gammaGlob; +real cs_vescGlob; +real va_vescGlob; +real k0_Glob; +real ParkerWind(real); + +const real kB{1.3807e-16}; +const real mp{1.6726e-24}; + +void MyClessThermalConductivity(DataBlock &data, const real t, std::vector> &userdefArr) { + IdefixArray4D Vc = data.hydro->Vc; + IdefixArray1D x1 = data.x[IDIR]; + IdefixArray1D x2 = data.x[JDIR]; + + IdefixArray3D kparArr = userdefArr.at(0); + IdefixArray3D knorArr = userdefArr.at(1); + IdefixArray3D clessAlpha = userdefArr.at(2); + IdefixArray3D clessBeta = userdefArr.at(3); + + real norm = mp*0.5/(udenGlob*uvelGlob*ulenGlob*kB); + real uTemp=0.5*uvelGlob*uvelGlob*mp/kB; + real k0 = k0_Glob*norm; + idefix_for("MyClessThConductivity",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + kparArr(k,j,i) = k0*pow(Vc(PRS,k,j,i)/Vc(RHO,k,j,i)*uTemp,2.5); + knorArr(k,j,i) = 0.; + clessAlpha(k,j,i) = (1.0-tanh(x1(i)-10))/2; + clessBeta(k,j,i) = -1.5; + }); +} + +// User-defined boundaries +void UserDefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) { + + DataBlock *data = hydro->data; + + if( (dir==IDIR) && (side == left)) { + IdefixArray4D Vc = hydro->Vc; + IdefixArray4D Vs = hydro->Vs; + IdefixArray1D x1 = data->x[IDIR]; + + real rc,vwind0; + real cs=cs_vescGlob*sqrt(2.); + real va_vesc = va_vescGlob; + real mu = va_vesc * sqrt(2.); + real PonRho; + + PonRho = cs*cs; + rc = 0.25 / (cs_vescGlob*cs_vescGlob); + vwind0 = ParkerWind(1./rc) * cs; + + hydro->boundary->BoundaryFor("UserDefBoundary", dir, side, + KOKKOS_LAMBDA (int k, int j, int i) { + real r = x1(i); + + Vc(RHO,k,j,i) = vwind0/(vwind0 * r * r); + Vc(PRS,k,j,i) = PonRho * Vc(RHO, k, j, i); + Vc(VX1,k,j,i) = vwind0; + Vc(VX2,k,j,i) = 0.0; + Vc(VX3,k,j,i) = 0.0; + Vc(BX1,k,j,i) = mu / (r*r); + Vc(BX2,k,j,i) = 0.0; + Vc(BX3,k,j,i) = 0.0; + + }); + } +} + +// Initialisation routine. Can be used to allocate +// Arrays or variables which are used later on +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + // Set the function for userdefboundary + data.hydro->EnrollUserDefBoundary(&UserDefBoundary); + data.hydro->bragThermalDiffusion->EnrollBragThermalDiffusivity(&MyClessThermalConductivity); + gammaGlob=input.Get("Hydro", "gamma", 0); + udenGlob=input.Get("Setup", "UNIT_DENSITY",0); + ulenGlob=input.Get("Setup", "UNIT_LENGTH",0); + uvelGlob=input.Get("Setup", "UNIT_VELOCITY",0); + cs_vescGlob=input.Get("Setup", "cs_vesc", 0); + va_vescGlob=input.Get("Setup", "va_vesc", 0); + k0_Glob = input.Get("Setup", "k0", 0); +} + +// This routine initialize the flow +// Note that data is on the device. +// One can therefore define locally +// a datahost and sync it, if needed +void Setup::InitFlow(DataBlock &data) { + // Create a host copy + DataBlockHost d(data); + + real r,rl; + real PonRho, vwind0, rc; + real cs=cs_vescGlob*sqrt(2.); + + + rc = 0.25 / (cs_vescGlob*cs_vescGlob); + vwind0 = ParkerWind(1./rc) * cs; + PonRho = cs*cs; + real mu = va_vescGlob * sqrt(2.); + + for(int k = 0; k < d.np_tot[KDIR] ; k++) { + for(int j = 0; j < d.np_tot[JDIR] ; j++) { + for(int i = 0; i < d.np_tot[IDIR] ; i++) { + r=d.x[IDIR](i); + + real vwind; + + vwind = ParkerWind(r/rc) * cs; + + d.Vc(RHO,k,j,i) = 1.0*vwind0/(vwind * r * r); + d.Vc(PRS,k,j,i) = PonRho * d.Vc(RHO, k, j, i); + d.Vc(VX1,k,j,i) = vwind; + d.Vc(VX2,k,j,i) = 0.0; + d.Vc(VX3,k,j,i) = 0.0; + + rl=d.xl[IDIR](i); // Radius on the left side of the cell + d.Vs(BX1s, k, j, i) = mu / (rl*rl); + + } + } + } + // Send it all, if needed + d.SyncToDevice(); +} + +// Analyse data to produce an output +void MakeAnalysis(DataBlock & data) { +} + +/**************************************************/ +real ParkerWind(real x) +/* Parker wind velocity in unit of iso sound speed + x = radius / critical radius. +**************************************************/ +{ + real v, f; + real vref; + + v = 1e-7; + f = v*v-2*log(v)-4/x-4*log(x)+3; + if (x>1) {v=10;} + while (fabs(f) > 1e-10){ + vref = v; + v = v - 0.5*f/(v-1/v); + while (v < 0){ + v = (v + vref)/2; + } + f = v*v-2*log(v)-4/x-4*log(x)+3; + } + return v; +} diff --git a/test/MHD/clessTDiffusion/testme.py b/test/MHD/clessTDiffusion/testme.py new file mode 100755 index 000000000..6ed5b1648 --- /dev/null +++ b/test/MHD/clessTDiffusion/testme.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) + +import pytools.idfx_test as tst + + +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini"] + + # loop on all the ini files for this test + name = "dump.0001.dmp" + for ini in inifiles: + test.run(inputFile=ini) + test.standardTest() + if test.init: + test.makeReference(filename=name) + test.nonRegressionTest(filename=name, tolerance=2e-15) + + +test=tst.idfxTest() + +if not test.all: + if(test.check): + test.checkOnly(filename="dump.0001.dmp") + else: + testMe(test) +else: + test.noplot = True + test.vectPot=False + test.single=False + test.reconstruction=2 + test.mpi=False + testMe(test) diff --git a/test/MHD/sphBragTDiffusion/idefix.ini b/test/MHD/sphBragTDiffusion/idefix.ini index 35ac54e8b..795ae7bb2 100644 --- a/test/MHD/sphBragTDiffusion/idefix.ini +++ b/test/MHD/sphBragTDiffusion/idefix.ini @@ -13,7 +13,7 @@ nstages 3 [Hydro] solver hlld gamma 1.4 -bragTDiffusion rkl nolimiter constant 50.0 +bragTDiffusion rkl nolimiter nosat constant 50.0 [Setup] amplitude 1e-4 diff --git a/test/utils/columnDensity/definitions.hpp b/test/utils/columnDensity/definitions.hpp new file mode 100644 index 000000000..3b581bd62 --- /dev/null +++ b/test/utils/columnDensity/definitions.hpp @@ -0,0 +1,4 @@ +#define COMPONENTS 3 +#define DIMENSIONS 3 + +#define GEOMETRY CARTESIAN diff --git a/test/utils/columnDensity/idefix.ini b/test/utils/columnDensity/idefix.ini new file mode 100644 index 000000000..5e80c8bce --- /dev/null +++ b/test/utils/columnDensity/idefix.ini @@ -0,0 +1,24 @@ +[Grid] +X1-grid 1 0.0 64 u 1.0 +X2-grid 1 0.0 32 u 1.0 +X3-grid 1 0.0 16 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 0.0 +nstages 2 + +[Hydro] +solver hllc +gamma 1.4 + +[Boundary] +X1-beg periodic +X1-end periodic +X2-beg periodic +X2-end periodic +X3-beg periodic +X3-end periodic + +[Output] +analysis 0.01 diff --git a/test/utils/columnDensity/setup.cpp b/test/utils/columnDensity/setup.cpp new file mode 100644 index 000000000..fe7c57e69 --- /dev/null +++ b/test/utils/columnDensity/setup.cpp @@ -0,0 +1,169 @@ +#include "idefix.hpp" +#include "setup.hpp" +#include "column.hpp" + + +Column *columnX1Left; +Column *columnX1Right; +Column *columnX2Left; +Column *columnX2Right; +Column *columnX3Left; +Column *columnX3Right; + +// Analyse data to check that column density works as expected +void Analysis(DataBlock & data) { + + DataBlockHost d(data); + + // Try the 4D array interface + columnX1Left->ComputeColumn(data.hydro->Vc,RHO); + columnX1Right->ComputeColumn(data.hydro->Vc,RHO); + + // Try the 3D array interface + IdefixArray3D rho("rho",data.np_tot[KDIR],data.np_tot[JDIR],data.np_tot[IDIR]); + auto Vc = data.hydro->Vc; + + idefix_for("init rho",data.beg[KDIR],data.end[KDIR],data.beg[JDIR],data.end[JDIR],data.beg[IDIR],data.end[IDIR], + KOKKOS_LAMBDA(int k, int j, int i) { + rho(k,j,i) = Vc(RHO,k,j,i); + }); + columnX2Left->ComputeColumn(rho); + columnX2Right->ComputeColumn(rho); + columnX3Left->ComputeColumn(rho); + columnX3Right->ComputeColumn(rho); + + IdefixArray3D columnDensityLeft, columnDensityRight; + IdefixArray3D::HostMirror columnDensityLeftHost, columnDensityRightHost; + // IDIR + columnDensityLeft = columnX1Left->GetColumn(); + columnDensityRight = columnX1Right->GetColumn(); + columnDensityLeftHost = Kokkos::create_mirror_view(columnDensityLeft); + columnDensityRightHost = Kokkos::create_mirror_view(columnDensityRight); + Kokkos::deep_copy(columnDensityLeftHost,columnDensityLeft); + Kokkos::deep_copy(columnDensityRightHost,columnDensityRight); + + real errMax = 0.0; + + // Because this particular setup has a constant density =1, we expect + // the column density from the left to match the cell right coordinates, + // and the column density from the right to match L-the cell left coordinate (where L is the box size) + // here we have L=1 + // This routine checks that all of the available column densities do match the expected values. + for(int k = data.beg[KDIR]; k < data.end[KDIR] ; k++) { + for(int j = data.beg[JDIR]; j < data.end[JDIR] ; j++) { + for(int i = data.beg[IDIR]; i < data.end[IDIR] ; i++) { + real err = std::fabs(columnDensityLeftHost(k,j,i)-d.xr[IDIR](i)); + if(err>errMax) errMax=err; + err = std::fabs(columnDensityRightHost(k,j,i)-(1-d.xl[IDIR](i))); + if(err>errMax) errMax=err; + } + } + } + idfx::cout << "Error on column density in IDIR=" << std::scientific << errMax << std::endl; + if(errMax>1e-14) { + IDEFIX_ERROR("Error above tolerance"); + } + // JDIR + columnDensityLeft = columnX2Left->GetColumn(); + columnDensityRight = columnX2Right->GetColumn(); + columnDensityLeftHost = Kokkos::create_mirror_view(columnDensityLeft); + columnDensityRightHost = Kokkos::create_mirror_view(columnDensityRight); + Kokkos::deep_copy(columnDensityLeftHost,columnDensityLeft); + Kokkos::deep_copy(columnDensityRightHost,columnDensityRight); + + errMax = 0.0; + + for(int k = data.beg[KDIR]; k < data.end[KDIR] ; k++) { + for(int j = data.beg[JDIR]; j < data.end[JDIR] ; j++) { + for(int i = data.beg[IDIR]; i < data.end[IDIR] ; i++) { + real err = std::fabs(columnDensityLeftHost(k,j,i)-d.xr[JDIR](j)); + if(err>errMax) errMax=err; + err = std::fabs(columnDensityRightHost(k,j,i)-(1-d.xl[JDIR](j))); + if(err>errMax) errMax=err; + } + } + } + idfx::cout << "Error on column density in JDIR=" << std::scientific << errMax << std::endl; + if(errMax>1e-14) { + IDEFIX_ERROR("Error above tolerance"); + } + // KDIR + columnDensityLeft = columnX3Left->GetColumn(); + columnDensityRight = columnX3Right->GetColumn(); + columnDensityLeftHost = Kokkos::create_mirror_view(columnDensityLeft); + columnDensityRightHost = Kokkos::create_mirror_view(columnDensityRight); + Kokkos::deep_copy(columnDensityLeftHost,columnDensityLeft); + Kokkos::deep_copy(columnDensityRightHost,columnDensityRight); + + errMax = 0.0; + + for(int k = data.beg[KDIR]; k < data.end[KDIR] ; k++) { + for(int j = data.beg[JDIR]; j < data.end[JDIR] ; j++) { + for(int i = data.beg[IDIR]; i < data.end[IDIR] ; i++) { + real err = std::fabs(columnDensityLeftHost(k,j,i)-d.xr[KDIR](k)); + if(err>errMax) errMax=err; + err = std::fabs(columnDensityRightHost(k,j,i)-(1-d.xl[KDIR](k))); + if(err>errMax) errMax=err; + } + } + } + idfx::cout << "Error on column density in KDIR=" << std::scientific << errMax << std::endl; + if(errMax>1e-14) { + IDEFIX_ERROR("Error above tolerance"); + } + +} + +void InternalBoundary(Fluid * hydro, const real t) { + IdefixArray4D Vc = hydro->Vc; + idefix_for("InternalBoundary",0,hydro->data->np_tot[KDIR], + 0,hydro->data->np_tot[JDIR], + 0,hydro->data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + // Cancel any motion that could be happening + Vc(VX1,k,j,i) = 0.0; + Vc(VX2,k,j,i) = 0.0; + Vc(VX3,k,j,i) = 0.0; + }); +} + + +// Initialisation routine. Can be used to allocate +// Arrays or variables which are used later on +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + output.EnrollAnalysis(&Analysis); + data.hydro->EnrollInternalBoundary(&InternalBoundary); + + columnX1Left = new Column(IDIR, 1, &data); + columnX1Right = new Column(IDIR, -1, &data); + columnX2Left = new Column(JDIR, 1, &data); + columnX2Right = new Column(JDIR, -1, &data); + columnX3Left = new Column(KDIR, 1, &data); + columnX3Right = new Column(KDIR, -1, &data); + // Initialise the output file +} + +// This routine initialize the flow +// Note that data is on the device. +// One can therefore define locally +// a datahost and sync it, if needed +void Setup::InitFlow(DataBlock &data) { + // Create a host copy + DataBlockHost d(data); + + + for(int k = 0; k < d.np_tot[KDIR] ; k++) { + for(int j = 0; j < d.np_tot[JDIR] ; j++) { + for(int i = 0; i < d.np_tot[IDIR] ; i++) { + + d.Vc(RHO,k,j,i) = 1.0; + d.Vc(VX1,k,j,i) = ZERO_F; + d.Vc(PRS,k,j,i) = 1.0; + + } + } + } + + // Send it all, if needed + d.SyncToDevice(); +} diff --git a/test/utils/columnDensity/testme.py b/test/utils/columnDensity/testme.py new file mode 100755 index 000000000..6f017d7d8 --- /dev/null +++ b/test/utils/columnDensity/testme.py @@ -0,0 +1,23 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) +import pytools.idfx_test as tst + +test=tst.idfxTest() + +test.configure() +test.compile() +# this test succeeds if it runs successfully +test.run() + +test.mpi = True +test.configure() +test.compile() +# this test succeeds if it runs successfully +test.run()