Merge pull request #5 from boutproject/v2.1.0

v2.1.0

.gitignore

@@ -59,8 +59,9 @@ BOUT.settings
 *.bmp
 *.pdf
 
-storm3d/storm
-*/make.config
+make.config
+*storm
 data_*
 *.txt
 storm2d/storm2d
+storm3d/storm_version.hxx

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "shared/BoutEquation"]
+	path = shared/BoutEquation
+	url = [email protected]:johnomotani/BoutEquation

build_tools/save_git_version.py

@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+
+from pathlib import Path
+from string import ascii_uppercase
+
+try:
+    from git import Repo
+except ImportError:
+    # print a return code to stdout so we can capture and test it in the makefile
+    print(11)
+    raise
+
+scriptdir = Path(__file__).parent
+
+repo = Repo(scriptdir.parent)
+
+diff = repo.git.diff()
+
+# Find a string to use for the delimiter
+delimiter = "_STORM_DIFF_"
+for letter in ascii_uppercase:
+    if ")" + delimiter not in diff:
+        delimiter_success = True
+        break
+    delimiter = "_STORM_DIFF_" + letter + "_"
+if ")" + delimiter in diff:
+    # print a return code to stdout so we can capture and test it in the makefile
+    print(12)
+    raise ValueError(
+        "save_git_version.py failed to find a delimiter that is not in the git diff"
+    )
+
+with open("storm_version.hxx", "w") as f:
+    f.write("constexpr auto storm_git_hash = \"")
+    f.write(repo.git.describe(abbrev=40, dirty=True, always=True, tags=True))
+    f.write("\";\n")
+    f.write("constexpr auto storm_git_diff = R\"" + delimiter + "(")
+    f.write(diff)
+    f.write(")" + delimiter + "\";\n")
+
+# print a return code to stdout so we can capture and test it in the makefile
+print(0)

storm3d/D-vpar.cxx

@@ -5,8 +5,8 @@
 #include <derivs.hxx>
 #include <bout/constants.hxx>
 
