diff --git a/src/calorically_perfect.cpp b/src/calorically_perfect.cpp index c182c14e6..aa7ef0805 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/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() { @@ -147,10 +153,15 @@ 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_; delete un_next_coeff_; @@ -167,6 +178,8 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() { // allocated in initializeSelf delete sfes_; delete sfec_; + delete vfes_; + delete vfec_; } void CaloricallyPerfectThermoChem::initializeSelf() { @@ -178,6 +191,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 +205,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 +255,19 @@ 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); + + gradT_.SetSize(vfes_truevsize); + gradT_ = 0.0; R0PM0_gf_.SetSpace(sfes_); @@ -471,6 +499,61 @@ 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_); + + // 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_); @@ -483,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_); @@ -517,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_); @@ -666,6 +771,17 @@ void CaloricallyPerfectThermoChem::step() { tmpR0_ = (dtP_ / Cp_); Ms_->AddMult(tmpR0_, resT_); + // Add streamwise stability to rhs + if (sw_stab_) { + // 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_); + } + // 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 +1138,56 @@ void CaloricallyPerfectThermoChem::screenValues(std::vector &values) { values[0] = thermo_pressure_ / ambient_pressure_; } } + +// void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { +void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &gradPhi, Vector &swDiff) { + (flow_interface_->velocity)->GetTrueDofs(tmpR0a_); + vel_gf_.SetFromTrueDofs(tmpR0a_); + + // 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_); + (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 *Reh = 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]); + 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; + } + */ +} diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp index 118cd3dd5..3f6f72e75 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,11 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { double hsolve_rtol_; double hsolve_atol_; + // 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) */ @@ -144,6 +150,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 +165,11 @@ class CaloricallyPerfectThermoChem : 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; @@ -174,10 +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_; @@ -186,6 +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; @@ -193,13 +214,18 @@ 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_; Vector NTn_, NTnm1_, NTnm2_; Vector Text_; Vector resT_; - Vector tmpR0_, tmpR0b_; + Vector tmpR0_, tmpR0a_, tmpR0b_, tmpR0c_; + Vector tmpR1_; + Vector swDiff_; + Vector gradT_; Vector Qt_; Vector rn_; @@ -233,7 +259,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 +283,8 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase { void computeExplicitTempConvectionOP(bool extrap); void computeQt(); void computeQtTO(); + // 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 5f539e21a..18204cfd1 100644 --- a/src/cases.cpp +++ b/src/cases.cpp @@ -204,7 +204,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/loMach.cpp b/src/loMach.cpp index 9f95926f7..c391b8379 100644 --- a/src/loMach.cpp +++ b/src/loMach.cpp @@ -166,9 +166,10 @@ 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_); + 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 { @@ -185,7 +186,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/lte_thermo_chem.cpp b/src/lte_thermo_chem.cpp index c41a48ebb..6a90ce58a 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,21 @@ 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); + + // 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); } LteThermoChem::~LteThermoChem() { @@ -187,11 +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_; @@ -208,6 +229,8 @@ LteThermoChem::~LteThermoChem() { // allocated in initializeSelf delete sfes_; delete sfec_; + delete vfes_; + delete vfec_; } void LteThermoChem::initializeSelf() { @@ -219,6 +242,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 +256,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 +344,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 +663,61 @@ 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_); + + // 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; @@ -684,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); @@ -905,6 +1016,17 @@ void LteThermoChem::step() { tmpR0_ = dtP_; Ms_->AddMult(tmpR0_, resT_); + // Add streamwise stability to rhs + if (sw_stab_) { + // 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_); + } + // Joule heating (and radiation sink) jh_form_->Update(); jh_form_->Assemble(); @@ -1231,3 +1353,52 @@ void LteThermoChem::computeQt() { Qt_.Neg(); Qt_gf_.SetFromTrueDofs(Qt_); } + +// void LteThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) { +void LteThermoChem::streamwiseDiffusion(Vector &gradPhi, Vector &swDiff) { + (flow_interface_->velocity)->GetTrueDofs(tmpR0a_); + vel_gf_.SetFromTrueDofs(tmpR0a_); + + // 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_); + (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 *Reh = 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]); + 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; + } +} diff --git a/src/lte_thermo_chem.hpp b/src/lte_thermo_chem.hpp index 42e5114bf..ff634b3a7 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_; @@ -91,10 +92,39 @@ 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_; + 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 +161,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 +182,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; @@ -183,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; @@ -195,14 +235,20 @@ class LteThermoChem final : public ThermoChemModelBase { ParBilinearForm *LQ_form_ = nullptr; 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; @@ -212,13 +258,18 @@ 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_; Vector NTn_, NTnm1_, NTnm2_; Vector Text_; Vector resT_; - Vector tmpR0_, tmpR0b_; + Vector tmpR0_, tmpR0a_, tmpR0b_, tmpR0c_; + Vector tmpR1_; + Vector swDiff_; + Vector gradT_; Vector Qt_; Vector rn_, rnm1_, rnm2_, rnm3_; @@ -241,7 +292,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 +313,8 @@ class LteThermoChem final : public ThermoChemModelBase { void computeExplicitTempConvectionOP(); void computeQt(); void updateHistory(); + // 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 9b57c4d99..cbec0575e 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 4e458f9eb..bbae2c392 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,16 @@ 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); + + 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() { @@ -224,6 +236,9 @@ Tomboulides::~Tomboulides() { delete gradU_gf_; delete gradV_gf_; delete gradW_gf_; + delete Reh_gf_; + delete tmpR0_gf_; + delete tmpR1_gf_; delete epsi_gf_; delete vfes_; delete vfec_; @@ -235,6 +250,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_); @@ -244,21 +264,18 @@ 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_); - - 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_); // u_next_rad_comp_gf_ = new ParGridFunction(pfes_, *u_next_gf_); + Reh_gf_ = new ParGridFunction(pfes_); + tmpR0_gf_ = new ParGridFunction(pfes_); + tmpR1_gf_ = new ParGridFunction(vfes_); + if (axisym_) { pp_div_rad_comp_gf_ = new ParGridFunction(pfes_); u_next_rad_comp_gf_ = new ParGridFunction(pfes_); @@ -280,6 +297,7 @@ void Tomboulides::initializeSelf() { *p_gf_ = 0.0; *resp_gf_ = 0.0; *mu_total_gf_ = 0.0; + *Reh_gf_ = 0.0; *epsi_gf_ = 0.0; @@ -288,11 +306,13 @@ void Tomboulides::initializeSelf() { *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_) { @@ -302,10 +322,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); @@ -341,7 +362,11 @@ void Tomboulides::initializeSelf() { utheta_next_vec_.SetSize(pfes_truevsize); } - tmpR0_.SetSize(sfes_truevsize); + swDiff_vec_.SetSize(vfes_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); @@ -793,6 +818,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_); @@ -1067,6 +1095,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 { @@ -1434,6 +1463,26 @@ 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_; + } + // Add ustar/dt contribution pp_div_vec_ += ustar_vec_; @@ -1461,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(); @@ -1539,6 +1586,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); } @@ -1576,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(); @@ -1801,3 +1849,80 @@ void Tomboulides::evaluateVelocityGradient() { gradV_gf_->SetFromTrueDofs(gradV_); 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 &gradPhi, Vector &swDiff) { + // compute streamwise gradient of input field + // 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_); + D_op_->Mult(tmpR1_, swDiff); + + (thermo_interface_->density)->GetTrueDofs(rho_vec_); + gridScale_gf_->GetTrueDofs(tmpR0_); + 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 *Reh = tmpR0c_.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]; + 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; + } + */ +} diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp index 1463b426e..f02339f18 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,9 @@ 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; mfem::ParGridFunction *epsi_gf_ = nullptr; /// Pressure FEM objects and fields @@ -239,6 +246,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; @@ -356,7 +365,11 @@ 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 tmpR0c_; mfem::Vector tmpR1_; mfem::Vector gradU_; mfem::Vector gradV_; @@ -392,7 +405,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; @@ -430,6 +444,17 @@ 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); + + /** + * @brief Computes element convective Reynolds number + */ + void computeReh(); + /** * @brief Compute turbulent dissipation using average u */ diff --git a/src/utils.cpp b/src/utils.cpp index 3d83de4c0..fad303164 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -768,6 +768,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); @@ -776,6 +777,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++) { @@ -865,7 +867,6 @@ 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); @@ -887,6 +888,7 @@ void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu) { int ldof = vdofs[j]; gu(ldof + 1 * nSize) += vals2[j]; } + if (dim == 3) { for (int j = 0; j < vdofs.Size(); j++) { int ldof = vdofs[j]; @@ -994,6 +996,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); @@ -1002,6 +1005,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++) { @@ -1102,6 +1106,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); @@ -1110,6 +1115,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++) { @@ -1181,6 +1187,156 @@ 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 &u, ParGridFunction &swGrad) { + /* + 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); + + const double *vel = u.HostRead(); + double *gPhi = swGrad.HostReadWrite(); + + int Sdof = u.Size() / dim; + 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]; + 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[0] > unitNorm[1] * unitNorm[1]) { + maxInd = 0; + } else { + maxInd = 1; + } + if (dim == 3 && unitNorm[maxInd] * unitNorm[maxInd] < unitNorm[2] * unitNorm[2]) { + maxInd = 2; + } + 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 /= std::sqrt(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 coeff + DenseMatrix swM(dim, dim); + swM = 0.0; + swM(0, 0) = 1.0; + + /* + std::cout << " " << endl; + for (int i = 0; i < dim_; i++) { + for (int j = 0; j < dim_; j++) { + std::cout << M(i,j) << " " ; + } + std::cout << endl; + } + */ + + // 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++) { + 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); + } + } + } + } + + // 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; + 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]; + } +} + +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; + } +} + void makeContinuous(ParGridFunction &u) { FiniteElementSpace *fes = u.FESpace(); diff --git a/src/utils.hpp b/src/utils.hpp index aeee83526..af3e6d78e 100644 --- a/src/utils.hpp +++ b/src/utils.hpp @@ -177,6 +177,10 @@ void scalarGrad3DV(FiniteElementSpace *fes, FiniteElementSpace *vfes, Vector u, void makeContinuous(ParGridFunction &u); bool copyFile(const char *SRC, const char *DEST); +// 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 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 5137f5874..a1edb2df4 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 = @@ -120,7 +121,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 \ @@ -180,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 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..b5af314cc --- /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 --delta=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