reinstated dtFactor and added Po read/write

src/calorically_perfect.hpp

@@ -102,7 +102,8 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   double static_rho_;
 
   /// pressure-related, closed-system thermo pressure changes
-  double ambient_pressure_, thermo_pressure_, system_mass_;
+  double ambient_pressure_, system_mass_;
+  // double thermo_pressure_;
   double dtP_;
 
   // Initial temperature value (if constant IC)
@@ -241,6 +242,10 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   /// Return a pointer to the current total thermal diffusivity ParGridFunction.
   ParGridFunction *GetCurrentThermalDiv() { return &Qt_gf_; }
 
+  /// Return thermodynamic pressure for restarts
+  // double GetCurrentThermoPressure() { return thermo_pressure_; }
+  // double SetCurrentThermoPressure(double Po) { thermo_pressure_ = Po;}
+
   /// Rotate entries in the time step and solution history arrays.
   void UpdateTimestepHistory(double dt);
 

src/loMach.cpp

@@ -217,6 +217,9 @@ void LoMachSolver::initialize() {
   // If restarting, read restart files
   if (loMach_opts_.io_opts_.enable_restart_) {
     restart_files_hdf5("read");
+    if (thermoPressure_ > 0.0) {
+      thermo_->SetThermoPressure(thermoPressure_);
+    }
   }
 
   // Exchange interface information
@@ -347,6 +350,8 @@ void LoMachSolver::solveStep() {
 
   // restart files
   if (iter % loMach_opts_.output_frequency_ == 0 && iter != 0) {
+    thermoPressure_ = thermo_->GetThermoPressure();
+
     // Write restart file!
     restart_files_hdf5("write");
 
@@ -405,15 +410,15 @@ void LoMachSolver::solve() {
 
 void LoMachSolver::updateTimestep() {
   // minimum timestep to not waste comp time
-  double dtMin = 1.0e-9;
+  double dtMin = loMach_opts_.ts_opts_.minimum_dt_;
+  double dtMax = loMach_opts_.ts_opts_.maximum_dt_;
+  double dtFactor = loMach_opts_.ts_opts_.factor_dt_;
 
   double Umax_lcl = 1.0e-12;
   double max_speed = Umax_lcl;
   double Umag;
   // int Sdof = sfes_->GetNDofs();
   // int Sdof = (turbModel_->getGridScale())->Size();
-  // TODO(trevilo): Let user set dtFactor
-  double dtFactor = 1.0;
   auto dataU = flow_->getCurrentVelocity()->HostRead();
 
   // come in divided by order
@@ -471,7 +476,6 @@ void LoMachSolver::updateTimestep() {
   // double dtInst = max(dtInst_conv, dtInst_visc);
 
   double dtInst = dtInst_conv;
-
   double &dt = temporal_coeff_.dt;
   if (dtInst > dt) {
     dt = dt * (1.0 + dtFactor);
@@ -480,6 +484,10 @@ void LoMachSolver::updateTimestep() {
     dt = dtInst;
   }
 
+  if (dt > dtMax) {
+    dt = dtMax;
+  }
+
   if (dt < dtMin) {
     if (rank0_) {
       std::cout << " => Timestep running away!!!" << endl;
@@ -587,6 +595,12 @@ void LoMachSolver::setTimestep() {
     MPI_Allreduce(&Umax_lcl, &max_speed, 1, MPI_DOUBLE, MPI_MAX, pmesh_->GetComm());
     double dtInst = CFL * hmin / (max_speed * (double)order);
     temporal_coeff_.dt = dtInst;
+
+    const double dt_initial = loMach_opts_.ts_opts_.initial_dt_;
+    if ((dt_initial < dtInst) && !(loMach_opts_.io_opts_.enable_restart_)) {
+      temporal_coeff_.dt = dt_initial;
+    }
+
     std::cout << "dt from setTimestep: " << temporal_coeff_.dt << " max_speed: " << max_speed << endl;
   }
 }

src/loMach.hpp

@@ -145,6 +145,9 @@ class LoMachSolver : public TPS::Solver {
   // min/max element size
   double hmin, hmax;
 
+  // for restarts
+  double thermoPressure_ = -1.0;
+
   // domain extent
   double xmin_, ymin_, zmin_;
   double xmax_, ymax_, zmax_;

src/loMachIO.cpp

@@ -58,6 +58,8 @@ void LoMachSolver::write_restart_files_hdf5(hid_t file, bool serialized_write) {
     h5_save_attribute(file, "order", order);
     // spatial dimension
     h5_save_attribute(file, "dimension", dim_);
+    // thermodynamic pressure
+    h5_save_attribute(file, "Po", thermoPressure_);
 
     // code revision
 #ifdef BUILD_VERSION
@@ -111,6 +113,7 @@ void LoMachSolver::read_restart_files_hdf5(hid_t file, bool serialized_read) {
     h5_read_attribute(file, "time", temporal_coeff_.time);
     h5_read_attribute(file, "dt", temporal_coeff_.dt);
     h5_read_attribute(file, "order", read_order);
+    h5_read_attribute(file, "Po", thermoPressure_);
   }
 
   if (serialized_read) {

src/loMach_options.cpp

@@ -101,6 +101,9 @@ void LoMachTemporalOptions::read(TPS::Tps *tps, std::string prefix) {
 
   tps->getInput((basename + "/enableConstantTimestep").c_str(), enable_constant_dt_, false);
   tps->getInput((basename + "/dt_fixed").c_str(), constant_dt_, -1.0);
-
+  tps->getInput((basename + "/dt_initial").c_str(), initial_dt_, 1.0e-8);
+  tps->getInput((basename + "/dt_min").c_str(), minimum_dt_, 1.0e-12);
+  tps->getInput((basename + "/dt_max").c_str(), maximum_dt_, 1.0);
   tps->getInput((basename + "/bdfOrder").c_str(), bdf_order_, 3);
+  tps->getInput((basename + "/dtFactor").c_str(), factor_dt_, 0.01);
 }

src/loMach_options.hpp

@@ -72,6 +72,10 @@ class LoMachTemporalOptions {
   bool enable_constant_dt_;
 
   double cfl_;
+  double initial_dt_;
+  double minimum_dt_;
+  double maximum_dt_;
+  double factor_dt_;
   double constant_dt_;
 
   int bdf_order_;

src/thermo_chem_base.hpp

@@ -77,6 +77,8 @@ class ThermoChemModelBase {
   const spongeToThermoChem *sponge_interface_;
   const extDataToThermoChem *extData_interface_;
 
+  double thermo_pressure_;
+
  public:
   /// Destructor
   virtual ~ThermoChemModelBase() {}
@@ -159,6 +161,10 @@ class ThermoChemModelBase {
 
   void initializeFromExtData(extDataToThermoChem *extData) { extData_interface_ = extData; }
   const extDataToThermoChem *getExtDataInterface() const { return extData_interface_; }
+
+  /// Return thermodynamic pressure for restarts
+  double GetThermoPressure() { return thermo_pressure_; }
+  void SetThermoPressure(double &Po) { thermo_pressure_ = Po; }
 };
 
 /**

test/lomach-flow.test

@@ -107,7 +107,7 @@ setup() {
     [[ ${status} -eq 0 ]]
 
     # check solutions are the same!
-    h5diff -r restart_output-lid.sol.h5 restart_output-lid.sol.2.h5
+    h5diff -r --delta=1e-12 restart_output-lid.sol.h5 restart_output-lid.sol.2.h5
 
     # delete intermediate files
     rm -f restart_output-lid*.h5
@@ -118,22 +118,22 @@ setup() {
     run $EXE --runFile $RUNFILE_PIPE
     [[ ${status} -eq 0 ]]
     test -s restart_output-pipe-lam.sol.h5
-    h5diff -r --delta=1e-16 restart_output-pipe-lam.sol.h5 $REF_SOLN_PIPE /velocity
-    h5diff -r --delta=1e-16 restart_output-pipe-lam.sol.h5 $REF_SOLN_PIPE /swirl
+    h5diff -r --delta=1e-12 restart_output-pipe-lam.sol.h5 $REF_SOLN_PIPE /velocity
+    h5diff -r --delta=1e-12 restart_output-pipe-lam.sol.h5 $REF_SOLN_PIPE /swirl
 }
 
 @test "[$TEST] check that Taylor-Couette (axisymmetric, with swirl) runs and verify result" {
     run $EXE --runFile $RUNFILE_TC
     [[ ${status} -eq 0 ]]
     test -s restart_output-tc.sol.h5
-    h5diff -r --delta=1e-16 restart_output-tc.sol.h5 $REF_SOLN_TC /velocity
-    h5diff -r --delta=1e-16 restart_output-tc.sol.h5 $REF_SOLN_TC /swirl
+    h5diff -r --delta=1e-12 restart_output-tc.sol.h5 $REF_SOLN_TC /velocity
+    h5diff -r --delta=1e-12 restart_output-tc.sol.h5 $REF_SOLN_TC /swirl
 }
 
 @test "[$TEST] check pipe (axisymmetric) flow with algebraic rans model" {
     run $EXE --runFile $RUNFILE_PIPE_ARANS
     [[ ${status} -eq 0 ]]
     test -s restart_output-pipe-arans.sol.h5
-    h5diff -r --delta=1e-16 restart_output-pipe-arans.sol.h5 $REF_SOLN_PIPE_ARANS /velocity
-    h5diff -r --delta=1e-16 restart_output-pipe-arans.sol.h5 $REF_SOLN_PIPE_ARANS /swirl
+    h5diff -r --delta=1e-12 restart_output-pipe-arans.sol.h5 $REF_SOLN_PIPE_ARANS /velocity
+    h5diff -r --delta=1e-12 restart_output-pipe-arans.sol.h5 $REF_SOLN_PIPE_ARANS /swirl
 }