V2.0.03 (#213)

* std::filesystem compatibility (#210)

* Improve and debug slices (#211)

* CLN: cleanup redefinition of idfx::randm in planet tests (#212)

* Improve lookup table (#214)

* fix slice VTK outputs when they are at a boundary (#215)

* use SlopeLimiter class's PPM implementation (#218)

* fix an overflow in dump restart routines (#219)

* fix a bug identified in gcc 9.3.0  (#220)

* add documentation for global IdefixArrays (#221)

---------

Co-authored-by: Clément Robert <cr52@protonmail.com>

Author: glesur

CHANGELOG.md

@@ -4,6 +4,16 @@ 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.0.3] 2023-11-24
+### Changed
+- fix and improve slice VTK outputs (#211)
+- improve std::filesystem compatibility with old stdlibc++ (#210)
+- cleanup planet test (#212)
+- allow the use of IdefixHostArrays as inputs for the LookupTable class (#214)
+- use Felker & Stone (2018) PPM scheme when using HIGH_ORDER_FARGO in place of Colella & Woodward (1984) (#218).
+- fix a bug that could result in segfaults when reading a restart dump, in particular in 1D problems with MPI (#219).
+- fix a bug preventing compilation on gcc 9.3.0 (#220)
+
 ## [2.0.2] 2023-11-6
 ### Changed
 - fixed compilation issues with gcc 8.x (#205)

CMakeLists.txt

@@ -4,9 +4,9 @@ set (CMAKE_CXX_STANDARD 17)
 
 set(Idefix_VERSION_MAJOR 2)
 set(Idefix_VERSION_MINOR 0)
-set(Idefix_VERSION_PATCH 01)
+set(Idefix_VERSION_PATCH 03)
 
-project (idefix VERSION 1.0.0)
+project (idefix VERSION 2.0.03)
 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)

doc/source/faq.rst

@@ -39,6 +39,13 @@ The compilation stops while compiling Kokkos with ``/usr/include/stdlib.h(58): e
 When using the Intel compiler on a Mac Intel, I get a linking error involving the ``SharedAllocationRecordIvvE18t_tracking_enabledE`` symbol.
   This is a known bug of the Intel Mac compiler with Kokkos. Apparently Intel has decided not to fix it. Check the issue on the `Kokkos git page <https://github.com/kokkos/kokkos/issues/1959>`_.
 
+I get an error during *Idefix* link that says there are undefined symbols to std::filesystem
+  This is a known bug/limitation of the stdlibc++ provided with gcc8 that does not include C++17 filesystem extensions.
+  While *Idefix* auto-detects gcc8 when it is used as the main compiler, it still misses the cases when another compiler
+  (like Clang or Intel) is used with gcc8 as a backend.
+  You should try to clear up CMakeCache.txt and explicitely add the required link library when calling cmake as in
+  ``LDFLAGS=-lstdc++fs cmake $IDEFIX_DIR ...```
+
 Execution
 ---------
 

doc/source/reference/setup.cpp.rst

@@ -395,3 +395,80 @@ User-defined analysis
 
 User-defined analysis and outputs can be coded in the ``setup.cpp`` file. Follow the
 guidelines in :ref:`output`.
+
+I need a global IdefixArray for my Setup
+-----------------------------------------
+
+There are situation where you will need one or several global IdefixArrays that can be accessed from different
+functions, e.g. the ``Initflow`` method and the user-defined boundary conditions.
+
+It is important to understand that IdefixArrays (equivalent to ``Kokkos::view`` that are references to memory chunks)
+are automatically dealocated when all of the IdefixArrays refereing to that memory chunk have been deleted. This deletion happens either
+implicitly (by a closing scope) in which case the objects contained in the scope are all deleted automatically,
+or explicitly (through a new/delete pair).
+
+If you define an IdefixArray in the global scope, it is deleted when the program terminates. Hence deallocation should happen then.
+Except that, according to `Kokkos documentation <https://kokkos.github.io/kokkos-core-wiki/ProgrammingGuide/Initialization.html#finalization>`_,  we
+need to call ``Kokkos::finalize`` before the program terminates and this ``finalize`` should be done once all of
+the Kokkos objects have been deleted (including IdefixArray). While *Idefix* makes sure that all of its objects (including the user's ``Setup``) are being deleted before calling ``finalize``,
+a simple IdefixArray in the global scope will not be explicitely deleted, and will typically lead to the following error:
+
+.. code-block:: bash
+
+  terminate called after throwing an instance of 'std::runtime_error'
+  what():  Kokkos allocation "MyAwesomeArray" is being deallocated after Kokkos::finalize was called
+
+The way to avoid this is to explicitely delete the object when you don't need it anymore. The cleanest way to do this for a setup is to define a "container" class,
+containing all of the arrays you will need in the global scope, and just have a global pointer to an instance of this class, that you eventually delete
+(and which deletes all of the arrays it contains automatically). More explicitely:
+
+#. start with a declaration of a class container (that we name MyGlobalClass in this example) and a global pointer to a class instance (note that you can put as many arrays as you want in the class)
+
+    .. code-block:: c++
+
+      // Class declaration
+      class MyGlobalClass {
+      public:
+        // Class constructor
+        MyGlobalClass(DataBlock &data) {
+        //allocate some memory for the array the class contains
+          this->array1 = IdefixArray3D<real>("MyAwesomeArray",data.np_tot[KDIR], data.np_tot[JDIR], data.np_tot[IDIR]);
+        }
+
+        // array1, member of the class
+        IdefixArray3D<real> array1;
+      };
+
+      // A global class instance named "myGlobals"
+      MyGlobalClass *myGlobals;
+
+#. initialise your global object in the Setup constructor (this will aumatically allocate the array it contains thanks to the class constructor we have defined):
+
+    .. code-block:: c++
+
+      Setup::Setup(....) {
+        ...
+        myGlobals = new MyGlobalClass(data);
+        ...
+      }
+
+#. to avoid the error message above, don't forget to delete the object on exit in the Setup destructor
+
+    .. code-block:: c++
+
+      Setup::~Setup(....) {
+        ...
+        delete myGlobals;
+        ...
+      }
+
+#. and finally, use your array when you need it:
+
+    .. code-block:: c++
+
+      MyXXXXFunction(....) {
+        // Shallow copy the global array
+        IdefixArray3D<real> array = myGlobals->array1;
+        // Do stuff
+        ....
+      }

src/dataBlock/fargo.hpp

@@ -14,8 +14,10 @@
   #include "mpi.hpp"
 #endif
 
-// Forward class hydro declaration
 #include "physics.hpp"
+#include "slopeLimiter.hpp"
+
+// Forward class hydro declaration
 template <typename Phys> class Fluid;
 using Hydro = Fluid<DefaultPhysics>;
 class DataBlock;
@@ -85,15 +87,6 @@ class Fargo {
 #endif
 
 #ifdef HIGH_ORDER_FARGO
-KOKKOS_FORCEINLINE_FUNCTION real PPMLim(real dvp, real dvm) {
-  if(dvp*dvm >0.0) {
-    real dqc = 0.5*(dvp+dvm);
-    real d2q = 2.0*( fabs(dvp) < fabs(dvm) ? dvp : dvm);
-    return( fabs(d2q) < fabs(dqc) ? d2q : dqc);
-  }
-  return(ZERO_F);
-}
-
 KOKKOS_INLINE_FUNCTION real FargoFlux(const IdefixArray4D<real> &Vin, int n, int k, int j, int i,
                                       int so, int ds, int sbeg, real eps,
                                       bool haveDomainDecomposition) {
@@ -108,51 +101,36 @@ KOKKOS_INLINE_FUNCTION real FargoFlux(const IdefixArray4D<real> &Vin, int n, int
   int som2 = som1-1;
   if(!haveDomainDecomposition && (som2-sbeg< 0 )) som2 = som2+ds;
 
-  real dqm2,dqm1,dqp1,dqp2, q0,qm1, qp1;
+  real q0,qm1, qp1, qm2, qp2;
   #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
     q0 = Vin(n,k,so,i);
     qm1 = Vin(n,k,som1,i);
     qp1 = Vin(n,k,sop1,i);
-    dqm2 = qm1 - Vin(n,k,som2,i);
-    dqm1 = q0 - qm1;
-    dqp1 = qp1 - q0;
-    dqp2 = Vin(n,k,sop2,i) - qp1;
+    qm2 = Vin(n,k,som2,i);
+    qp2 = Vin(n,k,sop2,i);
   #elif GEOMETRY == SPHERICAL
     q0 = Vin(n,so,j,i);
     qm1 = Vin(n,som1,j,i);
     qp1 = Vin(n,sop1,j,i);
-    dqm2 = qm1 - Vin(n,som2,j,i);
-    dqm1 = q0 - qm1;
-    dqp1 = qp1 - q0;
-    dqp2 = Vin(n,sop2,j,i) - qp1;
+    qm2 = Vin(n,som2,j,i);
+    qp2 = Vin(n,sop2,j,i);
   #endif
-    // slope limited values around the reference point
-    real dqlm = PPMLim(dqm1,dqm2);
-    real dql0 = PPMLim(dqp1,dqm1);
-    real dqlp = PPMLim(dqp2,dqp1);
-
-    real dqp = 0.5 * dqp1 - (dqlp - dql0) / 6.0;
-    real dqm = -0.5 * dqm1 - (dql0 - dqlm) / 6.0;
-
-    if(dqp*dqm>0.0) {
-       dqp = dqm = 0.0;
-    } else {
-      if(FABS(dqp) >= 2.0*FABS(dqm)) dqp = -2.0*dqm;
-      if(FABS(dqm) >= 2.0*FABS(dqp)) dqm = -2.0*dqp;
-    }
 
-    real qp = q0 + dqp;
-    real qm = q0 + dqm;
+  real qp, qm;
+  SlopeLimiter<>::getPPMStates(qm2,qm1,q0,qp1,qp2, qm, qp);
 
-    real dqc = dqp - dqm;
-    real d2q = dqp + dqm;
+  real dqp = qp-q0;
+  real dqm = qm-q0;
 
-    real F;
-    if(eps > 0.0) {
-      F = eps*(qp - 0.5*eps*(dqc + d2q*(3.0 - 2.0*eps)));
-    } else {
-      F = eps*(qm - 0.5*eps*(dqc - d2q*(3.0 + 2.0*eps)));
-    }
+  real dqc = dqp - dqm;
+  real d2q = dqp + dqm;
+
+  real F;
+  if(eps > 0.0) {
+    F = eps*(qp - 0.5*eps*(dqc + d2q*(3.0 - 2.0*eps)));
+  } else {
+    F = eps*(qm - 0.5*eps*(dqc - d2q*(3.0 + 2.0*eps)));
+  }
 
   return(F);
 }

src/fluid/braginskii/bragThermalDiffusion.hpp

@@ -73,12 +73,12 @@ class BragThermalDiffusion {
 
 template <typename Phys>
 BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid<Phys> *hydroin):
-                            Vc{hydroin->Vc},
-                            Vs{hydroin->Vs},
-                            dMax{hydroin->dMax},
-                            eos{hydroin->eos.get()},
-                            data{hydroin->data},
-                            status{hydroin->bragThermalDiffusionStatus} {
+                            Vc(hydroin->Vc),
+                            Vs(hydroin->Vs),
+                            dMax(hydroin->dMax),
+                            eos(hydroin->eos.get()),
+                            data(hydroin->data),
+                            status(hydroin->bragThermalDiffusionStatus) {
   idfx::pushRegion("BragThermalDiffusion::BragThermalDiffusion");
 
   if(input.CheckEntry("Hydro","bragTDiffusion")>=0) {

src/fluid/braginskii/bragViscosity.hpp

@@ -67,10 +67,10 @@ class BragViscosity {
 
 template<typename Phys>
 BragViscosity::BragViscosity(Input &input, Grid &grid, Fluid<Phys> *hydroin):
-                      Vc{hydroin->Vc},
-                      Vs{hydroin->Vs},
-                      dMax{hydroin->dMax},
-                      status{hydroin->bragViscosityStatus} {
+                      Vc(hydroin->Vc),
+                      Vs(hydroin->Vs),
+                      dMax(hydroin->dMax),
+                      status(hydroin->bragViscosityStatus) {
   idfx::pushRegion("BragViscosity::BragViscosity");
   // Save the parent hydro object
   this->data = hydroin->data;

src/fluid/thermalDiffusion.hpp

@@ -64,11 +64,11 @@ class ThermalDiffusion {
 
 template <typename Phys>
 ThermalDiffusion::ThermalDiffusion(Input &input, Grid &grid, Fluid<Phys> *hydroin):
-                            Vc{hydroin->Vc},
-                            dMax{hydroin->dMax},
-                            eos{hydroin->eos.get()},
-                            data{hydroin->data},
-                            status{hydroin->thermalDiffusionStatus} {
+                            Vc(hydroin->Vc),
+                            dMax(hydroin->dMax),
+                            eos(hydroin->eos.get()),
+                            data(hydroin->data),
+                            status(hydroin->thermalDiffusionStatus) {
   idfx::pushRegion("ThermalDiffusion::ThermalDiffusion");
 
   if(status.status == Constant) {

src/fluid/viscosity.hpp

@@ -66,9 +66,9 @@ class Viscosity {
 
 template<typename Phys>
 Viscosity::Viscosity(Input &input, Grid &grid, Fluid<Phys> *hydroin):
-                      Vc{hydroin->Vc},
-                      dMax{hydroin->dMax},
-                      status{hydroin->viscosityStatus} {
+                      Vc(hydroin->Vc),
+                      dMax(hydroin->dMax),
+                      status(hydroin->viscosityStatus) {
   idfx::pushRegion("Viscosity::Viscosity");
   // Save the parent hydro object
   this->data = hydroin->data;