Merge pull request #8 from boutproject/v2.3.0-rc

STORM v2.3.0

.gitmodules

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

CMakeLists.txt

@@ -8,7 +8,7 @@ find_package(bout++
              REQUIRED)
 
 add_executable(storm
-               storm3d/storm.cxx storm3d/boundaries.cxx storm3d/initialise.cxx storm3d/monitors.cxx storm3d/operators.cxx storm3d/utilities.cxx shared/fast_output.cxx shared/BoutEquation/equation.cxx storm3d/neutral-rates.cxx storm3d/D-vpar.cxx storm3d/neutral-model.cxx)
+               storm3d/storm.cxx storm3d/boundaries.cxx storm3d/initialise.cxx storm3d/monitors.cxx storm3d/operators.cxx storm3d/utilities.cxx shared/BoutFastOutput/fast_output.cxx shared/BoutEquation/equation.cxx storm3d/neutral-rates.cxx storm3d/D-vpar.cxx storm3d/neutral-model.cxx)
 
 target_link_libraries(storm
                       PRIVATE bout++::bout++)

storm3d/D-vpar.cxx

@@ -6,7 +6,7 @@
 #include <bout/constants.hxx>
 
 NeutralDVpar::NeutralDVpar(Solver *solver, Options &options, Datafile &dump)
-    : NeutralModel(options, dump), my_output_monitor(this) {
+    : NeutralModel(options, dump), minmaxmean_monitor(*this) {
 
   mesh = bout::globals::mesh;
   coordinates_centre = mesh->getCoordinates(CELL_CENTRE);
@@ -157,7 +157,10 @@ NeutralDVpar::NeutralDVpar(Solver *solver, Options &options, Datafile &dump)
   }
   // Add this object to the global list, so it can be found by timestepmonitor_func
   neutral_dvpar_instances.push_back(this);
-  solver->addMonitor(&my_output_monitor);
+
+  if (monitor_minmaxmean) {
+    solver->addMonitor(&minmaxmean_monitor, MonitorPosition::BACK);
+  }
 }
 
 // 3D model, diffusive in X-Z, fluid in Y
@@ -337,23 +340,6 @@ int NeutralDVpar::timestepmonitor(BoutReal simtime) {
     updaterates = false;
   }
 
-  if( monitor_minmaxmean ) {
-    if(minmaxmean_timelast < 0.){
-      minmaxmean_timelast = simtime;
-    } else if (simtime - minmaxmean_timelast > 5.) {
-      minmaxmean_timelast = simtime;
-      printMinMaxMean();
-    }
-  }
-
-  return 0;
-}
-
-int NeutralDVpar::OutputMonitor::call(Solver* UNUSED(solver), BoutReal UNUSED(simtime),
-                                      int UNUSED(iter), int UNUSED(NOUT)) {
-  if (neutral_dvpar->monitor_minmaxmean) {
-    neutral_dvpar->printMinMaxMean();
-  }
   return 0;
 }
 

storm3d/D-vpar.hxx

@@ -75,7 +75,6 @@ private:
 
   bool updaterates;
   BoutReal monitor_timelast = -1.;
-  BoutReal minmaxmean_timelast = -1.;
   int timestepmonitor(BoutReal simtime);
 
   // Need global list of all NeutralDVpar instances, to iterate over in the
@@ -102,15 +101,32 @@ private:
     return fromFieldAligned( Div_par(var_aligned, outloc) , "RGN_NOBNDRY");
   }
 
