From 39407ac79a0e1be596939b199b09c709eef3b89f Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Fri, 28 Jun 2024 14:52:39 -0700 Subject: [PATCH 01/12] explicit streamwise-diffusion addition. untested but building --- src/calorically_perfect.cpp | 90 +++++++++++++++++++++++++++++++++- src/calorically_perfect.hpp | 23 ++++++++- src/loMach.cpp | 6 ++- src/tomboulides.cpp | 70 +++++++++++++++++++++++++- src/tomboulides.hpp | 20 +++++++- src/utils.cpp | 97 +++++++++++++++++++++++++++++++++++++ src/utils.hpp | 2 +- 7 files changed, 299 insertions(+), 9 deletions(-) diff --git a/src/calorically_perfect.cpp b/src/calorically_perfect.cpp index c182c14e6..d07324157 100644 --- a/src/calorically_perfect.cpp +++ b/src/calorically_perfect.cpp @@ -58,10 +58,12 @@ MFEM_HOST_DEVICE double Sutherland(const double T, const double mu_star, const d } CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, - temporalSchemeCoefficients &time_coeff, TPS::Tps *tps) - : tpsP_(tps), pmesh_(pmesh), time_coeff_(time_coeff) { + temporalSchemeCoefficients &time_coeff, + ParGridFunction *gridScale, TPS::Tps *tps) + : tpsP_(tps), pmesh_(pmesh), dim_(pmesh->Dimension()), time_coeff_(time_coeff) { rank0_ = (pmesh_->GetMyRank() == 0); order_ = loMach_opts->order; + gridScale_gf_ = *gridScale; std::string visc_model; tpsP_->getInput("loMach/calperfect/viscosity-model", visc_model, std::string("sutherland")); @@ -132,6 +134,10 @@ CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, tpsP_->getInput("loMach/calperfect/msolve-atol", mass_inverse_atol_, default_atol_); tpsP_->getInput("loMach/calperfect/msolve-max-iter", mass_inverse_max_iter_, max_iter_); tpsP_->getInput("loMach/calperfect/msolve-verbosity", mass_inverse_pl_, pl_solve_); + + tps->getInput("loMach/tomboulides/streamwise-stabilization", sw_stab_, false); + tps->getInput("loMach/tomboulides/Reh_offset", re_offset_, 1.0); + tps->getInput("loMach/tomboulides/Reh_factor", re_factor_, 0.01); } CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() { @@ -151,6 +157,7 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() { delete MsRho_form_; delete Ms_form_; delete At_form_; + delete D_form_; delete rhou_coeff_; delete rhon_next_coeff_; delete un_next_coeff_; @@ -167,6 +174,8 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() { // allocated in initializeSelf delete sfes_; delete sfec_; + delete vfes_; + delete vfec_; } void CaloricallyPerfectThermoChem::initializeSelf() { @@ -178,6 +187,9 @@ void CaloricallyPerfectThermoChem::initializeSelf() { sfec_ = new H1_FECollection(order_); sfes_ = new ParFiniteElementSpace(pmesh_, sfec_); + vfec_ = new H1_FECollection(order_, dim_); + vfes_ = new ParFiniteElementSpace(pmesh_, vfec_, dim_); + // Check if fully periodic mesh if (!(pmesh_->bdr_attributes.Size() == 0)) { temp_ess_attr_.SetSize(pmesh_->bdr_attributes.Max()); @@ -189,6 +201,7 @@ void CaloricallyPerfectThermoChem::initializeSelf() { if (rank0_) grvy_printf(ginfo, "CaloricallyPerfectThermoChem paces constructed...\n"); int sfes_truevsize = sfes_->GetTrueVSize(); + int vfes_truevsize = vfes_->GetTrueVSize(); Qt_.SetSize(sfes_truevsize); Qt_ = 0.0; @@ -238,8 +251,16 @@ void CaloricallyPerfectThermoChem::initializeSelf() { visc_gf_.SetSpace(sfes_); visc_gf_ = 0.0; + vel_gf_.SetSpace(vfes_); + tmpR1_gf_.SetSpace(vfes_); + tmpR0_gf_.SetSpace(sfes_); + + swDiff_.SetSize(sfes_truevsize); tmpR0_.SetSize(sfes_truevsize); + tmpR0a_.SetSize(sfes_truevsize); tmpR0b_.SetSize(sfes_truevsize); + tmpR0c_.SetSize(sfes_truevsize); + tmpR1_.SetSize(vfes_truevsize); R0PM0_gf_.SetSpace(sfes_); @@ -471,6 +492,24 @@ void CaloricallyPerfectThermoChem::initializeOperators() { Ms_form_->Assemble(); Ms_form_->FormSystemMatrix(empty, Ms_); + // Divergence operator + D_form_ = new ParMixedBilinearForm(vfes_, sfes_); + VectorDivergenceIntegrator *vd_mblfi; + // if (axisym_) { + // vd_mblfi = new VectorDivergenceIntegrator(radius_coeff); + // } else { + vd_mblfi = new VectorDivergenceIntegrator(); + // } + if (numerical_integ_) { + vd_mblfi->SetIntRule(&ir_nli); + } + D_form_->AddDomainIntegrator(vd_mblfi); + if (partial_assembly_) { + D_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL); + } + D_form_->Assemble(); + D_form_->FormRectangularSystemMatrix(empty, empty, D_op_); + // mass matrix with rho MsRho_form_ = new ParBilinearForm(sfes_); auto *msrho_blfi = new MassIntegrator(*rho_coeff_); @@ -666,6 +705,12 @@ void CaloricallyPerfectThermoChem::step() { tmpR0_ = (dtP_ / Cp_); Ms_->AddMult(tmpR0_, resT_); + // Add streamwise stability to rhs + if (sw_stab_) { + streamwiseDiffusion(Tn_, swDiff_); + resT_.Add(1.0, swDiff_); + } + // Add natural boundary terms here later // NB: adiabatic natural BC is handled, but don't have ability to impose non-zero heat flux yet @@ -1022,3 +1067,44 @@ void CaloricallyPerfectThermoChem::screenValues(std::vector &values) { values[0] = thermo_pressure_ / ambient_pressure_; } } + +void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { + // compute streamwise gradient of input field + tmpR0_gf_.SetFromTrueDofs(phi); + (flow_interface_->velocity)->GetTrueDofs(tmpR0a_); + vel_gf_.SetFromTrueDofs(tmpR0a_); + streamwiseGrad(tmpR0_gf_, vel_gf_, tmpR1_gf_); + + // divergence of sw-grad + tmpR1_gf_.GetTrueDofs(tmpR1_); + D_form_->Mult(tmpR1_, swDiff); + + gridScale_gf_.GetTrueDofs(tmpR0b_); + (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + + const double *rho = rn_.HostRead(); + const double *vel = tmpR0a_.HostRead(); + const double *del = tmpR0b_.HostRead(); + const double *mu = visc_.HostRead(); + const double *muT = tmpR0c_.HostRead(); + double *data = swDiff.HostReadWrite(); + + int Sdof = rn_.Size(); + for (int dof = 0; dof < Sdof; dof++) { + double Umag = 0.0; + for (int i = 0; i < dim_; i++) Umag += vel[i] * vel[i]; + Umag = std::sqrt(Umag); + + // element Re + double Re = Umag * del[dof] * rho[dof] / (mu[dof] + muT[dof]); + + // SUPG weight + double Csupg = 0.5 * (tanh(re_factor_ * Re - re_offset_) + 1.0); + + // streamwise diffusion coeff + double CswDiff = Csupg * Umag * del[dof] * rho[dof]; + + // scaled streamwise Laplacian + data[dof] *= CswDiff; + } +} diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp index 118cd3dd5..3851b2881 100644 --- a/src/calorically_perfect.hpp +++ b/src/calorically_perfect.hpp @@ -65,6 +65,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { // Mesh and discretization scheme info ParMesh *pmesh_ = nullptr; + int dim_; int order_; IntegrationRules gll_rules_; const temporalSchemeCoefficients &time_coeff_; @@ -110,6 +111,10 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { double hsolve_rtol_; double hsolve_atol_; + bool sw_stab_; + double re_offset_; + double re_factor_; + // Boundary condition info Array temp_ess_attr_; /**< List of patches with Dirichlet BC on temperature */ Array Qt_ess_attr_; /**< List of patches with Dirichlet BC on Q (thermal divergence) */ @@ -144,6 +149,10 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { // Scalar \f$H^1\f$ finite element space. ParFiniteElementSpace *sfes_ = nullptr; + // Vector fe collection and space + FiniteElementCollection *vfec_ = nullptr; + ParFiniteElementSpace *vfes_ = nullptr; + // Fields ParGridFunction Tnm1_gf_, Tnm2_gf_; ParGridFunction Tn_gf_, Tn_next_gf_, Text_gf_, resT_gf_; @@ -155,6 +164,11 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { ParGridFunction R0PM0_gf_; ParGridFunction Qt_gf_; + ParGridFunction tmpR0_gf_; + ParGridFunction tmpR1_gf_; + ParGridFunction vel_gf_; + ParGridFunction gridScale_gf_; + // ParGridFunction *buffer_tInlet_ = nullptr; GridFunctionCoefficient *temperature_bc_field_ = nullptr; @@ -178,6 +192,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { ParBilinearForm *Ht_form_ = nullptr; ParBilinearForm *Mq_form_ = nullptr; ParBilinearForm *LQ_form_ = nullptr; + ParMixedBilinearForm *D_form_ = nullptr; ParLinearForm *LQ_bdry_ = nullptr; OperatorHandle LQ_; @@ -186,6 +201,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { OperatorHandle Ms_; OperatorHandle MsRho_; OperatorHandle Mq_; + OperatorHandle D_op_; mfem::Solver *MsInvPC_ = nullptr; mfem::CGSolver *MsInv_ = nullptr; @@ -199,7 +215,9 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { Vector NTn_, NTnm1_, NTnm2_; Vector Text_; Vector resT_; - Vector tmpR0_, tmpR0b_; + Vector tmpR0_, tmpR0a_, tmpR0b_, tmpR0c_; + Vector tmpR1_; + Vector swDiff_; Vector Qt_; Vector rn_; @@ -233,7 +251,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { public: CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &timeCoeff, - TPS::Tps *tps); + ParGridFunction *gridScale, TPS::Tps *tps); virtual ~CaloricallyPerfectThermoChem(); // Functions overriden from base class @@ -257,6 +275,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { void computeExplicitTempConvectionOP(bool extrap); void computeQt(); void computeQtTO(); + void streamwiseDiffusion(Vector &phi, Vector &swDiff); /// Return a pointer to the current temperature ParGridFunction. ParGridFunction *GetCurrentTemperature() { return &Tn_gf_; } diff --git a/src/loMach.cpp b/src/loMach.cpp index e75f64987..80da539d4 100644 --- a/src/loMach.cpp +++ b/src/loMach.cpp @@ -165,7 +165,8 @@ void LoMachSolver::initialize() { if (loMach_opts_.thermo_solver == "constant-property") { thermo_ = new ConstantPropertyThermoChem(pmesh_, loMach_opts_.order, tpsP_); } else if (loMach_opts_.thermo_solver == "calorically-perfect") { - thermo_ = new CaloricallyPerfectThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_); + thermo_ = + new CaloricallyPerfectThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_); } else if (loMach_opts_.thermo_solver == "lte-thermo-chem") { thermo_ = new LteThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_); } else if (loMach_opts_.thermo_solver == "reacting-flow") { @@ -184,7 +185,8 @@ void LoMachSolver::initialize() { flow_ = new ZeroFlow(pmesh_, 1); } else if (loMach_opts_.flow_solver == "tomboulides") { // Tomboulides flow solver - flow_ = new Tomboulides(pmesh_, loMach_opts_.order, loMach_opts_.order, temporal_coeff_, tpsP_); + flow_ = new Tomboulides(pmesh_, loMach_opts_.order, loMach_opts_.order, temporal_coeff_, + (meshData_->getGridScale()), tpsP_); } else { // Unknown choice... die if (rank0_) { diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index 0932a1889..8020cef28 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -67,7 +67,8 @@ void Orthogonalize(Vector &v, const ParFiniteElementSpace *pfes) { v -= global_sum / static_cast(global_size); } -Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, TPS::Tps *tps) +Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, + ParGridFunction *gridScale, TPS::Tps *tps) : gll_rules(0, Quadrature1D::GaussLobatto), tpsP_(tps), pmesh_(pmesh), @@ -80,6 +81,7 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS rank0_ = (pmesh_->GetMyRank() == 0); axisym_ = false; nvel_ = dim_; + gridScale_gf_ = gridScale; // make sure there is room for BC attributes if (!(pmesh_->bdr_attributes.Size() == 0)) { @@ -129,6 +131,10 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS tps->getInput("loMach/tomboulides/hsolve-maxIters", hsolve_max_iter_, default_max_iter_); tps->getInput("loMach/tomboulides/msolve-maxIters", mass_inverse_max_iter_, default_max_iter_); } + + tps->getInput("loMach/tomboulides/streamwise-stabilization", sw_stab_, false); + tps->getInput("loMach/tomboulides/Reh_offset", re_offset_, 1.0); + tps->getInput("loMach/tomboulides/Reh_factor", re_factor_, 0.01); } Tomboulides::~Tomboulides() { @@ -224,6 +230,8 @@ Tomboulides::~Tomboulides() { delete gradU_gf_; delete gradV_gf_; delete gradW_gf_; + delete tmpR0_gf_; + delete tmpR1_gf_; delete vfes_; delete vfec_; delete sfes_; @@ -256,6 +264,9 @@ void Tomboulides::initializeSelf() { pp_div_rad_comp_gf_ = new ParGridFunction(pfes_, *pp_div_gf_); u_next_rad_comp_gf_ = new ParGridFunction(pfes_, *u_next_gf_); + tmpR0_gf_ = new ParGridFunction(sfes_); + tmpR1_gf_ = new ParGridFunction(vfes_); + if (axisym_) { utheta_gf_ = new ParGridFunction(pfes_); utheta_next_gf_ = new ParGridFunction(pfes_); @@ -333,7 +344,10 @@ void Tomboulides::initializeSelf() { utheta_next_vec_.SetSize(pfes_truevsize); } + swDiff_vec_.SetSize(vfes_truevsize); tmpR0_.SetSize(sfes_truevsize); + tmpR0a_.SetSize(sfes_truevsize); + tmpR0b_.SetSize(sfes_truevsize); tmpR1_.SetSize(vfes_truevsize); gradU_.SetSize(vfes_truevsize); @@ -1345,6 +1359,17 @@ void Tomboulides::step() { }); } + // Add streamwise stability to rhs + if (sw_stab_) { + for (int i = 0; i < dim_; i++) { + setScalarFromVector(u_vec_, i, &tmpR0a_); + streamwiseDiffusion(tmpR0a_, tmpR0b_); + setVectorFromScalar(tmpR0b_, i, &swDiff_vec_); + } + Mv_rho_inv_->Mult(swDiff_vec_, tmpR1_); + pp_div_vec_ += tmpR1_; + } + // Add ustar/dt contribution pp_div_vec_ += ustar_vec_; @@ -1436,6 +1461,9 @@ void Tomboulides::step() { // rho * vstar / dt term Mv_rho_op_->AddMult(ustar_vec_, resu_vec_); + // Add streamwise diffusion + if (sw_stab_) resu_vec_ += swDiff_vec_; + for (auto &vel_dbc : vel_dbcs_) { u_next_gf_->ProjectBdrCoefficient(*vel_dbc.coeff, vel_dbc.attr); } @@ -1686,3 +1714,43 @@ void Tomboulides::evaluateVelocityGradient() { gradV_gf_->SetFromTrueDofs(gradV_); gradW_gf_->SetFromTrueDofs(gradW_); } + +// f(Re_h) * Mu_sw * div(streamwiseGrad), i.e. this does not consider grad of f(Re_h) * Mu_sw +void Tomboulides::streamwiseDiffusion(Vector &phi, Vector &swDiff) { + // compute streamwise gradient of input field + tmpR0_gf_->SetFromTrueDofs(phi); + streamwiseGrad(*tmpR0_gf_, *u_curr_gf_, *tmpR1_gf_); + + // divergence of sw-grad + tmpR1_gf_->GetTrueDofs(tmpR1_); + D_form_->Mult(tmpR1_, swDiff); + + (thermo_interface_->density)->GetTrueDofs(rho_vec_); + gridScale_gf_->GetTrueDofs(tmpR0_); + mu_total_gf_->GetTrueDofs(mu_vec_); + + const double *rho = rho_vec_.HostRead(); + const double *del = tmpR0_.HostRead(); + const double *vel = u_vec_.HostRead(); + const double *mu = mu_vec_.HostRead(); + double *data = swDiff.HostReadWrite(); + + int Sdof = rho_vec_.Size(); + for (int dof = 0; dof < Sdof; dof++) { + double Umag = 0.0; + for (int i = 0; i < dim_; i++) Umag += vel[i] * vel[i]; + Umag = std::sqrt(Umag); + + // element Re + double Re = Umag * del[dof] * rho[dof] / mu[dof]; + + // SUPG weight + double Csupg = 0.5 * (tanh(re_factor_ * Re - re_offset_) + 1.0); + + // streamwise diffusion coeff + double CswDiff = Csupg * Umag * del[dof] * rho[dof]; + + // scaled streamwise Laplacian + data[dof] *= CswDiff; + } +} diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp index 8bc5d8584..7b408d3fe 100644 --- a/src/tomboulides.hpp +++ b/src/tomboulides.hpp @@ -163,6 +163,10 @@ class Tomboulides final : public FlowBase { double hsolve_rtol_; double hsolve_atol_; + bool sw_stab_; + double re_offset_; + double re_factor_; + // To use "numerical integration", quadrature rule must persist mfem::IntegrationRules gll_rules; @@ -225,6 +229,8 @@ class Tomboulides final : public FlowBase { mfem::ParGridFunction *gradW_gf_ = nullptr; // mfem::ParGridFunction *buffer_uInlet_ = nullptr; mfem::VectorGridFunctionCoefficient *velocity_field_ = nullptr; + mfem::ParGridFunction *tmpR0_gf_ = nullptr; + mfem::ParGridFunction *tmpR1_gf_ = nullptr; /// Pressure FEM objects and fields mfem::FiniteElementCollection *pfec_ = nullptr; @@ -238,6 +244,8 @@ class Tomboulides final : public FlowBase { mfem::ParGridFunction *utheta_next_gf_ = nullptr; mfem::ParGridFunction *u_next_rad_comp_gf_ = nullptr; + mfem::ParGridFunction *gridScale_gf_ = nullptr; + /// "total" viscosity, including fluid, turbulence, sponge mfem::ParGridFunction *mu_total_gf_ = nullptr; @@ -355,7 +363,10 @@ class Tomboulides final : public FlowBase { mfem::Vector resp_vec_; mfem::Vector p_vec_; mfem::Vector resu_vec_; + mfem::Vector swDiff_vec_; mfem::Vector tmpR0_; + mfem::Vector tmpR0a_; + mfem::Vector tmpR0b_; mfem::Vector tmpR1_; mfem::Vector gradU_; mfem::Vector gradV_; @@ -391,7 +402,8 @@ class Tomboulides final : public FlowBase { public: /// Constructor - Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, TPS::Tps *tps = nullptr); + Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, + ParGridFunction *gridScale, TPS::Tps *tps = nullptr); /// Destructor ~Tomboulides() final; @@ -429,6 +441,12 @@ class Tomboulides final : public FlowBase { */ void initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) const final; + /** + * @brief Computes f(Re_h) * |U|*h * div(M_sw*grad(phi)) where M_sw transforms the gradient into the + * streamwise direction + */ + void streamwiseDiffusion(Vector &phi, Vector &swDiff); + /// Advance void step() final; diff --git a/src/utils.cpp b/src/utils.cpp index 8dc7ab759..e8663e7f8 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1163,6 +1163,103 @@ bool copyFile(const char *SRC, const char *DEST) { return src && dest; } +void streamwiseGrad(ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad) { + // need to make this general... + int dim_ = 3; + + // compute gradient of input field + scalarGrad3D(phi, swGrad); + + const double *vel = u.HostRead(); + double *gPhi = swGrad.HostReadWrite(); + + int Sdof = phi.Size(); + for (int dof = 0; dof < Sdof; dof++) { + // streamwise coordinate system + Vector unitNorm; + Vector unitT1; + Vector unitT2; + unitNorm.SetSize(dim_); + unitT1.SetSize(dim_); + unitT2.SetSize(dim_); + + // streamwise direction + for (int i = 0; i < dim_; i++) unitNorm[i] = vel[dof + i * Sdof]; + double mod = 0.0; + for (int i = 0; i < dim_; i++) mod += unitNorm[i] * unitNorm[i]; + unitNorm /= mod; + double Umag = std::sqrt(mod); + + // tangent direction (not unique) + int maxInd, minusInd, plusInd; + if (unitNorm[0] > unitNorm[1]) { + maxInd = 0; + } else { + maxInd = 1; + } + if (dim_ == 3 && unitNorm[maxInd] < unitNorm[2]) { + maxInd = 2; + } + minusInd = (maxInd - 1) % dim_; + plusInd = (maxInd + 1) % dim_; + unitT1[minusInd] = -unitNorm[maxInd]; + unitT1[maxInd] = unitNorm[minusInd]; + unitT1[plusInd] = 0.0; + mod = 0.0; + for (int i = 0; i < dim_; i++) mod += unitT1[i] * unitT1[i]; + unitT1 /= mod; + + // t2 is then orthogonal to both normal & t1 + if (dim_ == 3) { + unitT2[0] = +(unitNorm[1] * unitT1[2] - unitNorm[2] * unitT1[1]); + unitT2[1] = -(unitNorm[0] * unitT1[2] - unitNorm[2] * unitT1[0]); + unitT2[2] = +(unitNorm[0] * unitT1[1] - unitNorm[1] * unitT1[0]); + } + + // transform from streamwise coords to global + DenseMatrix M(dim_, dim_); + for (int d = 0; d < dim_; d++) { + M(d, 0) = unitNorm[d]; + M(d, 1) = unitT1[d]; + if (dim_ == 3) M(d, 2) = unitT2[d]; + } + + // streamwise diffusion coeff + DenseMatrix swM(dim_, dim_); + swM = 0.0; + swM(0, 0) = 1.0; + + // muSWgbl = M_{im} muSw_{mn} M_{jn} or M*mu*M^T (with n,t1,t2 in columns of M) + DenseMatrix swMgbl(dim_, dim_); + swMgbl = 0.0; + for (int i = 0; i < dim_; i++) { + for (int j = 0; j < dim_; j++) { + for (int m = 0; m < dim_; m++) { + for (int n = 0; n < dim_; n++) { + swMgbl(i, j) += M(i, m) * M(j, n) * swM(m, n); + } + } + } + } + + // mu*gPhi + Vector tmp1; + tmp1.SetSize(dim_); + for (int i = 0; i < dim_; i++) tmp1[i] = gPhi[dof + i * Sdof]; + Vector tmp2; + tmp2.SetSize(dim_); + tmp2 = 0.0; + for (int i = 0; i < dim_; i++) { + for (int j = 0; j < dim_; j++) { + tmp2[i] += swMgbl(i, j) * tmp1[j]; + } + } + + // copy back to input vector gf + for (int i = 0; i < dim_; i++) gPhi[dof + i * Sdof] = tmp2[i]; + } +} + namespace mfem { GradientVectorGridFunctionCoefficient::GradientVectorGridFunctionCoefficient(const GridFunction *gf) : MatrixCoefficient((gf) ? gf->VectorDim() : 0) { diff --git a/src/utils.hpp b/src/utils.hpp index 2cd4f0d86..5174787b6 100644 --- a/src/utils.hpp +++ b/src/utils.hpp @@ -173,8 +173,8 @@ void vectorGrad3D(ParGridFunction &u, ParGridFunction &gu, ParGridFunction &gv, void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu); void vectorGrad3DV(FiniteElementSpace *fes, Vector u, Vector *gu, Vector *gv, Vector *gw); void scalarGrad3DV(FiniteElementSpace *fes, FiniteElementSpace *vfes, Vector u, Vector *gu); - bool copyFile(const char *SRC, const char *DEST); +void streamwiseGrad(ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad); /// Eliminate essential BCs in an Operator and apply to RHS. /// rename this to something sensible "ApplyEssentialBC" or something From f412352ca3045129eba21047db47e8612dafc414 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Mon, 1 Jul 2024 13:40:57 -0700 Subject: [PATCH 02/12] builds and runs, more evaluation necessary --- src/calorically_perfect.cpp | 14 +++--- src/calorically_perfect.hpp | 2 +- src/tomboulides.cpp | 4 +- src/utils.cpp | 98 ++++++++++++++++++++++++------------- src/utils.hpp | 2 +- 5 files changed, 74 insertions(+), 46 deletions(-) diff --git a/src/calorically_perfect.cpp b/src/calorically_perfect.cpp index d07324157..e6c49c3d9 100644 --- a/src/calorically_perfect.cpp +++ b/src/calorically_perfect.cpp @@ -63,7 +63,7 @@ CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, : tpsP_(tps), pmesh_(pmesh), dim_(pmesh->Dimension()), time_coeff_(time_coeff) { rank0_ = (pmesh_->GetMyRank() == 0); order_ = loMach_opts->order; - gridScale_gf_ = *gridScale; + gridScale_gf_ = gridScale; std::string visc_model; tpsP_->getInput("loMach/calperfect/viscosity-model", visc_model, std::string("sutherland")); @@ -135,9 +135,9 @@ CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, tpsP_->getInput("loMach/calperfect/msolve-max-iter", mass_inverse_max_iter_, max_iter_); tpsP_->getInput("loMach/calperfect/msolve-verbosity", mass_inverse_pl_, pl_solve_); - tps->getInput("loMach/tomboulides/streamwise-stabilization", sw_stab_, false); - tps->getInput("loMach/tomboulides/Reh_offset", re_offset_, 1.0); - tps->getInput("loMach/tomboulides/Reh_factor", re_factor_, 0.01); + tps->getInput("loMach/calperfect/streamwise-stabilization", sw_stab_, false); + tps->getInput("loMach/calperfect/Reh_offset", re_offset_, 1.0); + tps->getInput("loMach/calperfect/Reh_factor", re_factor_, 0.1); } CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() { @@ -1073,13 +1073,13 @@ void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDi tmpR0_gf_.SetFromTrueDofs(phi); (flow_interface_->velocity)->GetTrueDofs(tmpR0a_); vel_gf_.SetFromTrueDofs(tmpR0a_); - streamwiseGrad(tmpR0_gf_, vel_gf_, tmpR1_gf_); + streamwiseGrad(dim_, tmpR0_gf_, vel_gf_, tmpR1_gf_); // divergence of sw-grad tmpR1_gf_.GetTrueDofs(tmpR1_); - D_form_->Mult(tmpR1_, swDiff); + D_op_->Mult(tmpR1_, swDiff); - gridScale_gf_.GetTrueDofs(tmpR0b_); + gridScale_gf_->GetTrueDofs(tmpR0b_); (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); const double *rho = rn_.HostRead(); diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp index 3851b2881..7f276dbd8 100644 --- a/src/calorically_perfect.hpp +++ b/src/calorically_perfect.hpp @@ -167,7 +167,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { ParGridFunction tmpR0_gf_; ParGridFunction tmpR1_gf_; ParGridFunction vel_gf_; - ParGridFunction gridScale_gf_; + ParGridFunction *gridScale_gf_ = nullptr; // ParGridFunction *buffer_tInlet_ = nullptr; GridFunctionCoefficient *temperature_bc_field_ = nullptr; diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index 8020cef28..4dc42b906 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -1719,11 +1719,11 @@ void Tomboulides::evaluateVelocityGradient() { void Tomboulides::streamwiseDiffusion(Vector &phi, Vector &swDiff) { // compute streamwise gradient of input field tmpR0_gf_->SetFromTrueDofs(phi); - streamwiseGrad(*tmpR0_gf_, *u_curr_gf_, *tmpR1_gf_); + streamwiseGrad(dim_, *tmpR0_gf_, *u_curr_gf_, *tmpR1_gf_); // divergence of sw-grad tmpR1_gf_->GetTrueDofs(tmpR1_); - D_form_->Mult(tmpR1_, swDiff); + D_op_->Mult(tmpR1_, swDiff); (thermo_interface_->density)->GetTrueDofs(rho_vec_); gridScale_gf_->GetTrueDofs(tmpR0_); diff --git a/src/utils.cpp b/src/utils.cpp index e8663e7f8..3fc5ff147 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -869,7 +869,7 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { vals1(dof) = grad(0, 0); vals2(dof) = grad(0, 1); - vals3(dof) = grad(0, 2); + if (dim_ == 3) vals3(dof) = grad(0, 2); } // Accumulate values in all dofs, count the zones. @@ -881,9 +881,11 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { int ldof = vdofs[j]; gu(ldof + 1 * nSize) += vals2[j]; } - for (int j = 0; j < vdofs.Size(); j++) { - int ldof = vdofs[j]; - gu(ldof + 2 * nSize) += vals3[j]; + if (dim_ == 3) { + for (int j = 0; j < vdofs.Size(); j++) { + int ldof = vdofs[j]; + gu(ldof + 2 * nSize) += vals3[j]; + } } for (int j = 0; j < vdofs.Size(); j++) { @@ -1163,9 +1165,16 @@ bool copyFile(const char *SRC, const char *DEST) { return src && dest; } -void streamwiseGrad(ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad) { +void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad) { // need to make this general... - int dim_ = 3; + // int dim_ = 3; + + /* + std::cout << "maxInd minusInd plusInd" << endl; + std::cout << 0 << " " << ((0-1) % dim_ + dim_) % dim_ << " " << (0 + 1) % dim_ << endl; + std::cout << 1 << " " << ((1-1) % dim_ + dim_) % dim_ << " " << (1 + 1) % dim_ << endl; + std::cout << 2 << " " << ((2-1) % dim_ + dim_) % dim_ << " " << (2 + 1) % dim_ << endl; + */ // compute gradient of input field scalarGrad3D(phi, swGrad); @@ -1179,63 +1188,81 @@ void streamwiseGrad(ParGridFunction &phi, ParGridFunction &u, ParGridFunction &s Vector unitNorm; Vector unitT1; Vector unitT2; - unitNorm.SetSize(dim_); - unitT1.SetSize(dim_); - unitT2.SetSize(dim_); + unitNorm.SetSize(dim); + unitT1.SetSize(dim); + unitT2.SetSize(dim); // streamwise direction - for (int i = 0; i < dim_; i++) unitNorm[i] = vel[dof + i * Sdof]; + for (int i = 0; i < dim; i++) unitNorm[i] = vel[dof + i * Sdof]; double mod = 0.0; - for (int i = 0; i < dim_; i++) mod += unitNorm[i] * unitNorm[i]; - unitNorm /= mod; + for (int i = 0; i < dim; i++) mod += unitNorm[i] * unitNorm[i]; + mod = std::max(mod, 1.0e-18); double Umag = std::sqrt(mod); + unitNorm /= Umag; + + // check for zero-flow areas + if (Umag < 1.0e-8) { + for (int i = 0; i < dim; i++) gPhi[dof + i * Sdof] = 0.0; + continue; + } // tangent direction (not unique) int maxInd, minusInd, plusInd; - if (unitNorm[0] > unitNorm[1]) { + if (unitNorm[0] * unitNorm[0] > unitNorm[1] * unitNorm[1]) { maxInd = 0; } else { maxInd = 1; } - if (dim_ == 3 && unitNorm[maxInd] < unitNorm[2]) { + if (dim == 3 && unitNorm[maxInd] * unitNorm[maxInd] < unitNorm[2] * unitNorm[2]) { maxInd = 2; } - minusInd = (maxInd - 1) % dim_; - plusInd = (maxInd + 1) % dim_; + minusInd = ((maxInd - 1) % dim + dim) % dim; + plusInd = (maxInd + 1) % dim; + unitT1[minusInd] = -unitNorm[maxInd]; unitT1[maxInd] = unitNorm[minusInd]; unitT1[plusInd] = 0.0; mod = 0.0; - for (int i = 0; i < dim_; i++) mod += unitT1[i] * unitT1[i]; - unitT1 /= mod; + for (int i = 0; i < dim; i++) mod += unitT1[i] * unitT1[i]; + unitT1 /= std::sqrt(mod); // t2 is then orthogonal to both normal & t1 - if (dim_ == 3) { + if (dim == 3) { unitT2[0] = +(unitNorm[1] * unitT1[2] - unitNorm[2] * unitT1[1]); unitT2[1] = -(unitNorm[0] * unitT1[2] - unitNorm[2] * unitT1[0]); unitT2[2] = +(unitNorm[0] * unitT1[1] - unitNorm[1] * unitT1[0]); } // transform from streamwise coords to global - DenseMatrix M(dim_, dim_); - for (int d = 0; d < dim_; d++) { + DenseMatrix M(dim, dim); + for (int d = 0; d < dim; d++) { M(d, 0) = unitNorm[d]; M(d, 1) = unitT1[d]; - if (dim_ == 3) M(d, 2) = unitT2[d]; + if (dim == 3) M(d, 2) = unitT2[d]; } // streamwise diffusion coeff - DenseMatrix swM(dim_, dim_); + DenseMatrix swM(dim, dim); swM = 0.0; swM(0, 0) = 1.0; - // muSWgbl = M_{im} muSw_{mn} M_{jn} or M*mu*M^T (with n,t1,t2 in columns of M) - DenseMatrix swMgbl(dim_, dim_); - swMgbl = 0.0; + /* + std::cout << " " << endl; for (int i = 0; i < dim_; i++) { for (int j = 0; j < dim_; j++) { - for (int m = 0; m < dim_; m++) { - for (int n = 0; n < dim_; n++) { + std::cout << M(i,j) << " " ; + } + std::cout << endl; + } + */ + + // muSWgbl = M_{im} muSw_{mn} M_{jn} or M*mu*M^T (with n,t1,t2 in columns of M) + DenseMatrix swMgbl(dim, dim); + swMgbl = 0.0; + for (int i = 0; i < dim; i++) { + for (int j = 0; j < dim; j++) { + for (int m = 0; m < dim; m++) { + for (int n = 0; n < dim; n++) { swMgbl(i, j) += M(i, m) * M(j, n) * swM(m, n); } } @@ -1244,19 +1271,20 @@ void streamwiseGrad(ParGridFunction &phi, ParGridFunction &u, ParGridFunction &s // mu*gPhi Vector tmp1; - tmp1.SetSize(dim_); - for (int i = 0; i < dim_; i++) tmp1[i] = gPhi[dof + i * Sdof]; + tmp1.SetSize(dim); + for (int i = 0; i < dim; i++) tmp1[i] = gPhi[dof + i * Sdof]; + Vector tmp2; - tmp2.SetSize(dim_); - tmp2 = 0.0; - for (int i = 0; i < dim_; i++) { - for (int j = 0; j < dim_; j++) { + tmp2.SetSize(dim); + for (int i = 0; i < dim; i++) tmp2[i] = 0.0; + for (int i = 0; i < dim; i++) { + for (int j = 0; j < dim; j++) { tmp2[i] += swMgbl(i, j) * tmp1[j]; } } // copy back to input vector gf - for (int i = 0; i < dim_; i++) gPhi[dof + i * Sdof] = tmp2[i]; + for (int i = 0; i < dim; i++) gPhi[dof + i * Sdof] = tmp2[i]; } } diff --git a/src/utils.hpp b/src/utils.hpp index 5174787b6..603e0804b 100644 --- a/src/utils.hpp +++ b/src/utils.hpp @@ -174,7 +174,7 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu); void vectorGrad3DV(FiniteElementSpace *fes, Vector u, Vector *gu, Vector *gv, Vector *gw); void scalarGrad3DV(FiniteElementSpace *fes, FiniteElementSpace *vfes, Vector u, Vector *gu); bool copyFile(const char *SRC, const char *DEST); -void streamwiseGrad(ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad); +void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad); /// Eliminate essential BCs in an Operator and apply to RHS. /// rename this to something sensible "ApplyEssentialBC" or something From a931a3414078a408514cb61a8f9ed262f5cc0043 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Tue, 2 Jul 2024 07:09:52 -0700 Subject: [PATCH 03/12] adding lte support --- src/calorically_perfect.hpp | 1 + src/loMach.cpp | 2 +- src/lte_thermo_chem.cpp | 89 ++++++++++++++++++++++++++++++++++++- src/lte_thermo_chem.hpp | 26 ++++++++++- src/utils.cpp | 5 +-- 5 files changed, 114 insertions(+), 9 deletions(-) diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp index 7f276dbd8..769430336 100644 --- a/src/calorically_perfect.hpp +++ b/src/calorically_perfect.hpp @@ -111,6 +111,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { double hsolve_rtol_; double hsolve_atol_; + // streamwise-stabilization bool sw_stab_; double re_offset_; double re_factor_; diff --git a/src/loMach.cpp b/src/loMach.cpp index 80da539d4..d37828bd0 100644 --- a/src/loMach.cpp +++ b/src/loMach.cpp @@ -168,7 +168,7 @@ void LoMachSolver::initialize() { thermo_ = new CaloricallyPerfectThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_); } else if (loMach_opts_.thermo_solver == "lte-thermo-chem") { - thermo_ = new LteThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_); + thermo_ = new LteThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_); } else if (loMach_opts_.thermo_solver == "reacting-flow") { thermo_ = new ReactingFlow(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_); } else { diff --git a/src/lte_thermo_chem.cpp b/src/lte_thermo_chem.cpp index c41a48ebb..d057c19b1 100644 --- a/src/lte_thermo_chem.cpp +++ b/src/lte_thermo_chem.cpp @@ -68,10 +68,11 @@ static double sigmaTorchStartUp(const Vector &pos) { static FunctionCoefficient sigma_start_up(sigmaTorchStartUp); LteThermoChem::LteThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &time_coeff, - TPS::Tps *tps) - : tpsP_(tps), pmesh_(pmesh), time_coeff_(time_coeff) { + ParGridFunction *gridScale, TPS::Tps *tps) + : tpsP_(tps), pmesh_(pmesh), dim_(pmesh->Dimension()), time_coeff_(time_coeff) { rank0_ = (pmesh_->GetMyRank() == 0); order_ = loMach_opts->order; + gridScale_gf_ = gridScale; tps->getInput("loMach/axisymmetric", axisym_, false); @@ -166,6 +167,10 @@ LteThermoChem::LteThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, t tps->getInput("loMach/ltethermo/linear-solver-rtol", rtol_, 1e-12); tps->getInput("loMach/ltethermo/linear-solver-max-iter", max_iter_, 1000); tps->getInput("loMach/ltethermo/linear-solver-verbosity", pl_solve_, 0); + + tps->getInput("loMach/ltethermo/streamwise-stabilization", sw_stab_, false); + tps->getInput("loMach/ltethermo/Reh_offset", re_offset_, 1.0); + tps->getInput("loMach/ltethermo/Reh_factor", re_factor_, 0.1); } LteThermoChem::~LteThermoChem() { @@ -192,6 +197,7 @@ LteThermoChem::~LteThermoChem() { delete M_rho_Cp_form_; delete Ms_form_; delete At_form_; + delete D_form_; delete rho_Cp_u_coeff_; delete un_next_coeff_; delete kap_gradT_coeff_; @@ -208,6 +214,8 @@ LteThermoChem::~LteThermoChem() { // allocated in initializeSelf delete sfes_; delete sfec_; + delete vfes_; + delete vfec_; } void LteThermoChem::initializeSelf() { @@ -219,6 +227,9 @@ void LteThermoChem::initializeSelf() { sfec_ = new H1_FECollection(order_); sfes_ = new ParFiniteElementSpace(pmesh_, sfec_); + vfec_ = new H1_FECollection(order_, dim_); + vfes_ = new ParFiniteElementSpace(pmesh_, vfec_, dim_); + // Check if fully periodic mesh if (!(pmesh_->bdr_attributes.Size() == 0)) { temp_ess_attr_.SetSize(pmesh_->bdr_attributes.Max()); @@ -230,6 +241,7 @@ void LteThermoChem::initializeSelf() { if (rank0_) grvy_printf(ginfo, "LteThermoChem paces constructed...\n"); int sfes_truevsize = sfes_->GetTrueVSize(); + int vfes_truevsize = vfes_->GetTrueVSize(); Qt_.SetSize(sfes_truevsize); Qt_ = 0.0; @@ -317,8 +329,16 @@ void LteThermoChem::initializeSelf() { radiation_sink_.SetSize(sfes_truevsize); radiation_sink_ = 0.0; + vel_gf_.SetSpace(vfes_); + tmpR1_gf_.SetSpace(vfes_); + tmpR0_gf_.SetSpace(sfes_); + + swDiff_.SetSize(sfes_truevsize); tmpR0_.SetSize(sfes_truevsize); + tmpR0a_.SetSize(sfes_truevsize); tmpR0b_.SetSize(sfes_truevsize); + tmpR0c_.SetSize(sfes_truevsize); + tmpR1_.SetSize(vfes_truevsize); R0PM0_gf_.SetSpace(sfes_); @@ -628,6 +648,24 @@ void LteThermoChem::initializeOperators() { M_rho_form_->Assemble(); M_rho_form_->FormSystemMatrix(empty, M_rho_); + // Divergence operator + D_form_ = new ParMixedBilinearForm(vfes_, sfes_); + VectorDivergenceIntegrator *vd_mblfi; + // if (axisym_) { + // vd_mblfi = new VectorDivergenceIntegrator(radius_coeff); + // } else { + vd_mblfi = new VectorDivergenceIntegrator(); + // } + if (numerical_integ_) { + vd_mblfi->SetIntRule(&ir_nli); + } + D_form_->AddDomainIntegrator(vd_mblfi); + if (partial_assembly_) { + D_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL); + } + D_form_->Assemble(); + D_form_->FormRectangularSystemMatrix(empty, empty, D_op_); + // helmholtz Ht_form_ = new ParBilinearForm(sfes_); MassIntegrator *hmt_blfi; @@ -905,6 +943,12 @@ void LteThermoChem::step() { tmpR0_ = dtP_; Ms_->AddMult(tmpR0_, resT_); + // Add streamwise stability to rhs + if (sw_stab_) { + streamwiseDiffusion(Tn_, swDiff_); + resT_.Add(1.0, swDiff_); + } + // Joule heating (and radiation sink) jh_form_->Update(); jh_form_->Assemble(); @@ -1231,3 +1275,44 @@ void LteThermoChem::computeQt() { Qt_.Neg(); Qt_gf_.SetFromTrueDofs(Qt_); } + +void LteThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { + // compute streamwise gradient of input field + tmpR0_gf_.SetFromTrueDofs(phi); + (flow_interface_->velocity)->GetTrueDofs(tmpR0a_); + vel_gf_.SetFromTrueDofs(tmpR0a_); + streamwiseGrad(dim_, tmpR0_gf_, vel_gf_, tmpR1_gf_); + + // divergence of sw-grad + tmpR1_gf_.GetTrueDofs(tmpR1_); + D_op_->Mult(tmpR1_, swDiff); + + gridScale_gf_->GetTrueDofs(tmpR0b_); + (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + + const double *rho = rn_.HostRead(); + const double *vel = tmpR0a_.HostRead(); + const double *del = tmpR0b_.HostRead(); + const double *mu = visc_.HostRead(); + const double *muT = tmpR0c_.HostRead(); + double *data = swDiff.HostReadWrite(); + + int Sdof = rn_.Size(); + for (int dof = 0; dof < Sdof; dof++) { + double Umag = 0.0; + for (int i = 0; i < dim_; i++) Umag += vel[i] * vel[i]; + Umag = std::sqrt(Umag); + + // element Re + double Re = Umag * del[dof] * rho[dof] / (mu[dof] + muT[dof]); + + // SUPG weight + double Csupg = 0.5 * (tanh(re_factor_ * Re - re_offset_) + 1.0); + + // streamwise diffusion coeff + double CswDiff = Csupg * Umag * del[dof] * rho[dof]; + + // scaled streamwise Laplacian + data[dof] *= CswDiff; + } +} diff --git a/src/lte_thermo_chem.hpp b/src/lte_thermo_chem.hpp index 42e5114bf..2aaed5e61 100644 --- a/src/lte_thermo_chem.hpp +++ b/src/lte_thermo_chem.hpp @@ -79,6 +79,7 @@ class LteThermoChem final : public ThermoChemModelBase { // Mesh and discretization scheme info ParMesh *pmesh_ = nullptr; + int dim_; int order_; IntegrationRules gll_rules_; const temporalSchemeCoefficients &time_coeff_; @@ -95,6 +96,11 @@ class LteThermoChem final : public ThermoChemModelBase { int max_iter_; /**< Maximum number of linear solver iterations */ double rtol_ = 1e-12; /**< Linear solver relative tolerance */ + // streamwise-stabilization + bool sw_stab_; + double re_offset_; + double re_factor_; + // Boundary condition info Array temp_ess_attr_; /**< List of patches with Dirichlet BC on temperature */ Array Qt_ess_attr_; /**< List of patches with Dirichlet BC on Q (thermal divergence) */ @@ -131,6 +137,10 @@ class LteThermoChem final : public ThermoChemModelBase { // Scalar \f$H^1\f$ finite element space. ParFiniteElementSpace *sfes_ = nullptr; + // Vector fe collection and space + FiniteElementCollection *vfec_ = nullptr; + ParFiniteElementSpace *vfes_ = nullptr; + // Fields ParGridFunction Tnm1_gf_, Tnm2_gf_; ParGridFunction Tn_gf_, Tn_next_gf_, Text_gf_, resT_gf_; @@ -148,6 +158,11 @@ class LteThermoChem final : public ThermoChemModelBase { ParGridFunction R0PM0_gf_; ParGridFunction Qt_gf_; + ParGridFunction tmpR0_gf_; + ParGridFunction tmpR1_gf_; + ParGridFunction vel_gf_; + ParGridFunction *gridScale_gf_ = nullptr; + // ParGridFunction *buffer_tInlet_ = nullptr; GridFunctionCoefficient *temperature_bc_field_ = nullptr; @@ -195,6 +210,8 @@ class LteThermoChem final : public ThermoChemModelBase { ParBilinearForm *LQ_form_ = nullptr; ParLinearForm *LQ_bdry_ = nullptr; + ParMixedBilinearForm *D_form_ = nullptr; + OperatorHandle At_; OperatorHandle Ht_; OperatorHandle Ms_; @@ -203,6 +220,7 @@ class LteThermoChem final : public ThermoChemModelBase { OperatorHandle M_rho_Cp_; OperatorHandle M_rho_; OperatorHandle A_rho_; + OperatorHandle D_op_; mfem::Solver *MsInvPC_ = nullptr; mfem::CGSolver *MsInv_ = nullptr; @@ -218,7 +236,9 @@ class LteThermoChem final : public ThermoChemModelBase { Vector NTn_, NTnm1_, NTnm2_; Vector Text_; Vector resT_; - Vector tmpR0_, tmpR0b_; + Vector tmpR0_, tmpR0a_, tmpR0b_, tmpR0c_; + Vector tmpR1_; + Vector swDiff_; Vector Qt_; Vector rn_, rnm1_, rnm2_, rnm3_; @@ -241,7 +261,8 @@ class LteThermoChem final : public ThermoChemModelBase { ParGridFunction Tn_filtered_gf_; public: - LteThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &timeCoeff, TPS::Tps *tps); + LteThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &timeCoeff, + ParGridFunction *gridScale, TPS::Tps *tps); virtual ~LteThermoChem(); // Functions overriden from base class @@ -261,6 +282,7 @@ class LteThermoChem final : public ThermoChemModelBase { void computeExplicitTempConvectionOP(); void computeQt(); void updateHistory(); + void streamwiseDiffusion(Vector &phi, Vector &swDiff); /// Return a pointer to the current temperature ParGridFunction. ParGridFunction *GetCurrentTemperature() { return &Tn_gf_; } diff --git a/src/utils.cpp b/src/utils.cpp index 3fc5ff147..afb14cd23 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1166,9 +1166,6 @@ bool copyFile(const char *SRC, const char *DEST) { } void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad) { - // need to make this general... - // int dim_ = 3; - /* std::cout << "maxInd minusInd plusInd" << endl; std::cout << 0 << " " << ((0-1) % dim_ + dim_) % dim_ << " " << (0 + 1) % dim_ << endl; @@ -1241,7 +1238,7 @@ void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFu if (dim == 3) M(d, 2) = unitT2[d]; } - // streamwise diffusion coeff + // streamwise coeff DenseMatrix swM(dim, dim); swM = 0.0; swM(0, 0) = 1.0; From ff4cdb9dbd26e493201bfab7280d7393db6bdfb8 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 3 Jul 2024 15:32:10 -0700 Subject: [PATCH 04/12] small fix to scalarGrad3D --- src/utils.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils.cpp b/src/utils.cpp index afb14cd23..5843fa626 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -820,7 +820,7 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { // AccumulateAndCountZones. Array zones_per_vdof; - zones_per_vdof.SetSize(3 * (fes->GetVSize())); + zones_per_vdof.SetSize(fes->GetVSize()); zones_per_vdof = 0; gu = 0.0; @@ -860,7 +860,7 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { // Eval and GetVectorGradientHat. el->CalcDShape(tr->GetIntPoint(), dshape); grad_hat.SetSize(vdim, dim); - DenseMatrix loc_data_mat(loc_data.GetData(), elndofs, 1); + DenseMatrix loc_data_mat(loc_data.GetData(), elndofs, vdim); MultAtB(loc_data_mat, dshape, grad_hat); const DenseMatrix &Jinv = tr->InverseJacobian(); From 6abe5045cc55ced5494758525a31b6136e35f36e Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 3 Jul 2024 22:47:56 -0700 Subject: [PATCH 05/12] attempt to remove duplicate code with streamwise-stab. Added Re_h gf and calcs in flow which is exported to thermo-chem, moved supg-like section to utils, changed interface for streamwiseDiffusion routines to take the gradients directly. Commented out communications portions of curl calcs as they seem to be introducing partition imprinting. Fixed small bug in channel temp ic --- src/calorically_perfect.cpp | 101 ++++++++++++++++++++++++++++++---- src/calorically_perfect.hpp | 10 +++- src/cases.cpp | 2 +- src/lte_thermo_chem.cpp | 104 ++++++++++++++++++++++++++++++++---- src/lte_thermo_chem.hpp | 34 +++++++++++- src/split_flow_base.hpp | 4 ++ src/tomboulides.cpp | 90 +++++++++++++++++++++++++------ src/tomboulides.hpp | 7 +++ src/utils.cpp | 48 +++++++++++++++-- src/utils.hpp | 5 +- 10 files changed, 361 insertions(+), 44 deletions(-) diff --git a/src/calorically_perfect.cpp b/src/calorically_perfect.cpp index e6c49c3d9..b93ba7fd1 100644 --- a/src/calorically_perfect.cpp +++ b/src/calorically_perfect.cpp @@ -153,10 +153,14 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() { delete HtInvPC_; delete MsInv_; delete MsInvPC_; + delete Mv_inv_; + delete Mv_inv_pc_; delete Ht_form_; delete MsRho_form_; delete Ms_form_; + delete Mv_form_; delete At_form_; + delete G_form_; delete D_form_; delete rhou_coeff_; delete rhon_next_coeff_; @@ -262,6 +266,9 @@ void CaloricallyPerfectThermoChem::initializeSelf() { tmpR0c_.SetSize(sfes_truevsize); tmpR1_.SetSize(vfes_truevsize); + gradT_.SetSize(vfes_truevsize); + gradT_ = 0.0; + R0PM0_gf_.SetSpace(sfes_); rhoDt.SetSpace(sfes_); @@ -510,6 +517,43 @@ void CaloricallyPerfectThermoChem::initializeOperators() { D_form_->Assemble(); D_form_->FormRectangularSystemMatrix(empty, empty, D_op_); + // Gradient + G_form_ = new ParMixedBilinearForm(sfes_, vfes_); + // auto *g_mblfi = new GradientIntegrator(); + GradientIntegrator *g_mblfi; + // if (axisym_) { + // g_mblfi = new GradientIntegrator(radius_coeff); + // } else { + g_mblfi = new GradientIntegrator(); + //} + if (numerical_integ_) { + g_mblfi->SetIntRule(&ir_nli); + } + G_form_->AddDomainIntegrator(g_mblfi); + if (partial_assembly_) { + G_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL); + } + G_form_->Assemble(); + G_form_->FormRectangularSystemMatrix(empty, empty, G_op_); + + // Mass matrix for the vector (gradT) + Mv_form_ = new ParBilinearForm(vfes_); + VectorMassIntegrator *mv_blfi; + // if (axisym_) { + // mv_blfi = new VectorMassIntegrator(radius_coeff); + // } else { + mv_blfi = new VectorMassIntegrator(); + //} + if (numerical_integ_) { + mv_blfi->SetIntRule(&ir_nli); + } + Mv_form_->AddDomainIntegrator(mv_blfi); + if (partial_assembly_) { + Mv_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL); + } + Mv_form_->Assemble(); + Mv_form_->FormSystemMatrix(temp_ess_tdof_, Mv_); + // mass matrix with rho MsRho_form_ = new ParBilinearForm(sfes_); auto *msrho_blfi = new MassIntegrator(*rho_coeff_); @@ -522,6 +566,7 @@ void CaloricallyPerfectThermoChem::initializeOperators() { MsRho_form_->FormSystemMatrix(empty, MsRho_); if (rank0_) std::cout << "CaloricallyPerfectThermoChem MsRho operator set" << endl; + // Helmholtz Ht_form_ = new ParBilinearForm(sfes_); auto *hmt_blfi = new MassIntegrator(*rho_over_dt_coeff_); auto *hdt_blfi = new DiffusionIntegrator(*thermal_diff_total_coeff_); @@ -556,6 +601,27 @@ void CaloricallyPerfectThermoChem::initializeOperators() { MsInv_->SetAbsTol(mass_inverse_atol_); MsInv_->SetMaxIter(mass_inverse_max_iter_); + // Inverse (unweighted) mass operator (velocity space) + if (partial_assembly_) { + Vector diag_pa(vfes_->GetTrueVSize()); + Mv_form_->AssembleDiagonal(diag_pa); + Mv_inv_pc_ = new OperatorJacobiSmoother(diag_pa, empty); + } else { + Mv_inv_pc_ = new HypreSmoother(*Mv_.As()); + dynamic_cast(Mv_inv_pc_)->SetType(HypreSmoother::Jacobi, smoother_passes_); + dynamic_cast(Mv_inv_pc_)->SetSOROptions(smoother_relax_weight_, smoother_relax_omega_); + dynamic_cast(Mv_inv_pc_) + ->SetPolyOptions(smoother_poly_order_, smoother_poly_fraction_, smoother_eig_est_); + } + Mv_inv_ = new CGSolver(vfes_->GetComm()); + Mv_inv_->iterative_mode = false; + Mv_inv_->SetOperator(*Mv_); + Mv_inv_->SetPreconditioner(*Mv_inv_pc_); + Mv_inv_->SetPrintLevel(mass_inverse_pl_); + Mv_inv_->SetAbsTol(mass_inverse_atol_); + Mv_inv_->SetRelTol(mass_inverse_rtol_); + Mv_inv_->SetMaxIter(mass_inverse_max_iter_); + HtInvPC_ = new HypreSmoother(*Ht_.As()); dynamic_cast(HtInvPC_)->SetType(HypreSmoother::Jacobi, smoother_passes_); dynamic_cast(HtInvPC_)->SetSOROptions(hsmoother_relax_weight_, hsmoother_relax_omega_); @@ -707,7 +773,12 @@ void CaloricallyPerfectThermoChem::step() { // Add streamwise stability to rhs if (sw_stab_) { - streamwiseDiffusion(Tn_, swDiff_); + // compute temp gradient (only really needed for sw-stab) + G_op_->Mult(Text_, tmpR1_); + Mv_inv_->Mult(tmpR1_, gradT_); + + // streamwiseDiffusion(Tn_, swDiff_); + streamwiseDiffusion(gradT_, swDiff_); resT_.Add(1.0, swDiff_); } @@ -1068,25 +1139,35 @@ void CaloricallyPerfectThermoChem::screenValues(std::vector &values) { } } -void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { - // compute streamwise gradient of input field - tmpR0_gf_.SetFromTrueDofs(phi); +// void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { +void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &gradPhi, Vector &swDiff) { (flow_interface_->velocity)->GetTrueDofs(tmpR0a_); vel_gf_.SetFromTrueDofs(tmpR0a_); - streamwiseGrad(dim_, tmpR0_gf_, vel_gf_, tmpR1_gf_); + + // compute streamwise gradient of input field + // tmpR0_gf_.SetFromTrueDofs(phi); + // streamwiseGrad(dim_, tmpR0_gf_, vel_gf_, tmpR1_gf_); + + tmpR1_gf_.SetFromTrueDofs(gradPhi); + streamwiseGrad(dim_, vel_gf_, tmpR1_gf_); // divergence of sw-grad tmpR1_gf_.GetTrueDofs(tmpR1_); D_op_->Mult(tmpR1_, swDiff); gridScale_gf_->GetTrueDofs(tmpR0b_); - (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + //(turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + (flow_interface_->Reh)->GetTrueDofs(tmpR0c_); + upwindDiff(dim_, re_factor_, re_offset_, tmpR0a_, rn_, tmpR0b_, tmpR0c_, swDiff); + + /* const double *rho = rn_.HostRead(); const double *vel = tmpR0a_.HostRead(); const double *del = tmpR0b_.HostRead(); - const double *mu = visc_.HostRead(); - const double *muT = tmpR0c_.HostRead(); + //const double *mu = visc_.HostRead(); + //const double *muT = tmpR0c_.HostRead(); + const double *Reh = tmpR0c_.HostRead(); double *data = swDiff.HostReadWrite(); int Sdof = rn_.Size(); @@ -1096,7 +1177,8 @@ void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDi Umag = std::sqrt(Umag); // element Re - double Re = Umag * del[dof] * rho[dof] / (mu[dof] + muT[dof]); + //double Re = Umag * del[dof] * rho[dof] / (mu[dof] + muT[dof]); + double Re = Reh[dof]; // SUPG weight double Csupg = 0.5 * (tanh(re_factor_ * Re - re_offset_) + 1.0); @@ -1107,4 +1189,5 @@ void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDi // scaled streamwise Laplacian data[dof] *= CswDiff; } + */ } diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp index 769430336..3f6f72e75 100644 --- a/src/calorically_perfect.hpp +++ b/src/calorically_perfect.hpp @@ -189,11 +189,13 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { // operators and solvers ParBilinearForm *At_form_ = nullptr; ParBilinearForm *Ms_form_ = nullptr; + ParBilinearForm *Mv_form_ = nullptr; ParBilinearForm *MsRho_form_ = nullptr; ParBilinearForm *Ht_form_ = nullptr; ParBilinearForm *Mq_form_ = nullptr; ParBilinearForm *LQ_form_ = nullptr; ParMixedBilinearForm *D_form_ = nullptr; + ParMixedBilinearForm *G_form_ = nullptr; ParLinearForm *LQ_bdry_ = nullptr; OperatorHandle LQ_; @@ -202,7 +204,9 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { OperatorHandle Ms_; OperatorHandle MsRho_; OperatorHandle Mq_; + OperatorHandle Mv_; OperatorHandle D_op_; + OperatorHandle G_op_; mfem::Solver *MsInvPC_ = nullptr; mfem::CGSolver *MsInv_ = nullptr; @@ -210,6 +214,8 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { mfem::CGSolver *MqInv_ = nullptr; mfem::Solver *HtInvPC_ = nullptr; mfem::CGSolver *HtInv_ = nullptr; + mfem::Solver *Mv_inv_pc_ = nullptr; + mfem::CGSolver *Mv_inv_ = nullptr; // Vectors Vector Tn_, Tn_next_, Tnm1_, Tnm2_; @@ -219,6 +225,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { Vector tmpR0_, tmpR0a_, tmpR0b_, tmpR0c_; Vector tmpR1_; Vector swDiff_; + Vector gradT_; Vector Qt_; Vector rn_; @@ -276,7 +283,8 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { void computeExplicitTempConvectionOP(bool extrap); void computeQt(); void computeQtTO(); - void streamwiseDiffusion(Vector &phi, Vector &swDiff); + // void streamwiseDiffusion(Vector &phi, Vector &swDiff); + void streamwiseDiffusion(Vector &gradPhi, Vector &swDiff); /// Return a pointer to the current temperature ParGridFunction. ParGridFunction *GetCurrentTemperature() { return &Tn_gf_; } diff --git a/src/cases.cpp b/src/cases.cpp index 6fecd0af9..89529ef26 100644 --- a/src/cases.cpp +++ b/src/cases.cpp @@ -202,7 +202,7 @@ double temp_channel(const Vector &coords, double t) { double Tlo = 200.0; double y = coords(1); double temp; - temp = Tlo + (y + 0.5) * (Thi - Tlo); + temp = Tlo + 0.5 * (y + 0.1) * (Thi - Tlo); return temp; } diff --git a/src/lte_thermo_chem.cpp b/src/lte_thermo_chem.cpp index d057c19b1..482305692 100644 --- a/src/lte_thermo_chem.cpp +++ b/src/lte_thermo_chem.cpp @@ -168,6 +168,17 @@ LteThermoChem::LteThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, t tps->getInput("loMach/ltethermo/linear-solver-max-iter", max_iter_, 1000); tps->getInput("loMach/ltethermo/linear-solver-verbosity", pl_solve_, 0); + // not deleting above block to maintain backwards-compatability + tpsP_->getInput("loMach/ltethermo/hsolve-rtol", hsolve_rtol_, rtol_); + tpsP_->getInput("loMach/ltethermo/hsolve-atol", hsolve_atol_, default_atol_); + tpsP_->getInput("loMach/ltethermo/hsolve-max-iter", hsolve_max_iter_, max_iter_); + tpsP_->getInput("loMach/ltethermo/hsolve-verbosity", hsolve_pl_, pl_solve_); + + tpsP_->getInput("loMach/ltethermo/msolve-rtol", mass_inverse_rtol_, rtol_); + tpsP_->getInput("loMach/ltethermo/msolve-atol", mass_inverse_atol_, default_atol_); + tpsP_->getInput("loMach/ltethermo/msolve-max-iter", mass_inverse_max_iter_, max_iter_); + tpsP_->getInput("loMach/ltethermo/msolve-verbosity", mass_inverse_pl_, pl_solve_); + tps->getInput("loMach/ltethermo/streamwise-stabilization", sw_stab_, false); tps->getInput("loMach/ltethermo/Reh_offset", re_offset_, 1.0); tps->getInput("loMach/ltethermo/Reh_factor", re_factor_, 0.1); @@ -192,12 +203,16 @@ LteThermoChem::~LteThermoChem() { delete MsInvPC_; delete MrhoInv_; delete MrhoInvPC_; + delete Mv_inv_; + delete Mv_inv_pc_; delete Ht_form_; delete M_rho_form_; delete M_rho_Cp_form_; delete Ms_form_; + delete Mv_form_; delete At_form_; delete D_form_; + delete G_form_; delete rho_Cp_u_coeff_; delete un_next_coeff_; delete kap_gradT_coeff_; @@ -666,6 +681,43 @@ void LteThermoChem::initializeOperators() { D_form_->Assemble(); D_form_->FormRectangularSystemMatrix(empty, empty, D_op_); + // Gradient + G_form_ = new ParMixedBilinearForm(sfes_, vfes_); + // auto *g_mblfi = new GradientIntegrator(); + GradientIntegrator *g_mblfi; + if (axisym_) { + g_mblfi = new GradientIntegrator(radius_coeff); + } else { + g_mblfi = new GradientIntegrator(); + } + if (numerical_integ_) { + g_mblfi->SetIntRule(&ir_nli); + } + G_form_->AddDomainIntegrator(g_mblfi); + if (partial_assembly_) { + G_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL); + } + G_form_->Assemble(); + G_form_->FormRectangularSystemMatrix(empty, empty, G_op_); + + // Mass matrix for the vector (gradT) + Mv_form_ = new ParBilinearForm(vfes_); + VectorMassIntegrator *mv_blfi; + if (axisym_) { + mv_blfi = new VectorMassIntegrator(radius_coeff); + } else { + mv_blfi = new VectorMassIntegrator(); + } + if (numerical_integ_) { + mv_blfi->SetIntRule(&ir_nli); + } + Mv_form_->AddDomainIntegrator(mv_blfi); + if (partial_assembly_) { + Mv_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL); + } + Mv_form_->Assemble(); + Mv_form_->FormSystemMatrix(temp_ess_tdof_, Mv_); + // helmholtz Ht_form_ = new ParBilinearForm(sfes_); MassIntegrator *hmt_blfi; @@ -722,6 +774,27 @@ void LteThermoChem::initializeOperators() { MrhoInv_->SetRelTol(rtol_); MrhoInv_->SetMaxIter(max_iter_); + // Inverse (unweighted) mass operator (velocity space) + if (partial_assembly_) { + Vector diag_pa(vfes_->GetTrueVSize()); + Mv_form_->AssembleDiagonal(diag_pa); + Mv_inv_pc_ = new OperatorJacobiSmoother(diag_pa, empty); + } else { + Mv_inv_pc_ = new HypreSmoother(*Mv_.As()); + dynamic_cast(Mv_inv_pc_)->SetType(HypreSmoother::Jacobi, smoother_passes_); + dynamic_cast(Mv_inv_pc_)->SetSOROptions(smoother_relax_weight_, smoother_relax_omega_); + dynamic_cast(Mv_inv_pc_) + ->SetPolyOptions(smoother_poly_order_, smoother_poly_fraction_, smoother_eig_est_); + } + Mv_inv_ = new CGSolver(vfes_->GetComm()); + Mv_inv_->iterative_mode = false; + Mv_inv_->SetOperator(*Mv_); + Mv_inv_->SetPreconditioner(*Mv_inv_pc_); + Mv_inv_->SetPrintLevel(mass_inverse_pl_); + Mv_inv_->SetAbsTol(mass_inverse_atol_); + Mv_inv_->SetRelTol(mass_inverse_rtol_); + Mv_inv_->SetMaxIter(mass_inverse_max_iter_); + if (partial_assembly_) { Vector diag_pa(sfes_->GetTrueVSize()); Ht_form_->AssembleDiagonal(diag_pa); @@ -945,7 +1018,12 @@ void LteThermoChem::step() { // Add streamwise stability to rhs if (sw_stab_) { - streamwiseDiffusion(Tn_, swDiff_); + // compute temp gradient (only really needed for sw-stab) + G_op_->Mult(Text_, tmpR1_); + Mv_inv_->Mult(tmpR1_, gradT_); + + // streamwiseDiffusion(Tn_, swDiff_); + streamwiseDiffusion(gradT_, swDiff_); resT_.Add(1.0, swDiff_); } @@ -1276,25 +1354,32 @@ void LteThermoChem::computeQt() { Qt_gf_.SetFromTrueDofs(Qt_); } -void LteThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { - // compute streamwise gradient of input field - tmpR0_gf_.SetFromTrueDofs(phi); +// void LteThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { +void LteThermoChem::streamwiseDiffusion(Vector &gradPhi, Vector &swDiff) { (flow_interface_->velocity)->GetTrueDofs(tmpR0a_); vel_gf_.SetFromTrueDofs(tmpR0a_); - streamwiseGrad(dim_, tmpR0_gf_, vel_gf_, tmpR1_gf_); + + // compute streamwise gradient of input field + // tmpR0_gf_.SetFromTrueDofs(phi); + // streamwiseGrad(dim_, tmpR0_gf_, vel_gf_, tmpR1_gf_); + + tmpR1_gf_.SetFromTrueDofs(gradPhi); + streamwiseGrad(dim_, vel_gf_, tmpR1_gf_); // divergence of sw-grad tmpR1_gf_.GetTrueDofs(tmpR1_); D_op_->Mult(tmpR1_, swDiff); gridScale_gf_->GetTrueDofs(tmpR0b_); - (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + //(turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + (flow_interface_->Reh)->GetTrueDofs(tmpR0c_); const double *rho = rn_.HostRead(); const double *vel = tmpR0a_.HostRead(); const double *del = tmpR0b_.HostRead(); - const double *mu = visc_.HostRead(); - const double *muT = tmpR0c_.HostRead(); + // const double *mu = visc_.HostRead(); + // const double *muT = tmpR0c_.HostRead(); + const double *Reh = tmpR0c_.HostRead(); double *data = swDiff.HostReadWrite(); int Sdof = rn_.Size(); @@ -1304,7 +1389,8 @@ void LteThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { Umag = std::sqrt(Umag); // element Re - double Re = Umag * del[dof] * rho[dof] / (mu[dof] + muT[dof]); + // double Re = Umag * del[dof] * rho[dof] / (mu[dof] + muT[dof]); + double Re = Reh[dof]; // SUPG weight double Csupg = 0.5 * (tanh(re_factor_ * Re - re_offset_) + 1.0); diff --git a/src/lte_thermo_chem.hpp b/src/lte_thermo_chem.hpp index 2aaed5e61..ff634b3a7 100644 --- a/src/lte_thermo_chem.hpp +++ b/src/lte_thermo_chem.hpp @@ -92,10 +92,34 @@ class LteThermoChem final : public ThermoChemModelBase { bool axisym_ = false; // Linear-solver-related options + int smoother_poly_order_; + double smoother_poly_fraction_ = 0.75; + int smoother_eig_est_ = 10; + int smoother_passes_ = 1; + double smoother_relax_weight_ = 0.4; + double smoother_relax_omega_ = 1.0; + double hsmoother_relax_weight_ = 0.8; + double hsmoother_relax_omega_ = 0.1; + + // solver tolerance options int pl_solve_ = 0; /**< Verbosity level passed to mfem solvers */ int max_iter_; /**< Maximum number of linear solver iterations */ double rtol_ = 1e-12; /**< Linear solver relative tolerance */ + int default_max_iter_ = 1000; + double default_rtol_ = 1.0e-10; + double default_atol_ = 1.0e-12; + + int mass_inverse_pl_ = 0; + int mass_inverse_max_iter_; + double mass_inverse_rtol_; + double mass_inverse_atol_; + + int hsolve_pl_ = 0; + int hsolve_max_iter_; + double hsolve_rtol_; + double hsolve_atol_; + // streamwise-stabilization bool sw_stab_; double re_offset_; @@ -198,6 +222,7 @@ class LteThermoChem final : public ThermoChemModelBase { // operators and solvers ParBilinearForm *At_form_ = nullptr; ParBilinearForm *Ms_form_ = nullptr; + ParBilinearForm *Mv_form_ = nullptr; ParBilinearForm *M_rho_Cp_form_ = nullptr; ParBilinearForm *Ht_form_ = nullptr; @@ -211,16 +236,19 @@ class LteThermoChem final : public ThermoChemModelBase { ParLinearForm *LQ_bdry_ = nullptr; ParMixedBilinearForm *D_form_ = nullptr; + ParMixedBilinearForm *G_form_ = nullptr; OperatorHandle At_; OperatorHandle Ht_; OperatorHandle Ms_; + OperatorHandle Mv_; OperatorHandle Mq_; OperatorHandle LQ_; OperatorHandle M_rho_Cp_; OperatorHandle M_rho_; OperatorHandle A_rho_; OperatorHandle D_op_; + OperatorHandle G_op_; mfem::Solver *MsInvPC_ = nullptr; mfem::CGSolver *MsInv_ = nullptr; @@ -230,6 +258,8 @@ class LteThermoChem final : public ThermoChemModelBase { mfem::CGSolver *MrhoInv_ = nullptr; mfem::Solver *HtInvPC_ = nullptr; mfem::CGSolver *HtInv_ = nullptr; + mfem::Solver *Mv_inv_pc_ = nullptr; + mfem::CGSolver *Mv_inv_ = nullptr; // Vectors Vector Tn_, Tn_next_, Tnm1_, Tnm2_; @@ -239,6 +269,7 @@ class LteThermoChem final : public ThermoChemModelBase { Vector tmpR0_, tmpR0a_, tmpR0b_, tmpR0c_; Vector tmpR1_; Vector swDiff_; + Vector gradT_; Vector Qt_; Vector rn_, rnm1_, rnm2_, rnm3_; @@ -282,7 +313,8 @@ class LteThermoChem final : public ThermoChemModelBase { void computeExplicitTempConvectionOP(); void computeQt(); void updateHistory(); - void streamwiseDiffusion(Vector &phi, Vector &swDiff); + // void streamwiseDiffusion(Vector &phi, Vector &swDiff); + void streamwiseDiffusion(Vector &gradPhi, Vector &swDiff); /// Return a pointer to the current temperature ParGridFunction. ParGridFunction *GetCurrentTemperature() { return &Tn_gf_; } diff --git a/src/split_flow_base.hpp b/src/split_flow_base.hpp index 9ef8949f0..3634e3c86 100644 --- a/src/split_flow_base.hpp +++ b/src/split_flow_base.hpp @@ -49,6 +49,8 @@ struct flowToThermoChem { bool swirl_supported = false; const mfem::ParGridFunction *swirl = nullptr; + + const mfem::ParGridFunction *Reh = nullptr; }; struct flowToTurbModel { @@ -60,6 +62,8 @@ struct flowToTurbModel { const mfem::ParGridFunction *gradU = nullptr; const mfem::ParGridFunction *gradV = nullptr; const mfem::ParGridFunction *gradW = nullptr; + + const mfem::ParGridFunction *Reh = nullptr; }; class FlowBase { diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index 4dc42b906..5a8eb9f8e 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -230,6 +230,7 @@ Tomboulides::~Tomboulides() { delete gradU_gf_; delete gradV_gf_; delete gradW_gf_; + delete Reh_gf_; delete tmpR0_gf_; delete tmpR1_gf_; delete vfes_; @@ -242,6 +243,11 @@ void Tomboulides::initializeSelf() { // Initialize minimal state and interface vfec_ = new H1_FECollection(vorder_, dim_); vfes_ = new ParFiniteElementSpace(pmesh_, vfec_, dim_); + // sfec_ = new H1_FECollection(vorder_); + // sfes_ = new ParFiniteElementSpace(pmesh_, sfec_); + pfec_ = new H1_FECollection(porder_); + pfes_ = new ParFiniteElementSpace(pmesh_, pfec_); + u_curr_gf_ = new ParGridFunction(vfes_); u_next_gf_ = new ParGridFunction(vfes_); curl_gf_ = new ParGridFunction(vfes_); @@ -251,20 +257,15 @@ void Tomboulides::initializeSelf() { gradU_gf_ = new ParGridFunction(vfes_); gradV_gf_ = new ParGridFunction(vfes_); gradW_gf_ = new ParGridFunction(vfes_); - sfec_ = new H1_FECollection(vorder_); - sfes_ = new ParFiniteElementSpace(pmesh_, sfec_); - - pfec_ = new H1_FECollection(porder_); - pfes_ = new ParFiniteElementSpace(pmesh_, pfec_); p_gf_ = new ParGridFunction(pfes_); resp_gf_ = new ParGridFunction(pfes_); - mu_total_gf_ = new ParGridFunction(pfes_); pp_div_rad_comp_gf_ = new ParGridFunction(pfes_, *pp_div_gf_); u_next_rad_comp_gf_ = new ParGridFunction(pfes_, *u_next_gf_); - tmpR0_gf_ = new ParGridFunction(sfes_); + Reh_gf_ = new ParGridFunction(pfes_); + tmpR0_gf_ = new ParGridFunction(pfes_); tmpR1_gf_ = new ParGridFunction(vfes_); if (axisym_) { @@ -285,17 +286,20 @@ void Tomboulides::initializeSelf() { *p_gf_ = 0.0; *resp_gf_ = 0.0; *mu_total_gf_ = 0.0; + *Reh_gf_ = 0.0; if (axisym_) { *utheta_gf_ = 0.0; *utheta_next_gf_ = 0.0; } + // exports toThermoChem_interface_.velocity = u_next_gf_; if (axisym_) { toThermoChem_interface_.swirl_supported = true; toThermoChem_interface_.swirl = utheta_next_gf_; } + toThermoChem_interface_.Reh = Reh_gf_; toTurbModel_interface_.velocity = u_next_gf_; if (axisym_) { @@ -305,10 +309,11 @@ void Tomboulides::initializeSelf() { toTurbModel_interface_.gradU = gradU_gf_; toTurbModel_interface_.gradV = gradV_gf_; toTurbModel_interface_.gradW = gradW_gf_; + toTurbModel_interface_.Reh = Reh_gf_; // Allocate Vector storage const int vfes_truevsize = vfes_->GetTrueVSize(); - const int sfes_truevsize = sfes_->GetTrueVSize(); + // const int sfes_truevsize = sfes_->GetTrueVSize(); const int pfes_truevsize = pfes_->GetTrueVSize(); forcing_vec_.SetSize(vfes_truevsize); @@ -345,9 +350,10 @@ void Tomboulides::initializeSelf() { } swDiff_vec_.SetSize(vfes_truevsize); - tmpR0_.SetSize(sfes_truevsize); - tmpR0a_.SetSize(sfes_truevsize); - tmpR0b_.SetSize(sfes_truevsize); + tmpR0_.SetSize(pfes_truevsize); + tmpR0a_.SetSize(pfes_truevsize); + tmpR0b_.SetSize(pfes_truevsize); + tmpR0c_.SetSize(pfes_truevsize); tmpR1_.SetSize(vfes_truevsize); gradU_.SetSize(vfes_truevsize); @@ -799,6 +805,9 @@ void Tomboulides::initializeOperators() { } Mv_form_->Assemble(); Mv_form_->FormSystemMatrix(empty, Mv_op_); + // NOTE: should have dirichlet bc's here for inverse solve, but + // this operator is overloaded and not just used for vel... + // Mv_form_->FormSystemMatrix(vel_ess_tdof_, Mv_op_); // Mass matrix (density weighted) for the velocity Mv_rho_form_ = new ParBilinearForm(vfes_); @@ -1073,6 +1082,7 @@ void Tomboulides::initializeViz(mfem::ParaViewDataCollection &pvdc) const { if (axisym_) { pvdc.RegisterField("swirl", utheta_gf_); } + pvdc.RegisterField("Re_h", Reh_gf_); } void Tomboulides::initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) const { @@ -1360,12 +1370,21 @@ void Tomboulides::step() { } // Add streamwise stability to rhs + computeReh(); if (sw_stab_) { + /* for (int i = 0; i < dim_; i++) { setScalarFromVector(u_vec_, i, &tmpR0a_); streamwiseDiffusion(tmpR0a_, tmpR0b_); setVectorFromScalar(tmpR0b_, i, &swDiff_vec_); } + */ + streamwiseDiffusion(gradU_, tmpR0b_); + setVectorFromScalar(tmpR0b_, 0, &swDiff_vec_); + streamwiseDiffusion(gradV_, tmpR0b_); + setVectorFromScalar(tmpR0b_, 1, &swDiff_vec_); + streamwiseDiffusion(gradW_, tmpR0b_); + setVectorFromScalar(tmpR0b_, 2, &swDiff_vec_); Mv_rho_inv_->Mult(swDiff_vec_, tmpR1_); pp_div_vec_ += tmpR1_; } @@ -1715,11 +1734,42 @@ void Tomboulides::evaluateVelocityGradient() { gradW_gf_->SetFromTrueDofs(gradW_); } +void Tomboulides::computeReh() { + (thermo_interface_->density)->GetTrueDofs(rho_vec_); + gridScale_gf_->GetTrueDofs(tmpR0_); + mu_total_gf_->GetTrueDofs(mu_vec_); + // u_curr_gf_->GetTrueDofs(u_vec_); + u_next_gf_->GetTrueDofs(uext_vec_); + + const double *rho = rho_vec_.HostRead(); + const double *del = tmpR0_.HostRead(); + const double *vel = uext_vec_.HostRead(); + const double *mu = mu_vec_.HostRead(); + double *data = tmpR0c_.HostReadWrite(); + + int Sdof = tmpR0c_.Size(); + for (int dof = 0; dof < Sdof; dof++) { + // vel mag + double Umag = 0.0; + for (int i = 0; i < dim_; i++) Umag += vel[dof + i * Sdof] * vel[dof + i * Sdof]; + Umag = std::sqrt(Umag); + + // element Re + double Re = Umag * del[dof] * rho[dof] / mu[dof]; + data[dof] = Re; + } + Reh_gf_->SetFromTrueDofs(tmpR0c_); +} + // f(Re_h) * Mu_sw * div(streamwiseGrad), i.e. this does not consider grad of f(Re_h) * Mu_sw -void Tomboulides::streamwiseDiffusion(Vector &phi, Vector &swDiff) { +// void Tomboulides::streamwiseDiffusion(Vector &phi, Vector &swDiff) { +void Tomboulides::streamwiseDiffusion(Vector &gradPhi, Vector &swDiff) { // compute streamwise gradient of input field - tmpR0_gf_->SetFromTrueDofs(phi); - streamwiseGrad(dim_, *tmpR0_gf_, *u_curr_gf_, *tmpR1_gf_); + // tmpR0_gf_->SetFromTrueDofs(phi); + // streamwiseGrad(dim_, *tmpR0_gf_, *u_curr_gf_, *tmpR1_gf_); + + tmpR1_gf_->SetFromTrueDofs(gradPhi); + streamwiseGrad(dim_, *u_curr_gf_, *tmpR1_gf_); // divergence of sw-grad tmpR1_gf_->GetTrueDofs(tmpR1_); @@ -1727,12 +1777,16 @@ void Tomboulides::streamwiseDiffusion(Vector &phi, Vector &swDiff) { (thermo_interface_->density)->GetTrueDofs(rho_vec_); gridScale_gf_->GetTrueDofs(tmpR0_); - mu_total_gf_->GetTrueDofs(mu_vec_); + Reh_gf_->GetTrueDofs(tmpR0c_); + // mu_total_gf_->GetTrueDofs(mu_vec_); + upwindDiff(dim_, re_factor_, re_offset_, u_vec_, rho_vec_, tmpR0_, tmpR0c_, swDiff); + /* const double *rho = rho_vec_.HostRead(); const double *del = tmpR0_.HostRead(); const double *vel = u_vec_.HostRead(); - const double *mu = mu_vec_.HostRead(); + //const double *mu = mu_vec_.HostRead(); + const double *Reh = tmpR0c_.HostRead(); double *data = swDiff.HostReadWrite(); int Sdof = rho_vec_.Size(); @@ -1742,7 +1796,8 @@ void Tomboulides::streamwiseDiffusion(Vector &phi, Vector &swDiff) { Umag = std::sqrt(Umag); // element Re - double Re = Umag * del[dof] * rho[dof] / mu[dof]; + //double Re = Umag * del[dof] * rho[dof] / mu[dof]; + double Re = Reh[dof]; // SUPG weight double Csupg = 0.5 * (tanh(re_factor_ * Re - re_offset_) + 1.0); @@ -1753,4 +1808,5 @@ void Tomboulides::streamwiseDiffusion(Vector &phi, Vector &swDiff) { // scaled streamwise Laplacian data[dof] *= CswDiff; } + */ } diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp index 7b408d3fe..1010d7c8a 100644 --- a/src/tomboulides.hpp +++ b/src/tomboulides.hpp @@ -229,6 +229,7 @@ class Tomboulides final : public FlowBase { mfem::ParGridFunction *gradW_gf_ = nullptr; // mfem::ParGridFunction *buffer_uInlet_ = nullptr; mfem::VectorGridFunctionCoefficient *velocity_field_ = nullptr; + mfem::ParGridFunction *Reh_gf_ = nullptr; mfem::ParGridFunction *tmpR0_gf_ = nullptr; mfem::ParGridFunction *tmpR1_gf_ = nullptr; @@ -367,6 +368,7 @@ class Tomboulides final : public FlowBase { mfem::Vector tmpR0_; mfem::Vector tmpR0a_; mfem::Vector tmpR0b_; + mfem::Vector tmpR0c_; mfem::Vector tmpR1_; mfem::Vector gradU_; mfem::Vector gradV_; @@ -447,6 +449,11 @@ class Tomboulides final : public FlowBase { */ void streamwiseDiffusion(Vector &phi, Vector &swDiff); + /** + * @brief Computes element convective Reynolds number + */ + void computeReh(); + /// Advance void step() final; diff --git a/src/utils.cpp b/src/utils.cpp index 5843fa626..8bb9ac27a 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -763,6 +763,7 @@ void ComputeCurl3D(const ParGridFunction &u, ParGridFunction &cu) { // Communication + /* // Count the zones globally. GroupCommunicator &gcomm = u.ParFESpace()->GroupComm(); gcomm.Reduce(zones_per_vdof, GroupCommunicator::Sum); @@ -771,6 +772,7 @@ void ComputeCurl3D(const ParGridFunction &u, ParGridFunction &cu) { // Accumulate for all vdofs. gcomm.Reduce(cu.GetData(), GroupCommunicator::Sum); gcomm.Bcast(cu.GetData()); + */ // Compute means. for (int i = 0; i < cu.Size(); i++) { @@ -894,6 +896,7 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { } } + /* // Count the zones globally. GroupCommunicator &gcomm = u.ParFESpace()->GroupComm(); gcomm.Reduce(zones_per_vdof, GroupCommunicator::Sum); @@ -902,6 +905,7 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { // Accumulate for all vdofs. gcomm.Reduce(gu.GetData(), GroupCommunicator::Sum); gcomm.Bcast(gu.GetData()); + */ // Compute means. for (int dir = 0; dir < dim_; dir++) { @@ -985,6 +989,7 @@ void ComputeCurl2D(const ParGridFunction &u, ParGridFunction &cu, bool assume_sc // Communication. + /* // Count the zones globally. GroupCommunicator &gcomm = u.ParFESpace()->GroupComm(); gcomm.Reduce(zones_per_vdof, GroupCommunicator::Sum); @@ -993,6 +998,7 @@ void ComputeCurl2D(const ParGridFunction &u, ParGridFunction &cu, bool assume_sc // Accumulate for all vdofs. gcomm.Reduce(cu.GetData(), GroupCommunicator::Sum); gcomm.Bcast(cu.GetData()); + */ // Compute means. for (int i = 0; i < cu.Size(); i++) { @@ -1088,6 +1094,7 @@ void ComputeCurlAxi(const ParGridFunction &u, ParGridFunction &cu, bool assume_s // Communication. + /* // Count the zones globally. GroupCommunicator &gcomm = u.ParFESpace()->GroupComm(); gcomm.Reduce(zones_per_vdof, GroupCommunicator::Sum); @@ -1096,6 +1103,7 @@ void ComputeCurlAxi(const ParGridFunction &u, ParGridFunction &cu, bool assume_s // Accumulate for all vdofs. gcomm.Reduce(cu.GetData(), GroupCommunicator::Sum); gcomm.Bcast(cu.GetData()); + */ // Compute means. for (int i = 0; i < cu.Size(); i++) { @@ -1165,7 +1173,8 @@ bool copyFile(const char *SRC, const char *DEST) { return src && dest; } -void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad) { +// void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad) { +void streamwiseGrad(int dim, ParGridFunction &u, ParGridFunction &swGrad) { /* std::cout << "maxInd minusInd plusInd" << endl; std::cout << 0 << " " << ((0-1) % dim_ + dim_) % dim_ << " " << (0 + 1) % dim_ << endl; @@ -1174,12 +1183,12 @@ void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFu */ // compute gradient of input field - scalarGrad3D(phi, swGrad); + // scalarGrad3D(phi, swGrad); const double *vel = u.HostRead(); double *gPhi = swGrad.HostReadWrite(); - int Sdof = phi.Size(); + int Sdof = u.Size() / dim; for (int dof = 0; dof < Sdof; dof++) { // streamwise coordinate system Vector unitNorm; @@ -1253,7 +1262,7 @@ void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFu } */ - // muSWgbl = M_{im} muSw_{mn} M_{jn} or M*mu*M^T (with n,t1,t2 in columns of M) + // M_{im} swM_{mn} M_{jn} or M*"mu"*M^T (with n,t1,t2 in columns of M) DenseMatrix swMgbl(dim, dim); swMgbl = 0.0; for (int i = 0; i < dim; i++) { @@ -1266,11 +1275,12 @@ void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFu } } - // mu*gPhi + // copy grad into local vecotr Vector tmp1; tmp1.SetSize(dim); for (int i = 0; i < dim; i++) tmp1[i] = gPhi[dof + i * Sdof]; + // gradient in streamwise-direction Vector tmp2; tmp2.SetSize(dim); for (int i = 0; i < dim; i++) tmp2[i] = 0.0; @@ -1285,6 +1295,34 @@ void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFu } } +void upwindDiff(int dim, double re_factor, double re_offset, Vector &u_vec, Vector &rho_vec, Vector &del_vec, + Vector &Reh_vec, Vector &swDiff) { + const double *rho = rho_vec.HostRead(); + const double *del = del_vec.HostRead(); + const double *vel = u_vec.HostRead(); + const double *Reh = Reh_vec.HostRead(); + double *data = swDiff.HostReadWrite(); + + int Sdof = rho_vec.Size(); + for (int dof = 0; dof < Sdof; dof++) { + double Umag = 0.0; + for (int i = 0; i < dim; i++) Umag += vel[i] * vel[i]; + Umag = std::sqrt(Umag); + + // element Re + double Re = Reh[dof]; + + // SUPG weight + double Csupg = 0.5 * (tanh(re_factor * Re - re_offset) + 1.0); + + // streamwise diffusion coeff + double CswDiff = Csupg * Umag * del[dof] * rho[dof]; + + // scaled streamwise Laplacian + data[dof] *= CswDiff; + } +} + namespace mfem { GradientVectorGridFunctionCoefficient::GradientVectorGridFunctionCoefficient(const GridFunction *gf) : MatrixCoefficient((gf) ? gf->VectorDim() : 0) { diff --git a/src/utils.hpp b/src/utils.hpp index 603e0804b..8ae3877f8 100644 --- a/src/utils.hpp +++ b/src/utils.hpp @@ -174,7 +174,10 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu); void vectorGrad3DV(FiniteElementSpace *fes, Vector u, Vector *gu, Vector *gv, Vector *gw); void scalarGrad3DV(FiniteElementSpace *fes, FiniteElementSpace *vfes, Vector u, Vector *gu); bool copyFile(const char *SRC, const char *DEST); -void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad); +// void streamwiseGrad(int dim, ParGridFunction &phi, ParGridFunction &u, ParGridFunction &swGrad); +void streamwiseGrad(int dim, ParGridFunction &u, ParGridFunction &swGrad); +void upwindDiff(int dim, double re_factor, double re_offset, Vector &u_vec, Vector &rho_vec, Vector &del_vec, + Vector &Reh_vec, Vector &swDiff); /// Eliminate essential BCs in an Operator and apply to RHS. /// rename this to something sensible "ApplyEssentialBC" or something From ed1b8484faee24a851f8baa661c9850db9db83cd Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 10 Jul 2024 14:44:18 -0700 Subject: [PATCH 06/12] fix to scalargrad3d comm and reinstatment of curl calc comms --- src/utils.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/utils.cpp b/src/utils.cpp index 8bb9ac27a..8caf58eed 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -763,7 +763,7 @@ void ComputeCurl3D(const ParGridFunction &u, ParGridFunction &cu) { // Communication - /* + /**/ // Count the zones globally. GroupCommunicator &gcomm = u.ParFESpace()->GroupComm(); gcomm.Reduce(zones_per_vdof, GroupCommunicator::Sum); @@ -772,7 +772,7 @@ void ComputeCurl3D(const ParGridFunction &u, ParGridFunction &cu) { // Accumulate for all vdofs. gcomm.Reduce(cu.GetData(), GroupCommunicator::Sum); gcomm.Bcast(cu.GetData()); - */ + /**/ // Compute means. for (int i = 0; i < cu.Size(); i++) { @@ -896,16 +896,14 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { } } - /* // Count the zones globally. - GroupCommunicator &gcomm = u.ParFESpace()->GroupComm(); + GroupCommunicator &gcomm = gu.ParFESpace()->GroupComm(); gcomm.Reduce(zones_per_vdof, GroupCommunicator::Sum); gcomm.Bcast(zones_per_vdof); // Accumulate for all vdofs. gcomm.Reduce(gu.GetData(), GroupCommunicator::Sum); gcomm.Bcast(gu.GetData()); - */ // Compute means. for (int dir = 0; dir < dim_; dir++) { From 650e116c81c3e38000f877109c662e5125119658 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 24 Jul 2024 10:46:25 -0700 Subject: [PATCH 07/12] added error msg when stabalization requested and not 3d (add support later) and added test case --- src/tomboulides.cpp | 6 ++ test/.gitattributes | 2 + test/Makefile.am | 6 +- test/inputs/input.stabChan.ini | 96 +++++++++++++++++++ test/lomach-chan-stab.test | 30 ++++++ test/meshes/channel182p4_25x9p4_THIN.msh | 3 + test/ref_solns/stabChan/restart_output.sol.h5 | 3 + 7 files changed, 144 insertions(+), 2 deletions(-) create mode 100644 test/inputs/input.stabChan.ini create mode 100755 test/lomach-chan-stab.test create mode 100755 test/meshes/channel182p4_25x9p4_THIN.msh create mode 100644 test/ref_solns/stabChan/restart_output.sol.h5 diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index 5a8eb9f8e..3517eca96 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -135,6 +135,12 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS tps->getInput("loMach/tomboulides/streamwise-stabilization", sw_stab_, false); tps->getInput("loMach/tomboulides/Reh_offset", re_offset_, 1.0); tps->getInput("loMach/tomboulides/Reh_factor", re_factor_, 0.01); + + if ((sw_stab_ == true) && (dim_ != 3)) { + if (rank0_) std::cout << " ERROR: streamwise stabalization only implemented for 3D simulations. exiting..." << endl; + assert(false); + exit(1); + } } Tomboulides::~Tomboulides() { diff --git a/test/.gitattributes b/test/.gitattributes index 0c1274b58..f943ebe51 100644 --- a/test/.gitattributes +++ b/test/.gitattributes @@ -53,3 +53,5 @@ ref_solns/spongeBox/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text meshes/spongeBox.msh filter=lfs diff=lfs merge=lfs -text ref_solns/lequere-varmu/restart_output-lequere-varmu.sol.h5 filter=lfs diff=lfs merge=lfs -text ref_solns/lequere-varmu/reference-lequere-varmu.sol.h5 filter=lfs diff=lfs merge=lfs -text +meshes/channel182p4_25x9p4_THIN.msh filter=lfs diff=lfs merge=lfs -text +ref_solns/stabChan/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text diff --git a/test/Makefile.am b/test/Makefile.am index 9d1314f3d..decc49ce4 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -34,6 +34,7 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- ref_solns/aveLoMach/*.h5 \ ref_solns/reactBinDiff/*.h5 \ ref_solns/reactSingleRx/*.h5 \ + ref_solns/stabChan/*.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 \ @@ -49,7 +50,7 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- 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 aveLoMach.test \ - reactFlow-binDiff.test reactFlow-singleRx.test + reactFlow-binDiff.test reactFlow-singleRx.test lomach-chan-stab.test TESTS = vpath.sh XFAIL_TESTS = @@ -117,7 +118,8 @@ TESTS += cyl3d.test \ autoPeriodic.test \ aveLoMach.test \ reactFlow-binDiff.test \ - reactFlow-singleRx.test + reactFlow-singleRx.test \ + loMach-chan-stab.test if PYTHON_ENABLED TESTS += cyl3d.python.test \ diff --git a/test/inputs/input.stabChan.ini b/test/inputs/input.stabChan.ini new file mode 100644 index 000000000..f0f020c28 --- /dev/null +++ b/test/inputs/input.stabChan.ini @@ -0,0 +1,96 @@ +[solver] +type = loMach + +[loMach] +mesh = meshes/channel182p4_25x9p4_THIN.msh +order = 1 +nFilter = 0 +filterWeight = 1.0 +maxIters = 100 +outputFreq = 100 +fluid = dry_air +refLength = 1.0 +equation_system = navier-stokes +enablePressureForcing = True +pressureGrad = '0.00005372 0.0 0.0' +enableGravity = False +gravity = '0.0 0.0 0.0' +openSystem = False +sgsModel = none +flow-solver = tomboulides +thermo-solver = calorically-perfect + +[loMach/calperfect] +viscosity-model = sutherland +sutherland/mu0 = 1.68e-5 +sutherland/T0 = 273.0 +sutherland/S0 = 110.4 +numerical-integ = false +Prandtl = 0.72 +#ic = channel +hsolve-atol = 1.0e-12 +msolve-atol = 1.0e-12 +hsolve-rtol = 1.0e-10 +msolve-rtol = 1.0e-10 +hsolve-maxIters = 2000 +msolve-maxIters = 2000 +streamwise-stabilization = true +Reh_offset = 1.0 +Reh_factor = 0.01 + +[loMach/tomboulides] +ic = channel +numerical-integ = false +psolve-atol = 1.0e-14 +hsolve-atol = 1.0e-12 +msolve-atol = 1.0e-12 +psolve-rtol = 1.0e-10 +hsolve-rtol = 1.0e-10 +msolve-rtol = 1.0e-10 +psolve-maxIters = 2000 +hsolve-maxIters = 2000 +msolve-maxIters = 2000 +streamwise-stabilization = true +Reh_offset = 1.0 +Reh_factor = 0.001 + +[io] +outdirBase = output +#enableRestart = True +#restartMode = variableP + +[time] +integrator = curlcurl +cfl = 0.4 +dt_initial = 1.0e-4 +dtFactor = 0.01 +#dt_fixed = 1.0e-4 +bdfOrder = 2 + +[spongeMultiplier] +uniform = true +uniformMult = 1.0 + +[initialConditions] +rho = 1.0 +rhoU = 1.0 +rhoV = 0.0 +rhoW = 0.0 +temperature = 300.0 +pressure = 101325.0 + +[boundaryConditions/wall1] +patch = 4 +type = viscous_isothermal +temperature = 300 +velocity = '0.0 0.0 0.0' + +[boundaryConditions] +numWalls = 1 +numInlets = 0 +numOutlets = 0 + +[periodicity] +enablePeriodic = True +periodicX = True +periodicZ = True diff --git a/test/lomach-chan-stab.test b/test/lomach-chan-stab.test new file mode 100755 index 000000000..22580cb18 --- /dev/null +++ b/test/lomach-chan-stab.test @@ -0,0 +1,30 @@ +#!./bats +# -*- mode: sh -*- + +TEST="stabChan" +RUNFILE="inputs/input.stabChan.ini" +EXE="../src/tps" +RESTART="ref_solns/stabChan/restart_output.sol.h5" + +setup() { + SOLN_FILE=restart_output.sol.h5 + REF_FILE=ref_solns/stabChan/restart_output.sol.h5 + OUT_FILE=output_solns/restart_output_stabChan.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 + $EXE --runFile $RUNFILE + test -s $SOLN_FILE +} + +@test "[$TEST] verify tps output with input -> $RUNFILE" { + test -s $SOLN_FILE + test -s $REF_FILE + h5diff -r --relative=1e-10 $SOLN_FILE $REF_FILE /velocity +} diff --git a/test/meshes/channel182p4_25x9p4_THIN.msh b/test/meshes/channel182p4_25x9p4_THIN.msh new file mode 100755 index 000000000..ec72d1aa7 --- /dev/null +++ b/test/meshes/channel182p4_25x9p4_THIN.msh @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b8bba9f4101291f8cfce812f02d669bdbe7178b74674a08eb745bb96531cd9bb +size 673400 diff --git a/test/ref_solns/stabChan/restart_output.sol.h5 b/test/ref_solns/stabChan/restart_output.sol.h5 new file mode 100644 index 000000000..c8b8feb37 --- /dev/null +++ b/test/ref_solns/stabChan/restart_output.sol.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2737ebf44db149c10aab0c9ecc049448df6ae6a1e8ff3d77de5b3314408c0ab3 +size 150168 From 01fcfd5da00c9ca144098836454c11fdca3eaa24 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 24 Jul 2024 11:49:06 -0700 Subject: [PATCH 08/12] fix to epsi_gf_ space --- src/tomboulides.cpp | 10 +++------- src/utils.cpp | 2 +- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index 216f7cb82..bbae2c392 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -266,7 +266,7 @@ void Tomboulides::initializeSelf() { gradW_gf_ = new ParGridFunction(vfes_); p_gf_ = new ParGridFunction(pfes_); resp_gf_ = new ParGridFunction(pfes_); - epsi_gf_ = new ParGridFunction(sfes_); + epsi_gf_ = new ParGridFunction(pfes_); mu_total_gf_ = new ParGridFunction(pfes_); // pp_div_rad_comp_gf_ = new ParGridFunction(pfes_, *pp_div_gf_); @@ -1510,9 +1510,7 @@ void Tomboulides::step() { // systems. As a workaround, we copy instead. auto d_pp_div_rad = pp_div_rad_comp_gf_->Write(); auto d_pp_div = pp_div_gf_->Read(); - MFEM_FORALL(i, pp_div_rad_comp_gf_->Size(), { - d_pp_div_rad[i] = d_pp_div[i]; - }); + MFEM_FORALL(i, pp_div_rad_comp_gf_->Size(), { d_pp_div_rad[i] = d_pp_div[i]; }); } pp_div_rad_comp_gf_->HostRead(); @@ -1628,9 +1626,7 @@ void Tomboulides::step() { // about pp_div_rad_comp_gf_ above. auto d_u_next_rad = u_next_rad_comp_gf_->Write(); auto d_u_next = u_next_gf_->Read(); - MFEM_FORALL(i, u_next_rad_comp_gf_->Size(), { - d_u_next_rad[i] = d_u_next[i]; - }); + MFEM_FORALL(i, u_next_rad_comp_gf_->Size(), { d_u_next_rad[i] = d_u_next[i]; }); } u_next_rad_comp_gf_->HostRead(); diff --git a/src/utils.cpp b/src/utils.cpp index 0a5a0fdb2..fad303164 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1336,7 +1336,7 @@ void upwindDiff(int dim, double re_factor, double re_offset, Vector &u_vec, Vect data[dof] *= CswDiff; } } - + void makeContinuous(ParGridFunction &u) { FiniteElementSpace *fes = u.FESpace(); From 31c597d90fe2414b8b4c2985d8916b9e9ac9adc5 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 24 Jul 2024 11:54:06 -0700 Subject: [PATCH 09/12] format checks --- src/calorically_perfect.cpp | 2 +- src/lte_thermo_chem.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/calorically_perfect.cpp b/src/calorically_perfect.cpp index b93ba7fd1..aa7ef0805 100644 --- a/src/calorically_perfect.cpp +++ b/src/calorically_perfect.cpp @@ -1156,7 +1156,7 @@ void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &gradPhi, Vector & D_op_->Mult(tmpR1_, swDiff); gridScale_gf_->GetTrueDofs(tmpR0b_); - //(turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + // (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); (flow_interface_->Reh)->GetTrueDofs(tmpR0c_); upwindDiff(dim_, re_factor_, re_offset_, tmpR0a_, rn_, tmpR0b_, tmpR0c_, swDiff); diff --git a/src/lte_thermo_chem.cpp b/src/lte_thermo_chem.cpp index 482305692..6a90ce58a 100644 --- a/src/lte_thermo_chem.cpp +++ b/src/lte_thermo_chem.cpp @@ -1371,7 +1371,7 @@ void LteThermoChem::streamwiseDiffusion(Vector &gradPhi, Vector &swDiff) { D_op_->Mult(tmpR1_, swDiff); gridScale_gf_->GetTrueDofs(tmpR0b_); - //(turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); + // (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_); (flow_interface_->Reh)->GetTrueDofs(tmpR0c_); const double *rho = rn_.HostRead(); From 10c7a689959fb25394561d024c84fac9b6f484d5 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 24 Jul 2024 16:32:04 -0700 Subject: [PATCH 10/12] removed test_tomboulides from tests makefile as this is no longer necessary and conflicts with new interface --- test/Makefile.am | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/Makefile.am b/test/Makefile.am index 322f67f85..66407e2d7 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -182,11 +182,11 @@ split_interface_test_LDADD = ../src/libtps.la split_interface_test_LDADD += $(HDF5_LIBS) split_interface_test_LDADD += $(GRVY_LIBS) -check_PROGRAMS += tomboulides_test -tomboulides_test_SOURCES = test_tomboulides.cpp -tomboulides_test_LDADD = ../src/libtps.la -tomboulides_test_LDADD += $(HDF5_LIBS) -tomboulides_test_LDADD += $(GRVY_LIBS) +#check_PROGRAMS += tomboulides_test +#tomboulides_test_SOURCES = test_tomboulides.cpp +#tomboulides_test_LDADD = ../src/libtps.la +#tomboulides_test_LDADD += $(HDF5_LIBS) +#tomboulides_test_LDADD += $(GRVY_LIBS) if !GPU_ENABLED From d4758cf018ab48692355a67db3b674f768792575 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 24 Jul 2024 19:41:06 -0700 Subject: [PATCH 11/12] typo fix in test makefile --- test/Makefile.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Makefile.am b/test/Makefile.am index 66407e2d7..a1edb2df4 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -122,7 +122,7 @@ TESTS += cyl3d.test \ aveLoMach.test \ reactFlow-binDiff.test \ reactFlow-singleRx.test \ - loMach-chan-stab.test + lomach-chan-stab.test if PYTHON_ENABLED TESTS += cyl3d.python.test \ From 467b6eaad839f123f6b3bbfe583c442f0db9c5f5 Mon Sep 17 00:00:00 2001 From: Sigfried Haering Date: Wed, 24 Jul 2024 19:42:15 -0700 Subject: [PATCH 12/12] reduced tol on new test case --- test/lomach-chan-stab.test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/lomach-chan-stab.test b/test/lomach-chan-stab.test index 22580cb18..b5af314cc 100755 --- a/test/lomach-chan-stab.test +++ b/test/lomach-chan-stab.test @@ -26,5 +26,5 @@ setup() { @test "[$TEST] verify tps output with input -> $RUNFILE" { test -s $SOLN_FILE test -s $REF_FILE - h5diff -r --relative=1e-10 $SOLN_FILE $REF_FILE /velocity + h5diff -r --delta=1e-10 $SOLN_FILE $REF_FILE /velocity }