Imex time stepping (#775)

* imex2 and vplb example

Author: mkstoyanov

CMakeLists.txt

@@ -504,13 +504,14 @@ if (${ASGARD_USE_PCH})
     target_precompile_headers (asgard_exe REUSE_FROM libasgard)
 endif ()
 
-set (_asgard_pdes continuity diffusion elliptic sinwav two_stream
+set (_asgard_pdes continuity diffusion elliptic sinwav two_stream vplb
                   mass_internal varcoeff_internal bound_internal relaxation_internal)
 
 foreach (_asgexe ${_asgard_pdes})
     add_executable (${_asgexe} "${CMAKE_CURRENT_SOURCE_DIR}/src/pde/${_asgexe}.cpp")
     target_link_libraries (${_asgexe} PUBLIC libasgard)
-    install(TARGETS ${_asgexe} EXPORT  "asgard-export" RUNTIME DESTINATION "share/asgard/pde")
+    install(TARGETS ${_asgexe} RUNTIME DESTINATION "share/asgard/pde")
+    install(FILES "${CMAKE_CURRENT_SOURCE_DIR}/src/pde/${_asgexe}.cpp" DESTINATION "share/asgard/examples")
 endforeach()
 
 #-------------------------------------------------------------------------------

doxygen/CMakeLists.txt

@@ -53,13 +53,14 @@ doxygen_add_docs(asgard_doxygen
                  src/asgard_reconstruct.hpp
                  src/asgard_dimension.hpp
                  src/pde/asgard_pde_base.hpp
-                 examples/continuity_2d.cpp
                  examples/inputs_1d.cpp
+                 examples/continuity_2d.cpp
                  src/pde/continuity.cpp
                  src/pde/diffusion.cpp
                  src/pde/elliptic.cpp
                  src/pde/sinwav.cpp
                  src/pde/two_stream.cpp
+                 src/pde/vplb.cpp
                  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/../
                  COMMENT "Building the ${PROJECT_NAME} documentation")
 

examples/CMakeLists.txt

@@ -20,3 +20,9 @@ if (asgard_PYTHON_FOUND)
 
     configure_file("${CMAKE_CURRENT_SOURCE_DIR}/inputs_1d.py"   "${CMAKE_CURRENT_BINARY_DIR}/inputs_1d.py" COPYONLY)
 endif()
+
+# additional examples following the same pattern
+foreach (_example continuity diffusion elliptic sinwav two_stream vplb)
+    add_executable(${_example}  ${_example}.cpp)
+    target_link_libraries(${_example}  asgard::asgard)
+endforeach()

examples/continuity_2d.cpp

@@ -12,9 +12,9 @@
 
 /*!
  * \ingroup asgard_examples
- * \addtogroup asgard_examples_continuity_2d Example 1, 2D continuity equation
+ * \addtogroup asgard_examples_continuity_2d Example: 2D continuity equation
  *
- * \par Example 1
+ * \par 2D continuity equation
  * Solves the continuity partial differential equation
  * \f[ \frac{\partial}{\partial t} f + \nabla \cdot f = s \f]
  * where both \b f and \b s are defined over the two dimensional domain

examples/inputs_1d.cpp

@@ -12,9 +12,9 @@
 
 /*!
  * \ingroup asgard_examples
- * \addtogroup asgard_examples_input_1d Example 2, Simple 1D equation
+ * \addtogroup asgard_examples_input_1d Example: Simple 1D equation
  *
- * \par Example 2
+ * \par Simple 1D equation
  * Solves the continuity partial differential equation
  * \f[ \frac{\partial}{\partial t} f + \frac{\partial}{\partial x} f = s \f]
  * where both \b f and \b s are defined over domain

python/asgard.py

@@ -54,6 +54,10 @@ def __del__(self):
 
     def __init__(self, filename, verbose = False):
         self.verbose = verbose
+        if filename == "::aux-filed":
+            self.recsol = None
+            return
+
         if self.verbose:
             print(' -- reading from: %s' % filename)
 
@@ -98,6 +102,9 @@ def __init__(self, filename, verbose = False):
 
                 self.num_dimensions = fdata['num_dims'][()]
 
+                self.num_position = fdata['num_pos'][()]
+                self.num_velocity = fdata['num_vel'][()]
+
                 self.cells = fdata['grid_indexes'][()]
                 self.time  = fdata['dtime_time'][()] # numeric time
 
@@ -114,6 +121,15 @@ def __init__(self, filename, verbose = False):
                     self.dimension_min[i] = drange[2 * i]
                     self.dimension_max[i] = drange[2 * i + 1]
 
+                num_aux = fdata['num_aux_fields'][()]
+                self.aux_fields = [None for i in range(num_aux)]
+                for i in range(num_aux):
+                    self.aux_fields[i] = {
+                        'name' : fdata[f"aux_field_{i}_name"][()].decode("utf-8"),
+                        'data' : fdata[f"aux_field_{i}_data"][()],
+                        'grid' : fdata[f"aux_field_{i}_grid"][()]
+                        }
+
         # for plotting purposes, say aways from the domain edges
         # rounding error at the edge may skew the plots
         self.eps = 1.E-6 * np.min(self.dimension_max - self.dimension_min)
@@ -158,6 +174,66 @@ def __init__(self, filename, verbose = False):
                                                         np.ctypeslib.as_ctypes(self.dimension_min.reshape(-1,)),
                                                         np.ctypeslib.as_ctypes(self.dimension_max.reshape(-1,)))
 
+    def get_aux_field(self, auxid):
+        assert isinstance(auxid, int) or isinstance(auxid, string), "auxid must be in int or a string"
+        if isinstance(auxid, int):
+            assert 0 <= auxid and auxid < len(self.aux_fields), f"the auxid {auxid} must point to a valid entry in the list with size {len(self.aux_fields)}"
+
+            idnum = auxid
+        else: # must be a string due to the assertion on top
+            idnum = -1
+            for i in range(len(self.aux_fields)):
+                if self.aux_fields[i].name == auxid:
+                    idnum = i
+                    break
+            assert idnum != -1, f"the auxid '{auxid}' is not in the list of aux-fields"
+
+        aux = pde_snapshot("::aux-filed", self.verbose)
+
+        aux.title    = self.aux_fields[idnum]['name']
+        aux.subtitle = "aux-field"
+        aux.degree   = self.degree
+        aux.state    = self.aux_fields[idnum]['data']
+        aux.cells    = self.aux_fields[idnum]['grid']
+
+        aux.default_view = self.default_view
+
+        aux.num_dimensions = self.num_dimensions
+        aux.num_position   = self.num_position
+        aux.num_velocity   = self.num_velocity
+        aux.num_cells      = aux.cells.shape[0] / aux.num_dimensions
+
+        aux.time  = self.time
+        aux.time  = self.time
+
+        aux.dimension_names = self.dimension_names
+        aux.dimension_min   = self.dimension_min
+        aux.dimension_max   = self.dimension_max
+
+        aux.aux_fields = []
+
+        aux.eps = self.eps
+
+        if self.state.dtype == np.float64:
+            aux.double_precision = True
+
+            aux.recsol = libasgard.asgard_make_dreconstruct_solution_v2(
+                self.num_dimensions, self.num_cells, np.ctypeslib.as_ctypes(aux.cells.reshape(-1,)),
+                self.degree, np.ctypeslib.as_ctypes(aux.state.reshape(-1,)))
+
+        else:
+            aux.double_precision = False
+
+            aux.recsol = libasgard.asgard_make_freconstruct_solution_v2(
+                self.num_dimensions, self.num_cells, np.ctypeslib.as_ctypes(aux.cells.reshape(-1,)),
+                self.degree, np.ctypeslib.as_ctypes(aux.state.reshape(-1,)))
+
+        libasgard.asgard_reconstruct_solution_setbounds(aux.recsol,
+                                                        np.ctypeslib.as_ctypes(self.dimension_min.reshape(-1,)),
+                                                        np.ctypeslib.as_ctypes(self.dimension_max.reshape(-1,)))
+
+        return aux
+
     def plot_data1d(self, dims, num_points = 32):
         '''
         Generates two 1d arrays (vals, x), so that vals is the values
@@ -301,8 +377,11 @@ def cell_centers(self):
     def __str__(self):
         s = "title: %s\n" % self.title
         if self.subtitle != "":
-            s += "   sub: %s\n" % self.subtitle
-        s += "  num-dimensions: %d\n" % self.num_dimensions
+            s += "        %s\n" % self.subtitle
+        if self.num_position == 0 and self.num_velocity == 0:
+            s += "  num-dimensions: %d\n" % self.num_dimensions
+        else: # have position/velocity dimensions
+            s += f"  num-dimensions: {self.num_dimensions}  ({self.num_position}x{self.num_velocity}v)\n"
         s += "  degree:         %d\n" % self.degree
         s += "  num-indexes:    %d\n" % self.num_cells
         s += "  state size:     %d\n" % self.state.size
@@ -374,6 +453,7 @@ def __str__(self):
         # we are plotting, consider extra options
         plotview = None
         savefig  = None
+        auxfield = None
         addgrid  = False
         if len(sys.argv) > 2:
             i = 2
@@ -385,13 +465,21 @@ def __str__(self):
                 elif sys.argv[i] == "-fig":
                     savefig = sys.argv[i + 1] if i + 1 < n else None
                     i += 2
+                elif sys.argv[i] == "-aux":
+                    auxfield = sys.argv[i + 1] if i + 1 < n else None
+                    i += 2
+                    assert auxfield is not None, "-aux requires an filed number"
+                    auxfield = int(auxfield)
                 elif sys.argv[i] == "-grid":
                     addgrid = True
                     i += 1
                 else:
                     savefig = sys.argv[i]
                     i += 1
 
+        if auxfield is not None:
+            shot = shot.get_aux_field(auxfield)
+
         asgplot.title(shot.title, fontsize = 'large')
 
         if shot.num_dimensions == 1:

src/asgard_coefficients_mats.hpp

@@ -379,7 +379,7 @@ void gen_diag_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
 //! moment over moment zero
 template<typename P, int multsign, pterm_dependence dep>
 void gen_diag_mom_cases(
-    legendre_basis<P> const &basis, P xleft, P xright, int level, int mindex,
+    legendre_basis<P> const &basis, int level, int mindex,
     std::vector<P> const &moms, block_diag_matrix<P> &coefficients)
 {
   static_assert(multsign == 1 or multsign == -1);
@@ -392,9 +392,6 @@ void gen_diag_mom_cases(
 
   coefficients.resize_and_zero(pdof * pdof, num_cells);
 
-  P const dx = (xright - xleft) / num_cells;
-  P const vscale = P{1} / dx; // volume scale
-
   size_t const wsize = [&]() -> size_t {
     if constexpr (dep == pterm_dependence::moment_divided_by_density)
       return num_quad * pdof + 2 * num_quad;

src/asgard_dimension.hpp

@@ -243,14 +243,11 @@ class pde_domain
 
 private:
   void check_init() {
-    if (num_dims_ < 1)
-      throw std::runtime_error("pde_domain created with zero or negative dimensions");
-    if (num_dims_ > max_num_dimensions)
-      throw std::runtime_error("pde_domain created with too many dimensions, max is 6D");
-    if (num_pos_ < 0)
-      throw std::runtime_error("pde_domain created with negative position dimensions");
-    if (num_vel_ < 0)
-      throw std::runtime_error("pde_domain created with negative velocity dimensions");
+    rassert(num_dims_ >= 1, "pde_domain created with zero or negative dimensions");
+    rassert(num_dims_ <= max_num_dimensions,
+            "pde_domain created with too many dimensions, max is 6D");
+    rassert(num_pos_ >= 0, "pde_domain created with negative position dimensions");
+    rassert(num_vel_ >= 0, "pde_domain created with negative velocity dimensions");
 
     if (num_pos_ == 0 and num_vel_ == 0) {
       for (int d : iindexof(num_dims_))
@@ -426,4 +423,32 @@ class separable_func
   scalar_func<P> time_func_;
 };
 
+/*!
+ * \ingroup asgard_discretization
+ * \brief Extra data-entry for plotting and post-processing
+ *
+ * In plotting and post-processing, it is sometime desirable to store
+ * additional data that sits on the sparse grid mesh, e.g.,
+ * deviation from a nominal state or initial condition.
+ * Since the data is defined on a sparse grid, it has to be accessed with
+ * the asgard::reconstruct_solution class (e.g., via python), but the data
+ * has to be saved/loaded in the asgard::discretization_manager
+ */
+template<typename P>
+struct aux_field_entry {
+  //! default constructor, creates and empty entry
+  aux_field_entry() = default;
+  //! constructor, set the name and data
+  aux_field_entry(std::string nm, std::vector<P> dat)
+      : name(std::move(nm)), data(std::move(dat))
+  {}
+  //! reference name for the field, should be unique
+  std::string name;
+  //! vector data
+  std::vector<P> data;
+  //! multi-indexes
+  std::vector<int> grid;
+};
+
+
 } // namespace asgard