diff --git a/src/algebraicSubgridModels.cpp b/src/algebraicSubgridModels.cpp index 3250a3d5f..e5643bae8 100644 --- a/src/algebraicSubgridModels.cpp +++ b/src/algebraicSubgridModels.cpp @@ -114,7 +114,7 @@ void AlgebraicSubgridModels::initializeSelf() { // gridScaleZ_.SetSize(sfes_truevsize); // resolution_gf_.SetSpace(sfes_); - if (verbose) grvy_printf(ginfo, "AlgebraicSubgridModels vectors and gf initialized...\n"); + if (rank0_) grvy_printf(ginfo, "AlgebraicSubgridModels vectors and gf initialized...\n"); // exports toFlow_interface_.eddy_viscosity = &subgridVisc_gf_; diff --git a/src/averaging.cpp b/src/averaging.cpp index 7ac3e88cb..5fe35b8cb 100644 --- a/src/averaging.cpp +++ b/src/averaging.cpp @@ -75,6 +75,9 @@ Averaging::Averaging(AveragingOptions &opts, std::string output_name) { mean_output_name_ = "mean_"; mean_output_name_.append(output_name); + + enable_mean_continuation_ = opts.enable_mean_continuation_; + zero_variances_ = opts.zero_variances_; } Averaging::~Averaging() { @@ -119,6 +122,7 @@ void Averaging::registerField(std::string name, const ParGridFunction *field_to_ // and store those fields in an AveragingFamily object that gets appended to the avg_families_ vector avg_families_.emplace_back(AveragingFamily(name, field_to_average, mean, vari, vari_start_index, vari_components)); + std::cout << "Fam size: " << avg_families_.size() << " with addition of " << name << endl; } void Averaging::initializeViz() { @@ -198,15 +202,37 @@ void Averaging::addSample(const int &iter, GasMixture *mixture) { assert(avg_families_.size() >= 1); + /* + if (!enable_mean_continuation_) { + ns_mean_ = 0; + enable_mean_continuation_ = true; + } + + if (zero_variances_) { + ns_vari_ = 0; + zero_variances_ = false; + } + */ + if (iter % sample_interval_ == 0 && iter >= step_start_mean_) { - if (iter == step_start_mean_ && rank0_) cout << "Starting mean calculation." << endl; + if (iter == step_start_mean_ && rank0_) cout << "Starting mean calculation at iter " << iter << endl; + + if (ns_mean_ == 0) { + for (size_t ifam = 0; ifam < avg_families_.size(); ifam++) { + if (avg_families_[ifam].mean_fcn_ != nullptr) { + *avg_families_[ifam].mean_fcn_ = 0.0; + } + } + } + // If we got here, then either this is our first time averaging + // or the variances have been restarted. Either way, valid to + // set variances to zero. if (ns_vari_ == 0) { - // If we got here, then either this is our first time averaging - // or the variances have been restarted. Either way, valid to - // set variances to zero. for (size_t ifam = 0; ifam < avg_families_.size(); ifam++) { - *avg_families_[ifam].vari_fcn_ = 0.0; + if (avg_families_[ifam].vari_fcn_ != nullptr) { + *avg_families_[ifam].vari_fcn_ = 0.0; + } } } @@ -251,18 +277,28 @@ void Averaging::addSampleInternal() { const double *d_inst = inst->Read(); double *d_mean = mean->ReadWrite(); - double *d_vari = vari->ReadWrite(); + double *d_vari = nullptr; + if (vari != nullptr) { + d_vari = vari->ReadWrite(); + } // Extract size information for use on device const int dof = mean->ParFESpace()->GetNDofs(); // dofs per scalar field const int neq = mean->ParFESpace()->GetVDim(); // number of scalar variables in mean field - const int d_vari_start = fam.vari_start_index_; - const int d_vari_components = fam.vari_components_; + int d_vari_start = 0; + int d_vari_components = 0; + if (vari != nullptr) { + d_vari_start = fam.vari_start_index_; + d_vari_components = fam.vari_components_; + } // Extract sample size information for use on device double d_ns_mean = (double)ns_mean_; - double d_ns_vari = (double)ns_vari_; + double d_ns_vari = 0; + if (vari != nullptr) { + d_ns_vari = (double)ns_vari_; + } // "Loop" over all dofs and update statistics MFEM_FORALL(n, dof, { @@ -296,7 +332,8 @@ void Averaging::addSampleInternal() { const double uinst_i = d_inst[n + i * dof]; const double umean_i = d_mean[n + i * dof]; delta_i = uinst_i - umean_i; - for (int j = d_vari_start + 1; j < d_vari_start + d_vari_components; j++) { + + for (int j = i + 1; j < d_vari_start + d_vari_components; j++) { const double uinst_j = d_inst[n + j * dof]; const double umean_j = d_mean[n + j * dof]; delta_j = uinst_j - umean_j; diff --git a/src/averaging.hpp b/src/averaging.hpp index 9ae4ece23..df37ddae5 100644 --- a/src/averaging.hpp +++ b/src/averaging.hpp @@ -177,6 +177,12 @@ class Averaging { /// Visualization directory (i.e., paraview dumped to mean_output_name_) std::string mean_output_name_; + /// flag to continue mean from restart + bool enable_mean_continuation_; + + /// flag to restart variances + bool zero_variances_; + /// mfem paraview data collection, used to write viz files ParaViewDataCollection *pvdc_ = nullptr; @@ -282,7 +288,10 @@ class Averaging { int GetSamplesMean() { return ns_mean_; } int GetSamplesRMS() { return ns_vari_; } int GetSamplesInterval() { return sample_interval_; } + bool ComputeMean() { return compute_mean_; } + bool ContinueMean() { return enable_mean_continuation_; } + bool RestartRMS() { return zero_variances_; } void SetSamplesMean(int &samples) { ns_mean_ = samples; } void SetSamplesRMS(int &samples) { ns_vari_ = samples; } diff --git a/src/calorically_perfect.cpp b/src/calorically_perfect.cpp index d4be8ef89..786cc37d8 100644 --- a/src/calorically_perfect.cpp +++ b/src/calorically_perfect.cpp @@ -47,6 +47,10 @@ using namespace mfem; using namespace mfem::common; +/// forward declarations +double temp_rt3d(const Vector &x, double t); +double temp_channel(const Vector &x, double t); + MFEM_HOST_DEVICE double Sutherland(const double T, const double mu_star, const double T_star, const double S_star) { const double T_rat = T / T_star; const double T_rat_32 = T_rat * sqrt(T_rat); @@ -239,6 +243,9 @@ void CaloricallyPerfectThermoChem::initializeSelf() { std::cout << "exports set..." << endl; } + tpsP_->getInput("loMach/calperfect/numerical-integ", numerical_integ_, true); + tpsP_->getInput("loMach/calperfect/over-integ", over_integrate_, false); + //----------------------------------------------------- // 2) Set the initial condition //----------------------------------------------------- @@ -256,9 +263,26 @@ void CaloricallyPerfectThermoChem::initializeSelf() { // 2) For restarts, this IC is overwritten by the restart field, // which is read later. - ConstantCoefficient t_ic_coef; - t_ic_coef.constant = T_ic_; - Tn_gf_.ProjectCoefficient(t_ic_coef); + tpsP_->getInput("loMach/calperfect/ic", ic_string_, std::string("")); + + // set IC if we have one at this point + if (!ic_string_.empty()) { + if (ic_string_ == "rt3D") { + if (rank0_) std::cout << "Setting rt3D IC..." << std::endl; + FunctionCoefficient t_excoeff(temp_rt3d); + t_excoeff.SetTime(0.0); + Tn_gf_.ProjectCoefficient(t_excoeff); + } else if (ic_string_ == "channel") { + if (rank0_) std::cout << "Setting channel IC..." << std::endl; + FunctionCoefficient t_excoeff(temp_channel); + t_excoeff.SetTime(0.0); + Tn_gf_.ProjectCoefficient(t_excoeff); + } + } else { + ConstantCoefficient t_ic_coef; + t_ic_coef.constant = T_ic_; + Tn_gf_.ProjectCoefficient(t_ic_coef); + } Tn_gf_.GetTrueDofs(Tn_); Tnm1_gf_.SetFromTrueDofs(Tn_); @@ -334,16 +358,13 @@ void CaloricallyPerfectThermoChem::initializeSelf() { Array attr_outlet(pmesh_->bdr_attributes.Max()); attr_outlet = 0; - // No code for this yet, so die if detected - // assert(numOutlets == 0); - - // But... outlets will just get homogeneous Neumann on T, so + // Outlets will just get homogeneous Neumann on T, so // basically need to do nothing. } // Wall BCs { - std::cout << "There are " << pmesh_->bdr_attributes.Max() << " boundary attributes!" << std::endl; + if (rank0_) std::cout << "There are " << pmesh_->bdr_attributes.Max() << " boundary attributes" << std::endl; Array attr_wall(pmesh_->bdr_attributes.Max()); attr_wall = 0; @@ -356,7 +377,7 @@ void CaloricallyPerfectThermoChem::initializeSelf() { tpsP_->getRequiredInput((basepath + "/type").c_str(), type); if (type == "viscous_isothermal") { - std::cout << "Adding patch = " << patch << " to isothermal wall list!" << std::endl; + if (rank0_) std::cout << "Adding patch = " << patch << " to isothermal wall list" << std::endl; attr_wall = 0; attr_wall[patch - 1] = 1; @@ -391,9 +412,17 @@ void CaloricallyPerfectThermoChem::initializeOperators() { // unsteady: p+p [+p] = 2p [3p] // convection: p+p+(p-1) [+p] = 3p-1 [4p-1] // diffusion: (p-1)+(p-1) [+p] = 2p-2 [3p-2] - const IntegrationRule &ir_i = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), 2 * order_ + 1); - const IntegrationRule &ir_nli = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), 4 * order_); - const IntegrationRule &ir_di = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), 3 * order_ - 1); + int nir_i = 2 * order_ - 1; + int nir_nli = 2 * order_ - 1; + int nir_di = 2 * order_ - 1; + if (over_integrate_) { + nir_i = 3 * order_ + 1; + nir_nli = 4 * order_; + nir_di = 3 * order_ - 1; + } + const IntegrationRule &ir_i = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), nir_i); + const IntegrationRule &ir_nli = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), nir_nli); + const IntegrationRule &ir_di = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), nir_di); if (rank0_) std::cout << "Integration rules set" << endl; // coefficients for operators @@ -733,6 +762,18 @@ void CaloricallyPerfectThermoChem::initializeViz(ParaViewDataCollection &pvdc) { pvdc.RegisterField("Qt", &Qt_gf_); } +void CaloricallyPerfectThermoChem::initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) { + if (average.ComputeMean()) { + // fields for averaging + average.registerField(std::string("temperature"), &Tn_gf_, false, 0, 1); + + // io init + io.registerIOFamily("Time-averaged temperature", "/meanTemp", average.GetMeanField(std::string("temperature")), + false, continuation, sfec_); + io.registerIOVar("/meanTemp", "", 0, true); + } +} + /** Update boundary conditions for Temperature, useful for ramping a dirichlet condition over time @@ -1426,3 +1467,37 @@ double temp_inlet(const Vector &coords, double t) { return temp; } #endif + +double temp_rt3d(const Vector &x, double t) { + double CC = 0.05; + double twoPi = 6.28318530718; + double yWidth = 0.1; + double yInt, dy, wt; + double temp, dT; + double Tlo = 100.0; + double Thi = 1500.0; + + yInt = std::cos(twoPi * x[0]) + std::cos(twoPi * x[2]); + yInt *= CC; + yInt += 4.0; + + dy = x[1] - yInt; + dT = Thi - Tlo; + + wt = 0.5 * (tanh(-dy / yWidth) + 1.0); + temp = Tlo + wt * dT; + + return temp; +} + +/// Used to set the channel IC +double temp_channel(const Vector &x, double t) { + double Thi = 400.0; + double Tlo = 200.0; + double temp = 300.0; + + // expects channel height (-1,1) + temp = Tlo + (Thi - Tlo) * 0.5 * (x(1) + 1.0); + + return temp; +} diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp index e1c04b714..c9e9e6b65 100644 --- a/src/calorically_perfect.hpp +++ b/src/calorically_perfect.hpp @@ -43,6 +43,7 @@ class Tps; #include +#include "averaging.hpp" #include "dirichlet_bc_helper.hpp" #include "io.hpp" #include "thermo_chem_base.hpp" @@ -70,10 +71,13 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { double dt_; double time_; + std::string ic_string_; + // Flags bool rank0_; /**< true if this is rank 0 */ bool partial_assembly_ = false; /**< Enable/disable partial assembly of forms. */ bool numerical_integ_ = true; /**< Enable/disable numerical integration rules of forms. */ + bool over_integrate_ = false; /**< Enable/disable numerical integration with higher order. */ bool constant_viscosity_ = false; /**< Enable/disable constant viscosity */ bool constant_density_ = false; /**< Enable/disable constant density */ bool domain_is_open_ = false; /**< true if domain is open */ @@ -215,6 +219,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { void step() final; void initializeIO(IODataOrganizer &io) final; void initializeViz(ParaViewDataCollection &pvdc) final; + void initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) final; void screenHeader(std::vector &header) const final; void screenValues(std::vector &values) final; diff --git a/src/loMach.cpp b/src/loMach.cpp index 2fdbdfadf..ee8979a2a 100644 --- a/src/loMach.cpp +++ b/src/loMach.cpp @@ -90,6 +90,8 @@ LoMachSolver::~LoMachSolver() { delete thermo_; delete sponge_; delete turbModel_; + delete avg_opts_; + delete average_; delete meshData_; // allocated in constructor @@ -186,6 +188,11 @@ void LoMachSolver::initialize() { exit(ERROR); } + // Instantiate averaging + avg_opts_ = new AveragingOptions(); + avg_opts_->read(tpsP_); + average_ = new Averaging(*avg_opts_, loMach_opts_.io_opts_.output_dir_); + // Initialize time marching coefficients. NB: dt will be reset // prior to time step, but must be initialized here in order to // avoid possible uninitialized usage when constructing operators in @@ -209,25 +216,6 @@ void LoMachSolver::initialize() { flow_->initializeSelf(); thermo_->initializeSelf(); - // Initialize restart read/write capability - flow_->initializeIO(ioData); - thermo_->initializeIO(ioData); - - const bool restart_serial = - (loMach_opts_.io_opts_.restart_serial_read_ || loMach_opts_.io_opts_.restart_serial_write_); - ioData.initializeSerial(rank0_, restart_serial, serial_mesh_, meshData_->getLocalToGlobalElementMap(), - &partitioning_); - // MPI_Barrier(groupsMPI->getTPSCommWorld()); - if (verbose) grvy_printf(ginfo, "ioData.init thingy...\n"); - - // 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 turbModel_->initializeFromThermoChem(&thermo_->toTurbModel_interface_); turbModel_->initializeFromFlow(&flow_->toTurbModel_interface_); @@ -246,9 +234,14 @@ void LoMachSolver::initialize() { turbModel_->setup(); turbModel_->initializeOperators(); thermo_->initializeOperators(); - // if(rank0_) {std::cout << "check: ops set..." << endl;} - // TODO(trevilo): Enable averaging. See note in loMach.hpp + // Initialize restart read/write capability + flow_->initializeIO(ioData); + thermo_->initializeIO(ioData); + + // Initialize statistics + flow_->initializeStats(*average_, ioData, average_->ContinueMean()); + thermo_->initializeStats(*average_, ioData, average_->ContinueMean()); // Initialize visualization pvdc_ = new ParaViewDataCollection(loMach_opts_.io_opts_.output_dir_, pmesh_); @@ -262,6 +255,20 @@ void LoMachSolver::initialize() { thermo_->initializeViz(*pvdc_); sponge_->initializeViz(*pvdc_); extData_->initializeViz(*pvdc_); + average_->initializeViz(); + + const bool restart_serial = + (loMach_opts_.io_opts_.restart_serial_read_ || loMach_opts_.io_opts_.restart_serial_write_); + ioData.initializeSerial(rank0_, restart_serial, serial_mesh_, meshData_->getLocalToGlobalElementMap(), + &partitioning_); + + // If restarting, read restart files + if (loMach_opts_.io_opts_.enable_restart_) { + restart_files_hdf5("read"); + if (thermoPressure_ > 0.0) { + thermo_->SetThermoPressure(thermoPressure_); + } + } } void LoMachSolver::UpdateTimestepHistory(double dt) { @@ -397,6 +404,13 @@ void LoMachSolver::solveStep() { } } + // averages + if (avg_opts_->sample_interval_ != 0) { + if (iter % avg_opts_->sample_interval_ == 0 && iter != 0) { + average_->addSample(iter, nullptr); + } + } + // restart files if (iter % loMach_opts_.output_frequency_ == 0 && iter != 0) { thermoPressure_ = thermo_->GetThermoPressure(); @@ -408,6 +422,7 @@ void LoMachSolver::solveStep() { pvdc_->SetCycle(iter); pvdc_->SetTime(temporal_coeff_.time); pvdc_->Save(); + average_->writeViz(iter, temporal_coeff_.time, avg_opts_->save_mean_history_); } // check for DIE @@ -443,6 +458,7 @@ void LoMachSolver::solveEnd() { pvdc_->SetTime(temporal_coeff_.time); if (rank0_ == true) std::cout << " Saving final step to paraview: " << iter << "... " << endl; pvdc_->Save(); + average_->writeViz(iter, temporal_coeff_.time, avg_opts_->save_mean_history_); MPI_Barrier(groupsMPI->getTPSCommWorld()); if (rank0_ == true) std::cout << " ...complete!" << endl; } @@ -651,7 +667,7 @@ void LoMachSolver::setTimestep() { temporal_coeff_.dt = dt_initial; } - std::cout << "dt from setTimestep: " << temporal_coeff_.dt << " max_speed: " << max_speed << endl; + if (rank0_) std::cout << "dt from setTimestep: " << temporal_coeff_.dt << " max_speed: " << max_speed << endl; } } diff --git a/src/loMach.hpp b/src/loMach.hpp index 4492423ba..7de29b8a6 100644 --- a/src/loMach.hpp +++ b/src/loMach.hpp @@ -61,6 +61,7 @@ class Tps; #include +#include "averaging.hpp" #include "externalData_base.hpp" #include "io.hpp" #include "loMach_options.hpp" @@ -127,6 +128,8 @@ class LoMachSolver : public TPS::PlasmaSolver { FlowBase *flow_ = nullptr; SpongeBase *sponge_ = nullptr; ExternalDataBase *extData_ = nullptr; + AveragingOptions *avg_opts_ = nullptr; + Averaging *average_ = nullptr; // Mesh and geometry related ParMesh *pmesh_ = nullptr; diff --git a/src/loMachIO.cpp b/src/loMachIO.cpp index 6ea124771..e2c7997af 100644 --- a/src/loMachIO.cpp +++ b/src/loMachIO.cpp @@ -61,6 +61,13 @@ void LoMachSolver::write_restart_files_hdf5(hid_t file, bool serialized_write) { // thermodynamic pressure h5_save_attribute(file, "Po", thermoPressure_); + if (average_->ComputeMean()) { + // samples meanUp + h5_save_attribute(file, "samplesMean", average_->GetSamplesMean()); + h5_save_attribute(file, "samplesRMS", average_->GetSamplesRMS()); + h5_save_attribute(file, "samplesInterval", average_->GetSamplesInterval()); + } + // code revision #ifdef BUILD_VERSION { @@ -116,6 +123,25 @@ void LoMachSolver::read_restart_files_hdf5(hid_t file, bool serialized_read) { if (H5Aexists(file, "Po")) { h5_read_attribute(file, "Po", thermoPressure_); } + if (average_->ComputeMean() && average_->ContinueMean()) { + int samplesMean, samplesRMS, intervals; + h5_read_attribute(file, "samplesMean", samplesMean); + h5_read_attribute(file, "samplesInterval", intervals); + average_->SetSamplesMean(samplesMean); + average_->SetSamplesInterval(intervals); + + int istatus = H5Aexists(file, "samplesRMS"); + if (istatus > 0) { + h5_read_attribute(file, "samplesRMS", samplesRMS); + } else { + samplesRMS = samplesMean; + } + + if (average_->RestartRMS() == true) { + samplesRMS = 0; + } + average_->SetSamplesRMS(samplesRMS); + } } if (serialized_read) { diff --git a/src/split_flow_base.hpp b/src/split_flow_base.hpp index 84032fd47..8241ce453 100644 --- a/src/split_flow_base.hpp +++ b/src/split_flow_base.hpp @@ -38,6 +38,7 @@ #include "tps_mfem_wrap.hpp" class IODataOrganizer; +class Averaging; struct thermoChemToFlow; struct turbModelToFlow; struct spongeToFlow; @@ -90,6 +91,7 @@ class FlowBase { virtual void initializeIO(IODataOrganizer &io) const {} virtual void initializeViz(mfem::ParaViewDataCollection &pvdc) const {} + virtual void initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) const {} /** * @brief Header strings for screen dump diff --git a/src/thermo_chem_base.hpp b/src/thermo_chem_base.hpp index 098642f9d..78477259c 100644 --- a/src/thermo_chem_base.hpp +++ b/src/thermo_chem_base.hpp @@ -42,6 +42,7 @@ class Tps; } class IODataOrganizer; +class Averaging; struct flowToThermoChem; struct turbModelToThermoChem; struct spongeToThermoChem; @@ -123,6 +124,11 @@ class ThermoChemModelBase { */ virtual void initializeViz(mfem::ParaViewDataCollection &pvdc) {} + /** + * @brief Hook to let averaging register fields and restart fields with the IODataOrganizer. + */ + virtual void initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) {} + /** * @brief Header strings for screen dump * diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index a8d419d09..b255c109f 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -46,7 +46,9 @@ using namespace mfem; /// forward declarations void vel_exact_tgv2d(const Vector &x, double t, Vector &u); +void vel_tgv2d_uniform(const Vector &x, double t, Vector &u); void vel_exact_pipe(const Vector &x, double t, Vector &u); +void vel_channel(const Vector &x, double t, Vector &u); static double radius(const Vector &pos) { return pos[0]; } FunctionCoefficient radius_coeff(radius); @@ -112,6 +114,7 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS // NOTE: this should default to false as it is generally not safe, but much of the // test results rely on it being true tps->getInput("loMach/tomboulides/numerical-integ", numerical_integ_, true); + tps->getInput("loMach/tomboulides/over-integ", over_integrate_, false); // Can't use numerical integration with axisymmetric b/c it // locates quadrature points on the axis, which can lead to @@ -386,10 +389,20 @@ void Tomboulides::initializeSelf() { // set IC if we have one at this point if (!ic_string_.empty()) { if (ic_string_ == "tgv2d") { - std::cout << "Setting tgv2d IC..." << std::endl; + if (rank0_) std::cout << "Setting tgv2d IC..." << std::endl; VectorFunctionCoefficient u_excoeff(2, vel_exact_tgv2d); u_excoeff.SetTime(0.0); u_curr_gf_->ProjectCoefficient(u_excoeff); + } else if (ic_string_ == "tgv2d_uniform") { + if (rank0_) std::cout << "Setting tgv2d+uniform IC..." << std::endl; + VectorFunctionCoefficient u_excoeff(2, vel_tgv2d_uniform); + u_excoeff.SetTime(0.0); + u_curr_gf_->ProjectCoefficient(u_excoeff); + } else if (ic_string_ == "channel") { + if (rank0_) std::cout << "Setting channel IC..." << std::endl; + VectorFunctionCoefficient u_excoeff(3, vel_channel); + u_excoeff.SetTime(0.0); + u_curr_gf_->ProjectCoefficient(u_excoeff); } } @@ -620,8 +633,18 @@ void Tomboulides::initializeOperators() { // Gauss-Lobatto quad pts correspond to the Gauss-Lobatto nodes. // For most terms this will result in an under-integration, but it // has the nice consequence that the mass matrix is diagonal. - const IntegrationRule &ir_ni_v = gll_rules.Get(vfes_->GetFE(0)->GetGeomType(), 2 * vorder_ - 1); - const IntegrationRule &ir_ni_p = gll_rules.Get(pfes_->GetFE(0)->GetGeomType(), 2 * porder_ - 1); + int nir_v = 2 * vorder_ - 1; + int nir_p = 2 * porder_ - 1; + if (over_integrate_) { + int nir_iv = 3 * vorder_ + 1; + int nir_nliv = 4 * vorder_; + int nir_ip = 3 * porder_ + 1; + int nir_nlip = 4 * porder_; + nir_v = std::max(nir_iv, nir_nliv); + nir_p = std::max(nir_ip, nir_nlip); + } + const IntegrationRule &ir_ni_v = gll_rules.Get(vfes_->GetFE(0)->GetGeomType(), nir_v); + const IntegrationRule &ir_ni_p = gll_rules.Get(pfes_->GetFE(0)->GetGeomType(), nir_p); // Empty array, use where we want operators without BCs Array empty; @@ -1032,6 +1055,44 @@ void Tomboulides::initializeViz(mfem::ParaViewDataCollection &pvdc) const { } } +void Tomboulides::initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) const { + if (average.ComputeMean()) { + // fields for averaging + average.registerField(std::string("velocity"), u_curr_gf_, true, 0, nvel_); + average.registerField(std::string("pressure"), p_gf_, false, 0, 1); + + // io init + io.registerIOFamily("Time-averaged velocity", "/meanVel", average.GetMeanField(std::string("velocity")), false, + continuation, vfec_); + io.registerIOVar("/meanVel", "", 0, true); + if (dim_ >= 2) io.registerIOVar("/meanVel", "", 1, true); + if (dim_ == 3) io.registerIOVar("/meanVel", "", 2, true); + + io.registerIOFamily("Time-averaged pressure", "/meanPres", average.GetMeanField(std::string("pressure")), false, + continuation, pfec_); + io.registerIOVar("/meanPres", "

", 0, true); + + // rms + io.registerIOFamily("RMS velocity fluctuation", "/rmsData", average.GetVariField(std::string("velocity")), false, + continuation, vfec_); + if (nvel_ == 3) { + io.registerIOVar("/rmsData", "uu", 0, true); + io.registerIOVar("/rmsData", "vv", 1, true); + io.registerIOVar("/rmsData", "ww", 2, true); + io.registerIOVar("/rmsData", "uv", 3, true); + io.registerIOVar("/rmsData", "uw", 4, true); + io.registerIOVar("/rmsData", "vw", 5, true); + } else if (nvel_ == 2) { + io.registerIOVar("/rmsData", "uu", 0, true); + io.registerIOVar("/rmsData", "vv", 1, true); + io.registerIOVar("/rmsData", "uv", 2, true); + } else { + // only nvel = 2 or 3 supported + assert(false); + } + } +} + void Tomboulides::step() { //------------------------------------------------------------------------ // The time step is split into 4 large chunks: @@ -1634,3 +1695,57 @@ void vel_exact_pipe(const Vector &x, double t, Vector &u) { u(0) = 0.0; u(1) = 2.0 * (1 - x[0] * x[0]); } + +/// Used to set the velocity IC with TG field and uniform +void vel_tgv2d_uniform(const Vector &x, double t, Vector &u) { + const double u0 = 1.0; + const double F = 0.1; + const double PI = 3.14159265359; + double twoPi = 2.0 * PI; + + u(0) = u0; + u(1) = 0.0; + + u(0) += +F * std::sin(twoPi * x[0]) * std::cos(twoPi * x[1]); + u(1) += -F * std::cos(twoPi * x[0]) * std::sin(twoPi * x[1]); +} + +/// Used to set the channel IC +void vel_channel(const Vector &x, double t, Vector &u) { + double PI = 3.14159265359; + double Lx = 25.0; + double Ly = 2.0; + double Lz = 9.4; + double Umean = 1.0; + double uInt = 0.1; + int nModes = 4; + double uM; + double ax, by, cz; + double AA, BB, CC; + double wall; + + // expects channel height (-1,1) + wall = (1.0 - std::pow(x(1), 8.0)); + u(0) = Umean * wall; + u(1) = 0.0; + u(2) = 0.0; + + for (int n = 1; n <= nModes; n++) { + ax = 4.0 * PI / Lx * (double)n; + by = 2.0 * PI / Ly * (double)n; + cz = 2.0 * PI / Lz * (double)n; + + AA = 1.0; + BB = 1.0; + CC = -(AA * ax + BB * by) / cz; + + uM = uInt / (double)n; + + u(0) += uM * AA * cos(ax * (x(0) + (double)(n - 1) * Umean)) * sin(by * x(1)) * + sin(cz * (x(2) + 0.5 * (double)(n - 1) * Umean)) * wall; + u(1) += uM * BB * sin(ax * (x(0) + (double)(n - 1) * Umean)) * cos(by * x(1)) * + sin(cz * (x(2) + 0.5 * (double)(n - 1) * Umean)) * wall; + u(2) += uM * CC * sin(ax * (x(0) + (double)(n - 1) * Umean)) * sin(by * x(1)) * + cos(cz * (x(2) + 0.5 * (double)(n - 1) * Umean)) * wall; + } +} diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp index 9b525a899..fc4953bd2 100644 --- a/src/tomboulides.hpp +++ b/src/tomboulides.hpp @@ -82,6 +82,7 @@ class Tomboulides final : public FlowBase { // Options // TODO(trevilo): hardcoded for testing. Need to set based on input file. bool numerical_integ_ = false; + bool over_integrate_ = false; bool partial_assembly_ = false; int pressure_solve_pl_ = 0; @@ -353,6 +354,11 @@ class Tomboulides final : public FlowBase { */ void initializeViz(mfem::ParaViewDataCollection &pvdc) const final; + /** + * @brief Initialize statistics outputs + */ + void initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) const final; + /// Advance void step() final; diff --git a/test/Makefile.am b/test/Makefile.am index 5714df0b7..1a4cd446d 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -31,6 +31,7 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- ref_solns/lequere-varmu/*.h5 \ ref_solns/taylor-couette/*.h5 \ ref_solns/pipe/*.h5 \ + ref_solns/aveLoMach/*.h5 \ vpath.sh die.sh count_gpus.sh sniff_mpirun.sh \ cyl3d.gpu.test cyl3d.mflow.gpu.test wedge.gpu.test \ averaging.gpu.test cyl3d.test cyl3d.gpu.python.test cyl3d.mflow.test cyl3d.dtconst.test \ @@ -45,7 +46,7 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- tabulated.test lte_mixture.test distance_fcn.test \ sgsSmag.test sgsSigma.test heatEq.test sponge.test plate.test pipe.mix.test lte2noneq-restart.test \ coupled-3d.interface.test plasma.axisym.test plasma.axisym.lte1d.test \ - lomach-flow.test lomach-lequere.test interpInlet.test sgsLoMach.test autoPeriodic.test + lomach-flow.test lomach-lequere.test interpInlet.test sgsLoMach.test autoPeriodic.test aveLoMach.test TESTS = vpath.sh XFAIL_TESTS = @@ -110,7 +111,8 @@ TESTS += cyl3d.test \ lomach-flow.test \ lomach-lequere.test \ interpInlet.test \ - autoPeriodic.test + autoPeriodic.test \ + aveLoMach.test if PYTHON_ENABLED TESTS += cyl3d.python.test \ diff --git a/test/aveLoMach.test b/test/aveLoMach.test new file mode 100755 index 000000000..31af36add --- /dev/null +++ b/test/aveLoMach.test @@ -0,0 +1,30 @@ +#!./bats +# -*- mode: sh -*- + +TEST="aveLoMach" +RUNFILE="inputs/input.aveLoMach.ini" +EXE="../src/tps" +RESTART="ref_solns/aveLoMach/restart_output.sol.h5" + +setup() { + SOLN_FILE=restart_output.sol.h5 + REF_FILE=ref_solns/aveLoMach/restart_output.sol.h5 +} + +@test "[$TEST] check for input file $RUNFILE" { + test -s $RUNFILE +} + +@test "[$TEST] run tps with input -> $RUNFILE" { + rm -rf output/* + rm -f $SOLN_FILE + run $EXE --runFile $RUNFILE + [[ ${status} -eq 0 ]] + test -s $SOLN_FILE +} + +@test "[$TEST] verify tps output with input -> $RUNFILE" { + test -s $SOLN_FILE + test -s $REF_FILE + h5diff -r --delta=1e-12 $SOLN_FILE $REF_FILE meanVel +} diff --git a/test/inputs/input.aveLoMach.ini b/test/inputs/input.aveLoMach.ini new file mode 100644 index 000000000..aa72f0148 --- /dev/null +++ b/test/inputs/input.aveLoMach.ini @@ -0,0 +1,83 @@ +[solver] +type = loMach + +[loMach] +flow-solver = tomboulides +thermo-solver = calorically-perfect +mesh = meshes/flatBox.msh +order = 1 +maxIters = 100 +outputFreq = 100 +fluid = dry_air +refLength = 1.0 +equation_system = navier-stokes +enablePressureForcing = False +enableGravity = False +openSystem = True +ambientPressure = 101326. +sgsModel = none +#sgsModel = smagorinsky +sgsModelConstant = 0.09 + +[loMach/calperfect] +#viscosity-model = constant +#density-model = constant +#constant-visc/mu0 = 1.552e-5 +#constant-density = 1.1839 +linear-solver-rtol = 1.0e-12 +linear-solver-max-iter = 2000 +viscosity-model = sutherland +sutherland/mu0 = 1.68e-5 +sutherland/T0 = 273.0 +sutherland/S0 = 110.4 +numerical-integ = false +#Prandtl = 0.72 +Prandtl = 0.01 + +[loMach/tomboulides] +linear-solver-rtol = 1.0e-12 +linear-solver-max-iter = 2000 +numerical-integ = false +ic = tgv2d_uniform + +[io] +outdirBase = output +#enableRestart = True +#restartMode = variableP + +[time] +integrator = curlcurl +cfl = 0.5 +#dt_fixed = 1.0e-4 +dt_initial = 1.0e-4 +dtFactor = 0.01 + +[initialConditions] +rho = 1.1839 +rhoU = 0. +rhoV = 0. +rhoW = 0. +temperature = 298.15 +pressure = 101325.0 + +[averaging] +enableContinuation = False +#enableContinuation = True +#restartRMS = True +saveMeanHist = True +startIter = 1 +sampleFreq = 10 + +[boundaryConditions] +numWalls = 0 +numInlets = 0 +numOutlets = 0 + +[periodicity] +enablePeriodic = True +periodicX = True +periodicY = True +periodicZ = True +#xTrans = 1.0 +#yTrans = 1.0 +#zTrans = 0.1 \ No newline at end of file diff --git a/test/inputs/input.heatedBox.ini b/test/inputs/input.heatedBox.ini index c3b68817a..80996ced4 100644 --- a/test/inputs/input.heatedBox.ini +++ b/test/inputs/input.heatedBox.ini @@ -27,6 +27,7 @@ viscosity-model = sutherland sutherland/mu0 = 1.68e-5 sutherland/T0 = 273.0 sutherland/S0 = 110.4 +numerical-integ = false [io] outdirBase = output diff --git a/test/ref_solns/aveLoMach/restart_output.sol.h5 b/test/ref_solns/aveLoMach/restart_output.sol.h5 new file mode 100644 index 000000000..238cd7c2e Binary files /dev/null and b/test/ref_solns/aveLoMach/restart_output.sol.h5 differ