From 3153d08dec18e84a4b42d7aca6b1113ff909a0d5 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 18 Sep 2025 17:21:04 -0500 Subject: [PATCH 1/6] Toward supg in Tomboulides --- src/loMach.cpp | 2 +- src/tomboulides.cpp | 127 +++++++++++++++++++++++++++++++++++++++++++- src/tomboulides.hpp | 30 ++++++++++- 3 files changed, 156 insertions(+), 3 deletions(-) diff --git a/src/loMach.cpp b/src/loMach.cpp index 521c8ad3c..64ce5ca97 100644 --- a/src/loMach.cpp +++ b/src/loMach.cpp @@ -186,7 +186,7 @@ void LoMachSolver::initialize() { flow_ = new ZeroFlow(pmesh_, 1, tpsP_); } 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 184df27b6..3902da648 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, + mfem::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)) { @@ -128,6 +130,9 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS tps->getInput("loMach/tomboulides/psolve-maxIters", pressure_solve_max_iter_, default_max_iter_); 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_); + + // artificial diffusion (SUPG) + tpsP_->getInput("loMach/tomboulides/streamwise-stabilization", sw_stab_, false); } } @@ -169,6 +174,8 @@ Tomboulides::~Tomboulides() { delete L_iorho_lor_; delete L_iorho_form_; + delete Mv_stab_form_; + delete swirl_var_viscosity_coeff_; delete utheta_vec_coeff_; delete rho_ur_ut_coeff_; @@ -204,6 +211,19 @@ Tomboulides::~Tomboulides() { delete iorho_coeff_; delete rho_coeff_; + delete umag_coeff_; + delete gscale_coeff_; + delete visc_inv_coeff_; + delete reh1_coeff_; + delete reh2_coeff_; + delete Reh_coeff_; + delete csupg_coeff_; + delete uw1_coeff_; + delete uw2_coeff_; + delete upwind_coeff_; + delete swdiff_coeff_; + delete supg_coeff_; + // objects allocated by initalizeSelf if (axisym_) delete gravity_vec_; delete utheta_next_gf_; @@ -343,6 +363,8 @@ void Tomboulides::initializeSelf() { tmpR0_.SetSize(sfes_truevsize); tmpR1_.SetSize(vfes_truevsize); + swDiff_vec_.SetSize(vfes_truevsize); + tmpR0b_.SetSize(pfes_truevsize); gradU_.SetSize(vfes_truevsize); gradV_.SetSize(vfes_truevsize); @@ -661,6 +683,45 @@ void Tomboulides::initializeOperators() { swirl_var_viscosity_coeff_ = new InnerProductCoefficient(*grad_mu_coeff_, *utheta_vec_coeff_); } + // artifical diffusion coefficients + if(sw_stab_) { + umag_coeff_ = new VectorMagnitudeCoefficient(*u_next_coeff_); + gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_); + visc_inv_coeff_ = new PowerCoefficient(*mu_coeff_, -1.0); + + // compute Reh + reh1_coeff_ = new ProductCoefficient(*rho_coeff_, *visc_inv_coeff_); + reh2_coeff_ = new ProductCoefficient(*reh1_coeff_, *gscale_coeff_); + Reh_coeff_ = new ProductCoefficient(*reh2_coeff_, *umag_coeff_); + + // Csupg + csupg_coeff_ = new TransformedCoefficient(Reh_coeff_, csupgFactor); + + if (axisym_) { + // compute upwind magnitude + uw1_coeff_ = new ProductCoefficient(*rad_rho_coeff_, *csupg_coeff_); + uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_); + upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_); + + // streamwise diffusion direction + swdiff_coeff_ = new TransformedMatrixVectorCoefficient(u_next_coeff_, &streamwiseTensor); + + supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_); + } + else { + // compute upwind magnitude + // dividing by rho anyway + // uw1_coeff_ = new ProductCoefficient(*rho_coeff_, *csupg_coeff_); + uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_); + upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_); + + // streamwise diffusion direction + swdiff_coeff_ = new TransformedMatrixVectorCoefficient(u_next_coeff_, &streamwiseTensor); + + supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_); + } + } + // Integration rules (only used if numerical_integ_ is true). When // this is the case, the quadrature degree set such that the // Gauss-Lobatto quad pts correspond to the Gauss-Lobatto nodes. @@ -688,6 +749,18 @@ void Tomboulides::initializeOperators() { L_iorho_blfi->SetIntRule(&ir_ni_p); } L_iorho_form_->AddDomainIntegrator(L_iorho_blfi); + + // DiffusionIntegrator *slio_blfi; + // if (sw_stab_) { + // // auto *slio_blfi = new DiffusionIntegrator(*supg_coeff_); + // slio_blfi = new DiffusionIntegrator(*supg_coeff_); + // // SUPG diffusion + // // if (numerical_integ_) { + // // slio_blfi->SetIntRule(&ir_di); + // // } + // L_iorho_form_->AddDomainIntegrator(slio_blfi); + // } + L_iorho_form_->Assemble(); L_iorho_form_->FormSystemMatrix(pres_ess_tdof_, L_iorho_op_); @@ -809,6 +882,19 @@ void Tomboulides::initializeOperators() { Mv_rho_form_->Assemble(); Mv_rho_form_->FormSystemMatrix(empty, Mv_rho_op_); + // Vector mass matrix for streamwise stability of the velocity gradients + if (sw_stab_) { + Mv_stab_form_ = new ParBilinearForm(vfes_); + VectorMassIntegrator *mvs_blfi; + mvs_blfi = new VectorMassIntegrator(*supg_coeff_); + if (numerical_integ_) { + mvs_blfi->SetIntRule(&ir_ni_v); + } + Mv_stab_form_->AddDomainIntegrator(mvs_blfi); + Mv_stab_form_->Assemble(); + Mv_stab_form_->FormSystemMatrix(empty, Mv_stab_op_); + } + // Inverse (unweighted) mass operator (velocity space) if (partial_assembly_) { Vector diag_pa(vfes_->GetTrueVSize()); @@ -886,6 +972,7 @@ void Tomboulides::initializeOperators() { Hv_form_ = new ParBilinearForm(vfes_); VectorMassIntegrator *hmv_blfi; VectorDiffusionIntegrator *hdv_blfi; + if (axisym_) { hmv_blfi = new VectorMassIntegrator(*rad_rho_over_dt_coeff_); hdv_blfi = new VectorDiffusionIntegrator(*rad_mu_coeff_); @@ -899,6 +986,17 @@ void Tomboulides::initializeOperators() { } Hv_form_->AddDomainIntegrator(hmv_blfi); Hv_form_->AddDomainIntegrator(hdv_blfi); + + VectorDiffusionIntegrator *shdv_blfi; + if (sw_stab_) { + // auto *shdv_blfi = new VectorDiffusionIntegrator(*supg_coeff_); + shdv_blfi = new VectorDiffusionIntegrator(*supg_coeff_); + // if (numerical_integ_) { + // shdv_blfi->SetIntRule(&ir_di); + // } + Hv_form_->AddDomainIntegrator(shdv_blfi); + } + if (axisym_) { auto *hfv_blfi = new VectorMassIntegrator(*visc_forcing_coeff_); Hv_form_->AddDomainIntegrator(hfv_blfi); @@ -995,6 +1093,13 @@ void Tomboulides::initializeOperators() { auto *hms_blfi = new MassIntegrator(*rad_rho_over_dt_coeff_); auto *hds_blfi = new DiffusionIntegrator(*rad_mu_coeff_); auto *hfs_blfi = new MassIntegrator(*mu_over_rad_coeff_); + + DiffusionIntegrator *shds_blfi; + if (sw_stab_) { + // auto *shds_blfi = new DiffusionIntegrator(*supg_coeff_); + shds_blfi = new DiffusionIntegrator(*supg_coeff_); + Hs_form_->AddDomainIntegrator(shds_blfi); + } Hs_form_->AddDomainIntegrator(hms_blfi); Hs_form_->AddDomainIntegrator(hds_blfi); @@ -1434,6 +1539,26 @@ void Tomboulides::step() { }); } + if (sw_stab_) { + /* + for (int i = 0; i < dim_; i++) { + setScalarFromVector(u_vec_, i, &tmpR0a_); + streamwiseDiffusion(tmpR0a_, tmpR0b_); + setVectorFromScalar(tmpR0b_, i, &swDiff_vec_); + } + */ + Mv_stab_op_->Mult(gradU_, tmpR0b_); + setVectorFromScalar(tmpR0b_, 0, &swDiff_vec_); + Mv_stab_op_->Mult(gradV_, tmpR0b_); + setVectorFromScalar(tmpR0b_, 1, &swDiff_vec_); + if (dim_ == 3) { + Mv_stab_op_->Mult(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_; diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp index 1463b426e..4a1703ff3 100644 --- a/src/tomboulides.hpp +++ b/src/tomboulides.hpp @@ -181,6 +181,9 @@ class Tomboulides final : public FlowBase { // reference here. const temporalSchemeCoefficients &coeff_; + // streamwise-stabilization + bool sw_stab_; + StopWatch sw_press_, sw_helm_; std::string ic_string_; @@ -234,6 +237,8 @@ class Tomboulides final : public FlowBase { mfem::ParGridFunction *resp_gf_ = nullptr; mfem::ParGridFunction *pp_div_rad_comp_gf_ = nullptr; + mfem::ParGridFunction *gridScale_gf_ = nullptr; + /// Swirl mfem::ParGridFunction *utheta_gf_ = nullptr; mfem::ParGridFunction *utheta_next_gf_ = nullptr; @@ -286,6 +291,19 @@ class Tomboulides final : public FlowBase { mfem::VectorArrayCoefficient *utheta_vec_coeff_ = nullptr; mfem::InnerProductCoefficient *swirl_var_viscosity_coeff_ = nullptr; + mfem::VectorMagnitudeCoefficient *umag_coeff_ = nullptr; + mfem::GridFunctionCoefficient *gscale_coeff_ = nullptr; + mfem::PowerCoefficient *visc_inv_coeff_ = nullptr; + mfem::ProductCoefficient *reh1_coeff_ = nullptr; + mfem::ProductCoefficient *reh2_coeff_ = nullptr; + mfem::ProductCoefficient *Reh_coeff_ = nullptr; + mfem::TransformedCoefficient *csupg_coeff_ = nullptr; + mfem::ProductCoefficient *uw1_coeff_ = nullptr; + mfem::ProductCoefficient *uw2_coeff_ = nullptr; + mfem::ProductCoefficient *upwind_coeff_ = nullptr; + mfem::TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr; + mfem::ScalarMatrixProductCoefficient *supg_coeff_ = nullptr; + // mfem "form" objects used to create operators mfem::ParBilinearForm *L_iorho_form_ = nullptr; // \int (1/\rho) \nabla \phi_i \cdot \nabla \phi_j mfem::ParLinearForm *forcing_form_ = nullptr; // \int \phi_i f @@ -309,6 +327,9 @@ class Tomboulides final : public FlowBase { mfem::ParLinearForm *rho_ur_ut_form_ = nullptr; mfem::ParLinearForm *swirl_var_viscosity_form_ = nullptr; + // streamwise stability + mfem::ParBilinearForm *Mv_stab_form_ = nullptr; + // mfem operator objects mfem::OperatorHandle L_iorho_op_; mfem::OperatorHandle Ms_op_; @@ -321,6 +342,9 @@ class Tomboulides final : public FlowBase { mfem::OperatorHandle Hs_op_; mfem::OperatorHandle As_op_; + // streamwise stability + mfem::OperatorHandle Mv_stab_op_; + // solver objects mfem::ParLORDiscretization *L_iorho_lor_ = nullptr; mfem::HypreBoomerAMG *L_iorho_inv_pc_ = nullptr; @@ -374,6 +398,9 @@ class Tomboulides final : public FlowBase { mfem::Vector utheta_m2_vec_; mfem::Vector utheta_next_vec_; + mfem::Vector tmpR0b_; + mfem::Vector swDiff_vec_; + // miscellaneous double volume_; mfem::ParLinearForm *mass_lform_ = nullptr; @@ -392,7 +419,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, + mfem::ParGridFunction *gridScale, TPS::Tps *tps = nullptr); /// Destructor ~Tomboulides() final; From d8fa14192f3c5f42cbd4b7423eb556e7321daca3 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 18 Sep 2025 17:28:22 -0500 Subject: [PATCH 2/6] fixed Tomboulides bug (no help in torch) --- src/tomboulides.cpp | 15 ++++++++++++--- src/tomboulides.hpp | 1 + 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index 3902da648..edab3710a 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -210,7 +210,7 @@ Tomboulides::~Tomboulides() { delete rho_over_dt_coeff_; delete iorho_coeff_; delete rho_coeff_; - + delete umag_coeff_; delete gscale_coeff_; delete visc_inv_coeff_; @@ -223,6 +223,7 @@ Tomboulides::~Tomboulides() { delete upwind_coeff_; delete swdiff_coeff_; delete supg_coeff_; + delete visc_coeff_; // objects allocated by initalizeSelf if (axisym_) delete gravity_vec_; @@ -685,9 +686,11 @@ void Tomboulides::initializeOperators() { // artifical diffusion coefficients if(sw_stab_) { + visc_coeff_ = new GridFunctionCoefficient(thermo_interface_->viscosity); umag_coeff_ = new VectorMagnitudeCoefficient(*u_next_coeff_); gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_); - visc_inv_coeff_ = new PowerCoefficient(*mu_coeff_, -1.0); + visc_inv_coeff_ = new PowerCoefficient(*visc_coeff_, -1.0); + // visc_inv_coeff_ = new PowerCoefficient(*mu_coeff_, -1.0); // compute Reh reh1_coeff_ = new ProductCoefficient(*rho_coeff_, *visc_inv_coeff_); @@ -1540,6 +1543,12 @@ void Tomboulides::step() { } if (sw_stab_) { + + // Update matrix + Array empty; + Mv_stab_form_->Update(); + Mv_stab_form_->Assemble(); + Mv_stab_form_->FormSystemMatrix(empty, Mv_stab_op_); /* for (int i = 0; i < dim_; i++) { setScalarFromVector(u_vec_, i, &tmpR0a_); @@ -1558,7 +1567,7 @@ void Tomboulides::step() { Mv_rho_inv_->Mult(swDiff_vec_, tmpR1_); pp_div_vec_ += tmpR1_; } - + // printf("%f\n", tmpR1_.Norml2()); // Add ustar/dt contribution pp_div_vec_ += ustar_vec_; diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp index 4a1703ff3..2a4016097 100644 --- a/src/tomboulides.hpp +++ b/src/tomboulides.hpp @@ -303,6 +303,7 @@ class Tomboulides final : public FlowBase { mfem::ProductCoefficient *upwind_coeff_ = nullptr; mfem::TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr; mfem::ScalarMatrixProductCoefficient *supg_coeff_ = nullptr; + mfem::GridFunctionCoefficient *visc_coeff_ = nullptr; // mfem "form" objects used to create operators mfem::ParBilinearForm *L_iorho_form_ = nullptr; // \int (1/\rho) \nabla \phi_i \cdot \nabla \phi_j From ee1de4dd910af8fd637e4ffa1cdc9210a64ae5ef Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 18 Sep 2025 18:05:22 -0500 Subject: [PATCH 3/6] Minor change in Tomboulides --- src/tomboulides.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index edab3710a..090ea11a4 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -649,6 +649,8 @@ void Tomboulides::initializeOperators() { S_poisson_coeff_ = new VectorSumCoefficient(*twoS_gradmu_coeff_, *gradmu_Qt_coeff_, 1.0, -2. / 3); S_mom_coeff_ = new VectorSumCoefficient(*graduT_gradmu_coeff_, *gradmu_Qt_coeff_, 1.0, -1.0); + u_next_coeff_ = new VectorGridFunctionCoefficient(u_next_gf_); + // Coefficients for axisymmetric if (axisym_) { rad_rho_coeff_ = new ProductCoefficient(radius_coeff, *rho_coeff_); @@ -669,7 +671,6 @@ void Tomboulides::initializeOperators() { ur_conv_forcing_coeff_ = new VectorArrayCoefficient(2); ur_conv_forcing_coeff_->Set(0, utheta2_coeff_); - u_next_coeff_ = new VectorGridFunctionCoefficient(u_next_gf_); rad_rhou_coeff_ = new ScalarVectorProductCoefficient(*rad_rho_coeff_, *u_next_coeff_); u_next_rad_coeff_ = new GridFunctionCoefficient(u_next_rad_comp_gf_); @@ -689,6 +690,7 @@ void Tomboulides::initializeOperators() { visc_coeff_ = new GridFunctionCoefficient(thermo_interface_->viscosity); umag_coeff_ = new VectorMagnitudeCoefficient(*u_next_coeff_); gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_); + visc_inv_coeff_ = new PowerCoefficient(*visc_coeff_, -1.0); // visc_inv_coeff_ = new PowerCoefficient(*mu_coeff_, -1.0); From 795e4819cc0b5fc9f54d37ce93d5c1b1a4e1d84f Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 18 Sep 2025 17:57:23 -0500 Subject: [PATCH 4/6] Style fixes for Tomboulides --- src/loMach.cpp | 3 ++- src/tomboulides.cpp | 32 +++++++++++++++----------------- src/tomboulides.hpp | 18 +++++++++--------- 3 files changed, 26 insertions(+), 27 deletions(-) diff --git a/src/loMach.cpp b/src/loMach.cpp index 64ce5ca97..befa822cf 100644 --- a/src/loMach.cpp +++ b/src/loMach.cpp @@ -186,7 +186,8 @@ void LoMachSolver::initialize() { flow_ = new ZeroFlow(pmesh_, 1, tpsP_); } else if (loMach_opts_.flow_solver == "tomboulides") { // Tomboulides flow solver - flow_ = new Tomboulides(pmesh_, loMach_opts_.order, loMach_opts_.order, temporal_coeff_, (meshData_->getGridScale()), 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 090ea11a4..6da6ff20d 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -67,8 +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, - mfem::ParGridFunction *gridScale, TPS::Tps *tps) +Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, + mfem::ParGridFunction *gridScale, TPS::Tps *tps) : gll_rules(0, Quadrature1D::GaussLobatto), tpsP_(tps), pmesh_(pmesh), @@ -130,7 +130,7 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS tps->getInput("loMach/tomboulides/psolve-maxIters", pressure_solve_max_iter_, default_max_iter_); 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_); - + // artificial diffusion (SUPG) tpsP_->getInput("loMach/tomboulides/streamwise-stabilization", sw_stab_, false); } @@ -210,7 +210,7 @@ Tomboulides::~Tomboulides() { delete rho_over_dt_coeff_; delete iorho_coeff_; delete rho_coeff_; - + delete umag_coeff_; delete gscale_coeff_; delete visc_inv_coeff_; @@ -220,7 +220,7 @@ Tomboulides::~Tomboulides() { delete csupg_coeff_; delete uw1_coeff_; delete uw2_coeff_; - delete upwind_coeff_; + delete upwind_coeff_; delete swdiff_coeff_; delete supg_coeff_; delete visc_coeff_; @@ -650,7 +650,7 @@ void Tomboulides::initializeOperators() { S_mom_coeff_ = new VectorSumCoefficient(*graduT_gradmu_coeff_, *gradmu_Qt_coeff_, 1.0, -1.0); u_next_coeff_ = new VectorGridFunctionCoefficient(u_next_gf_); - + // Coefficients for axisymmetric if (axisym_) { rad_rho_coeff_ = new ProductCoefficient(radius_coeff, *rho_coeff_); @@ -686,34 +686,33 @@ void Tomboulides::initializeOperators() { } // artifical diffusion coefficients - if(sw_stab_) { + if (sw_stab_) { visc_coeff_ = new GridFunctionCoefficient(thermo_interface_->viscosity); umag_coeff_ = new VectorMagnitudeCoefficient(*u_next_coeff_); gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_); visc_inv_coeff_ = new PowerCoefficient(*visc_coeff_, -1.0); // visc_inv_coeff_ = new PowerCoefficient(*mu_coeff_, -1.0); - + // compute Reh reh1_coeff_ = new ProductCoefficient(*rho_coeff_, *visc_inv_coeff_); reh2_coeff_ = new ProductCoefficient(*reh1_coeff_, *gscale_coeff_); Reh_coeff_ = new ProductCoefficient(*reh2_coeff_, *umag_coeff_); - + // Csupg csupg_coeff_ = new TransformedCoefficient(Reh_coeff_, csupgFactor); - + if (axisym_) { // compute upwind magnitude uw1_coeff_ = new ProductCoefficient(*rad_rho_coeff_, *csupg_coeff_); uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_); upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_); - + // streamwise diffusion direction swdiff_coeff_ = new TransformedMatrixVectorCoefficient(u_next_coeff_, &streamwiseTensor); - + supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_); - } - else { + } else { // compute upwind magnitude // dividing by rho anyway // uw1_coeff_ = new ProductCoefficient(*rho_coeff_, *csupg_coeff_); @@ -977,7 +976,7 @@ void Tomboulides::initializeOperators() { Hv_form_ = new ParBilinearForm(vfes_); VectorMassIntegrator *hmv_blfi; VectorDiffusionIntegrator *hdv_blfi; - + if (axisym_) { hmv_blfi = new VectorMassIntegrator(*rad_rho_over_dt_coeff_); hdv_blfi = new VectorDiffusionIntegrator(*rad_mu_coeff_); @@ -1098,7 +1097,7 @@ void Tomboulides::initializeOperators() { auto *hms_blfi = new MassIntegrator(*rad_rho_over_dt_coeff_); auto *hds_blfi = new DiffusionIntegrator(*rad_mu_coeff_); auto *hfs_blfi = new MassIntegrator(*mu_over_rad_coeff_); - + DiffusionIntegrator *shds_blfi; if (sw_stab_) { // auto *shds_blfi = new DiffusionIntegrator(*supg_coeff_); @@ -1545,7 +1544,6 @@ void Tomboulides::step() { } if (sw_stab_) { - // Update matrix Array empty; Mv_stab_form_->Update(); diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp index 2a4016097..a508115d9 100644 --- a/src/tomboulides.hpp +++ b/src/tomboulides.hpp @@ -294,14 +294,14 @@ class Tomboulides final : public FlowBase { mfem::VectorMagnitudeCoefficient *umag_coeff_ = nullptr; mfem::GridFunctionCoefficient *gscale_coeff_ = nullptr; mfem::PowerCoefficient *visc_inv_coeff_ = nullptr; - mfem::ProductCoefficient *reh1_coeff_ = nullptr; - mfem::ProductCoefficient *reh2_coeff_ = nullptr; - mfem::ProductCoefficient *Reh_coeff_ = nullptr; + mfem::ProductCoefficient *reh1_coeff_ = nullptr; + mfem::ProductCoefficient *reh2_coeff_ = nullptr; + mfem::ProductCoefficient *Reh_coeff_ = nullptr; mfem::TransformedCoefficient *csupg_coeff_ = nullptr; - mfem::ProductCoefficient *uw1_coeff_ = nullptr; + mfem::ProductCoefficient *uw1_coeff_ = nullptr; mfem::ProductCoefficient *uw2_coeff_ = nullptr; - mfem::ProductCoefficient *upwind_coeff_ = nullptr; - mfem::TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr; + mfem::ProductCoefficient *upwind_coeff_ = nullptr; + mfem::TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr; mfem::ScalarMatrixProductCoefficient *supg_coeff_ = nullptr; mfem::GridFunctionCoefficient *visc_coeff_ = nullptr; @@ -328,7 +328,7 @@ class Tomboulides final : public FlowBase { mfem::ParLinearForm *rho_ur_ut_form_ = nullptr; mfem::ParLinearForm *swirl_var_viscosity_form_ = nullptr; - // streamwise stability + // streamwise stability mfem::ParBilinearForm *Mv_stab_form_ = nullptr; // mfem operator objects @@ -343,7 +343,7 @@ class Tomboulides final : public FlowBase { mfem::OperatorHandle Hs_op_; mfem::OperatorHandle As_op_; - // streamwise stability + // streamwise stability mfem::OperatorHandle Mv_stab_op_; // solver objects @@ -421,7 +421,7 @@ class Tomboulides final : public FlowBase { public: /// Constructor Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, - mfem::ParGridFunction *gridScale, TPS::Tps *tps = nullptr); + mfem::ParGridFunction *gridScale = nullptr, TPS::Tps *tps = nullptr); /// Destructor ~Tomboulides() final; From 5cf8bbbca4a2c9a23b55e48be7f5f83065a6a078 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 18 Sep 2025 21:38:34 -0500 Subject: [PATCH 5/6] fixed pressure artificial diffusion laplace operator --- src/tomboulides.cpp | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp index 6da6ff20d..53f9206a0 100644 --- a/src/tomboulides.cpp +++ b/src/tomboulides.cpp @@ -1549,19 +1549,16 @@ void Tomboulides::step() { Mv_stab_form_->Update(); Mv_stab_form_->Assemble(); Mv_stab_form_->FormSystemMatrix(empty, Mv_stab_op_); - /* - for (int i = 0; i < dim_; i++) { - setScalarFromVector(u_vec_, i, &tmpR0a_); - streamwiseDiffusion(tmpR0a_, tmpR0b_); - setVectorFromScalar(tmpR0b_, i, &swDiff_vec_); - } - */ - Mv_stab_op_->Mult(gradU_, tmpR0b_); + + Mv_stab_op_->Mult(gradU_, tmpR1_); + D_op_->Mult(tmpR1_, tmpR0b_); setVectorFromScalar(tmpR0b_, 0, &swDiff_vec_); - Mv_stab_op_->Mult(gradV_, tmpR0b_); + Mv_stab_op_->Mult(gradV_, tmpR1_); + D_op_->Mult(tmpR1_, tmpR0b_); setVectorFromScalar(tmpR0b_, 1, &swDiff_vec_); if (dim_ == 3) { - Mv_stab_op_->Mult(gradW_, tmpR0b_); + Mv_stab_op_->Mult(gradW_, tmpR1_); + D_op_->Mult(tmpR1_, tmpR0b_); setVectorFromScalar(tmpR0b_, 2, &swDiff_vec_); } Mv_rho_inv_->Mult(swDiff_vec_, tmpR1_); From 8cb944c2c7b125d882e593602c696601f7579b5a Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 18 Sep 2025 22:07:50 -0500 Subject: [PATCH 6/6] a stabilized test case (channel flow) --- .gitattributes | 1 + test/.gitattributes | 1 + test/Makefile.am | 9 ++- test/inputs/input.stabChan.ini | 77 +++++++++++++++++++ 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, 121 insertions(+), 3 deletions(-) create mode 100644 test/inputs/input.stabChan.ini create mode 100755 test/lomach-chan-stab.test create mode 100644 test/meshes/channel182p4_25x9p4_THIN.msh create mode 100644 test/ref_solns/stabChan/restart_output.sol.h5 diff --git a/.gitattributes b/.gitattributes index f98bb4510..84d62c6cd 100644 --- a/.gitattributes +++ b/.gitattributes @@ -30,3 +30,4 @@ test/ref_solns/aveLoMach/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -te test/ref_solns/reactNitrogen/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text test/ref_solns/supg/restart_output_wstab1T.sol.h5 filter=lfs diff=lfs merge=lfs -text test/ref_solns/supg/restart_output_lte_supg.sol.h5 filter=lfs diff=lfs merge=lfs -text +test/ref_solns/stabChan/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text diff --git a/test/.gitattributes b/test/.gitattributes index 83d6aca25..c68e89fe1 100644 --- a/test/.gitattributes +++ b/test/.gitattributes @@ -54,3 +54,4 @@ 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 ref_solns/collisions/Oliver.fit.h5 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 879ab4596..260bb6609 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -39,6 +39,7 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- ref_solns/radDecay/*.h5 \ ref_solns/reactNitrogen/*.h5 \ ref_solns/supg/*.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 \ @@ -52,10 +53,11 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- gradient.test coupled-3d.test coupled-3d.splitcomm.test coupled-3d.py-loop.test \ tabulated.test lte_mixture.test distance_fcn.test \ sgsSmag.test sgsSigma.test heatEq.test sponge.test plate.test pipe.mix.test lte2noneq-restart.test \ - coupled-3d.interface.test plasma.axisym.test plasma.axisym.lte1d.test \ + 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-table.test radiativeDecay.test \ - reactFlow-nitrogen.test argon_minimal.test argon_minimal.binary.test supg-scalars.test + reactFlow-nitrogen.test argon_minimal.test argon_minimal.binary.test supg-scalars.test \ + lomach-chan-stab.test TESTS = vpath.sh XFAIL_TESTS = @@ -130,7 +132,8 @@ TESTS += cyl3d.test \ reactFlow-table.test \ radiativeDecay.test \ reactFlow-nitrogen.test \ - supg-scalars.test + supg-scalars.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..cdfab8b28 --- /dev/null +++ b/test/inputs/input.stabChan.ini @@ -0,0 +1,77 @@ +[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 +streamwise-stabilization = true + +[loMach/tomboulides] +ic = channel +numerical-integ = false +streamwise-stabilization = true + +[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 \ No newline at end of file diff --git a/test/lomach-chan-stab.test b/test/lomach-chan-stab.test new file mode 100755 index 000000000..40798288c --- /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 +} \ No newline at end of file diff --git a/test/meshes/channel182p4_25x9p4_THIN.msh b/test/meshes/channel182p4_25x9p4_THIN.msh new file mode 100644 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