-  class OutputMonitor : public Monitor {
+  class MinMaxMeanMonitor : public Monitor {
   public:
-    OutputMonitor(NeutralDVpar* ndvp_in, BoutReal timestep=-1)
-      : Monitor(timestep), neutral_dvpar(ndvp_in) {};
-    int call(Solver *solver, BoutReal simtime, int iter, int NOUT) override;
+    MinMaxMeanMonitor(NeutralDVpar& ndvp_in)
+        : neutral_dvpar(ndvp_in) {
+
+      BoutReal timestep = Options::root()["TIMESTEP"]
+                              .doc("Output time step size").withDefault(1.0);
+
+      BoutReal minmaxmean_timestep = 5.0;
+      if (timestep > minmaxmean_timestep) {
+        // Find integer fraction of timestep that's close to minmaxmean_timestep
+        int n = std::round(timestep / minmaxmean_timestep);
+
+        // Set MinMaxMeanMonitor 'timestep' to something close to minmaxmean_timestep
+        setTimestep(timestep / n);
+      }
+    };
+    int call(Solver* UNUSED(solver), BoutReal UNUSED(simtime), int UNUSED(iter),
+             int UNUSED(NOUT)) override {
+      neutral_dvpar.printMinMaxMean();
+      return 0;
+    }
   private:
-    NeutralDVpar* neutral_dvpar;
+    NeutralDVpar& neutral_dvpar;
   };
-  OutputMonitor my_output_monitor;
+  MinMaxMeanMonitor minmaxmean_monitor;
 
   // Print minimum, maximum and mean of evolving variables. Used when
   // monitor_minmaxmean=true

storm3d/makefile

@@ -27,8 +27,11 @@ endif
 ifeq ($(wildcard ../shared/BoutEquation/equation.hxx),)
   $(error "BoutEquation header missing - try running `git submodule update --init --recursive`")
 endif
+ifeq ($(wildcard ../shared/BoutFastOutput/fast_output.hxx),)
+  $(error "FastOutput header missing - try running `git submodule update --init --recursive`")
+endif
 
-SOURCEC = storm.cxx boundaries.cxx initialise.cxx monitors.cxx operators.cxx utilities.cxx ../shared/fast_output.cxx ../shared/BoutEquation/equation.cxx neutral-rates.cxx D-vpar.cxx neutral-model.cxx
+SOURCEC = storm.cxx boundaries.cxx initialise.cxx monitors.cxx operators.cxx utilities.cxx ../shared/BoutFastOutput/fast_output.cxx ../shared/BoutEquation/equation.cxx neutral-rates.cxx D-vpar.cxx neutral-model.cxx
 TARGET = storm
 
 save_git_version_status := $(shell ../build_tools/save_git_version.py)

storm3d/monitors.cxx

@@ -56,16 +56,15 @@ void STORM::printMinMaxMean() {
     output.write("min(psi) = %e, max(psi) = %e, mean(psi) = %e\n",
       min(psi,true,"RGN_NOBNDRY"),max(psi,true,"RGN_NOBNDRY"),mean(psi,true,"RGN_NOBNDRY"));
   }
+
+  output.write("\n");
 }
 
 int STORM::outputMonitor(BoutReal UNUSED(simtime), int UNUSED(iteration), int UNUSED(nout)) {
   if (hydrodynamic) {
     phisolver_1d();
     phi.applyBoundary();
   }
-  if (monitor_minmaxmean) {
-    printMinMaxMean();
-  }
 
   return 0;
 }

storm3d/storm.cxx

@@ -404,6 +404,7 @@ int STORM::init(bool restarting) {
   n_stag.setLocation(CELL_YLOW) ;
   S_stag.setLocation(CELL_YLOW) ;
   phi_stag.setLocation(CELL_YLOW) ;  
+  nu_parallel0.setLocation(CELL_YLOW) ;
   nu_parallel.setLocation(CELL_YLOW) ;
   qpar_aligned.setLocation(CELL_YLOW);
   logT_stag.setLocation(CELL_YLOW);
@@ -818,6 +819,12 @@ int STORM::init(bool restarting) {
     nu_parallel0 = 0.1;
   }
 
+  // Add monitor for printing min, max and mean of evolving variables and their time
+  // derivatives
+  if (monitor_minmaxmean) {
+    solver->addMonitor(&minmaxmean_monitor, MonitorPosition::BACK);
+  }
+
   /////////////////////////////////////////////////////////
   // Neutral models
 
@@ -1308,8 +1315,6 @@ int STORM::rhs(BoutReal time) {
     if (diff_perp_V > 0.) {
       electron_momentum_equation["perpendicular_diffusion"] = diff_perp_V*Delp2(V);
     }
-
-    set_lower_ddt_zero(chiV) ;
   }
 
   if (!isothermal){
@@ -1365,7 +1370,7 @@ int STORM::rhs(BoutReal time) {
     }
 
     if (S_in_peq) {
-      electron_pressure_equation["S*V^2/3/mu/p"] = S*SQ(V_centre)/(3.*mu*p);
+      electron_pressure_equation["S*V^2_3_mu_p"] = S*SQ(V_centre)/(3.*mu*p);
     }
 
     //Curvature terms for electron temperature