-NeutralDVpar::NeutralDVpar(Solver *solver, Options &options)
-    : NeutralModel(options), my_output_monitor(this) {
+NeutralDVpar::NeutralDVpar(Solver *solver, Options &options, Datafile &dump)
+    : NeutralModel(options, dump), my_output_monitor(this) {
 
   mesh = bout::globals::mesh;
   coordinates_centre = mesh->getCoordinates(CELL_CENTRE);
@@ -165,8 +165,7 @@ void NeutralDVpar::update(const Field3D& n, const Field3D& n_stag, const Field3D
                           const Field3D& V, const Field3D& T, const Field3D& T_stag,
                           BoutReal time) {
   TRACE("NeutralDVpar::update");
-  ddt(lognn) = 0.;
-  ddt(nvn) = 0.;
+  rhs_counter++;
 
   Tn = Tn0/T0;
   Tn_stag = Tn0/T0;
@@ -278,24 +277,31 @@ void NeutralDVpar::update(const Field3D& n, const Field3D& n_stag, const Field3D
   }
   Field3D nn_vperp_GradPerpvn = nn_stag*vnperp_GradPerpvn;
 
-  ddt(lognn) += - vn_centre*myGrad_par(lognn_aligned, CELL_CENTRE)
-                - myDiv_par(vn_aligned, CELL_CENTRE)
-                + Dn*LaplPerpnn/nn
-                + GradPerpDn_GradPerpnn/nn
-                + S/nn
-                + Snn/nn;
-
-  ddt(nvn) += - vn*myGrad_par(nvn_aligned, CELL_YLOW)
-              - nvn*myDiv_par(vn_aligned, CELL_YLOW)
-              - Tn_stag*myGrad_par(nn_aligned, CELL_YLOW)
-              - vn*myInterp_to(div_vnperpnn_aligned, CELL_YLOW)
-              - nn_vperp_GradPerpvn
-              + myDiv_par(Dn_nn_Gradpar_vn_aligned, CELL_YLOW)
-              + Fi 
-              + Fe/mu
-              + Snvn;
-
-  if (munvn > 0.) ddt(nvn) += munvn*Delp2(nvn);
+  neutral_density_equation["parallel_advection"] =
+    - vn_centre*myGrad_par(lognn_aligned, CELL_CENTRE)
+    - myDiv_par(vn_aligned, CELL_CENTRE);
+  neutral_density_equation["perpendicular_diffusion"] =
+    Dn*LaplPerpnn/nn + GradPerpDn_GradPerpnn/nn;
+  neutral_density_equation["ionisation_recombination"] = S/nn;
+  neutral_density_equation["neutral_source"] = Snn/nn;
+
+  neutral_momentum_equation["parallel_advection"] =
+    - vn*myGrad_par(nvn_aligned, CELL_YLOW)
+    - nvn*myDiv_par(vn_aligned, CELL_YLOW);
+  neutral_momentum_equation["density_gradient"] =
+    - Tn_stag*myGrad_par(nn_aligned, CELL_YLOW);
+  neutral_momentum_equation["perpendicular_diffusion"] =
+    - vn*myInterp_to(div_vnperpnn_aligned, CELL_YLOW)
+    - nn_vperp_GradPerpvn;
+  neutral_momentum_equation["parallel_diffusion"] =
+    myDiv_par(Dn_nn_Gradpar_vn_aligned, CELL_YLOW);
+  neutral_momentum_equation["ion_friction"] = Fi;
+  neutral_momentum_equation["electron_friction"] = Fe/mu;
+  neutral_momentum_equation["neutral_momentum_source"] = Snvn;
+
+  if (munvn > 0.) {
+    neutral_momentum_equation["neutral_viscosity"] = munvn*Delp2(nvn);
+  }
 
   if (mesh->hasBndryLowerY()) {
     int j = mesh->ystart;
@@ -494,14 +500,18 @@ void NeutralDVpar::recycleFluxes(BoutReal time) {
 
   if (bndry_neutrals == 2 || bndry_neutrals == 3) {
     // Sources to recycle density and momentum in the first cell 
+    auto neutral_density_recycling =
+      neutral_density_equation["recycling_source"].localAccessor();
+    auto neutral_momentum_recycling =
+      neutral_momentum_equation["recycling_source"].localAccessor();
     j = mesh->ystart;
     for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) {
       i = r.ind;
       BoutReal dd = 1./( sqrt(coordinates_stag->g_22(i,j)) * coordinates_stag->dy(i,j) );
       for (int k = 0; k < mesh->LocalNz; k++) {
         BoutReal dnndt = recycled_lower(i,k)*dd;
-        ddt(lognn)(i,j,k) += dnndt/nn(i,j,k);
-        ddt(nvn)(i,j+1,k) += dnndt * vnth;
+        neutral_density_recycling(i,j,k) += dnndt/nn(i,j,k);
+        neutral_momentum_recycling(i,j+1,k) += dnndt * vnth;
       }
     }
 
@@ -511,8 +521,8 @@ void NeutralDVpar::recycleFluxes(BoutReal time) {
       BoutReal dd = 1./( sqrt(coordinates_stag->g_22(i,j+1)) * coordinates_stag->dy(i,j+1) );
       for (int k = 0; k < mesh->LocalNz; k++) {
         BoutReal dnndt = recycled_upper(i,k)*dd;
-        ddt(lognn)(i,j,k) += dnndt/nn(i,j,k);
-        ddt(nvn)(i,j,k)   -= dnndt * vnth;
+        neutral_density_recycling(i,j,k) += dnndt/nn(i,j,k);
+        neutral_momentum_recycling(i,j,k) -= dnndt * vnth;
       }
     }
   }

storm3d/D-vpar.hxx

@@ -7,12 +7,13 @@
 #define __NEUTRAL_DVPAR_H__
 
 #include "neutral-model.hxx"
+#include "../shared/BoutEquation/equation.hxx"
 #include <invert_laplace.hxx>
 #include <interpolation.hxx>
 
 class NeutralDVpar : public NeutralModel {
 public:
-  NeutralDVpar(Solver *solver, Options &options);
+  NeutralDVpar(Solver *solver, Options &options, Datafile &dump);
   ~NeutralDVpar() {}
 
   /// Update plasma quantities
@@ -27,6 +28,12 @@ public:
 
   void precon(BoutReal t, BoutReal gamma, BoutReal delta); 
 private:
+  int rhs_counter = 0;
+  Equation neutral_density_equation{lognn, "lognn", Options::root()["save_equations"],
+                                    dump, rhs_counter};
+  Equation neutral_momentum_equation{nvn, "nvn", Options::root()["save_equations"], dump,
+                                     rhs_counter};
+
   Field3D nn, nvn;
   Field3D Tn;
   Field3D lognn;