diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp index 5900ccf13..da5ae85d0 100644 --- a/src/calorically_perfect.hpp +++ b/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); diff --git a/src/loMach.cpp b/src/loMach.cpp index 26e91d11a..9d84893dc 100644 --- a/src/loMach.cpp +++ b/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; } } diff --git a/src/loMach.hpp b/src/loMach.hpp index d5a9aca07..d7db0bf16 100644 --- a/src/loMach.hpp +++ b/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_; diff --git a/src/loMachIO.cpp b/src/loMachIO.cpp index 489dc6768..24709f0c1 100644 --- a/src/loMachIO.cpp +++ b/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) { diff --git a/src/loMach_options.cpp b/src/loMach_options.cpp index e5cc56e8c..a78c5df3b 100644 --- a/src/loMach_options.cpp +++ b/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); } diff --git a/src/loMach_options.hpp b/src/loMach_options.hpp index b94f1b1a0..31ec965b3 100644 --- a/src/loMach_options.hpp +++ b/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_; diff --git a/src/thermo_chem_base.hpp b/src/thermo_chem_base.hpp index 9f73e502f..c51b4074e 100644 --- a/src/thermo_chem_base.hpp +++ b/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; } }; /** diff --git a/test/lomach-flow.test b/test/lomach-flow.test index ab82b5081..2da68c84e 100755 --- a/test/lomach-flow.test +++ b/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 }