From acdb51c473b1360047f378006fe790470a9f3897 Mon Sep 17 00:00:00 2001
From: Alexandre Gangloff <alexandre.gangloff@cea.fr>
Date: Mon, 24 Feb 2025 14:06:23 +0100
Subject: [PATCH] Add time second order extension for coupling method
---
.../Order2LocalDtHyperelasticSolver.cpp | 105 ++++++++++++------
1 file changed, 72 insertions(+), 33 deletions(-)
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index a98d525a1..0dfe808d2 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -1911,6 +1911,67 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
}
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ compute_sub_iterations_midpoint_method(const SolverType& solver_type,
+ const size_t& law,
+ const double& CFL,
+ const double& dt,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr,
+ std::vector<NodeId> map2,
+ const size_t& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma,
+ const double& p_inf) const
+ {
+ std::shared_ptr sub_m = getCommonMesh({rho, u});
+ std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_u = u;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_E = E;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_CG = CG;
+
+ double sum_dt = 0;
+ while(sum_dt < dt){
+ const auto& [sigma, aL, aT, p] = _apply_state_law(law,sub_rho,sub_u,sub_E,sub_CG,chi_solid,lambda,mu,gamma,p_inf);
+
+ const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
+ const std::shared_ptr<const DiscreteFunctionVariant>& c =
+ std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+
+ double sub_dt = CFL * hyperelastic_dt(c);
+ if(sum_dt + sub_dt > dt){
+ sub_dt = dt - sum_dt;
+ }
+
+ auto [sub_ur, sub_Fjr] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+ auto [sub_ur_o1, sub_Fjr_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+
+ auto [sub_m_app, sub_rho_app, sub_u_app, sub_E_app, sub_CG_app] = apply_fluxes(sub_dt/2., sub_rho, sub_u, sub_E, sub_CG, sub_ur, sub_ur_o1, sub_Fjr, sub_Fjr_o1);
+
+ const auto& [sigma_app, aL_app, aT_app, p_app] = _apply_state_law(law,sub_rho_app,sub_u_app,sub_E_app,sub_CG_app,chi_solid,lambda,mu,gamma,p_inf);
+
+ auto [sub_ur_app, sub_Fjr_app] = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, p_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+ auto [sub_ur_o1_app, sub_Fjr_o1_app] = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+
+ std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = order2_apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_rho_app, sub_CG_app, sub_ur_app, sub_ur_o1_app, sub_Fjr_app, sub_Fjr_o1_app);
+
+ sum_dt += sub_dt;
+ }
+
+ return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
+ }
+
//std::tuple<std::shared_ptr<const MeshVariant>,
// std::shared_ptr<const DiscreteFunctionVariant>,
// std::shared_ptr<const DiscreteFunctionVariant>,
@@ -2050,51 +2111,29 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
- //std::shared_ptr m1_app = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
- //std::shared_ptr<const DiscreteFunctionVariant> rho1_app = rho1;
- //std::shared_ptr<const DiscreteFunctionVariant> u1_app = u1;
- //std::shared_ptr<const DiscreteFunctionVariant> E1_app = E1;
- //std::shared_ptr<const DiscreteFunctionVariant> CG1_app = CG1;
-
- //std::shared_ptr m2_app = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
- //std::shared_ptr<const DiscreteFunctionVariant> rho2_app = rho2;
- //std::shared_ptr<const DiscreteFunctionVariant> u2_app = u2;
- //std::shared_ptr<const DiscreteFunctionVariant> E2_app = E2;
- //std::shared_ptr<const DiscreteFunctionVariant> CG2_app = CG2;
- std::shared_ptr<const DiscreteFunctionVariant> aL2_app = aL2;
- std::shared_ptr<const DiscreteFunctionVariant> aT2_app = aT2;
- std::shared_ptr<const DiscreteFunctionVariant> sigma2_app = sigma2;
- std::shared_ptr<const DiscreteFunctionVariant> p2_app = p2;
-
-
auto [map1, map2] = _computeMapping(m1,m2,bc_descriptor_list1,bc_descriptor_list2);
- auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,p1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,p2,bc_descriptor_list2,map1,map2);
- auto [ur1_o1, Fjr1_o1, ur2_o1, Fjr2_o1, CR_ur_o1, CR_Fjr_o1] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,bc_descriptor_list2,map1,map2);
-
- //auto [m1_app, rho1_app, u1_app, E1_app, CG1_app] = apply_fluxes(dt1/2., rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);
-
+ //Préparation au calcul des lois d'états
size_t chi_solid = 0;
if((mu_s!=0) or (lambda_s!=0)){
chi_solid = 1;
}
const size_t chi_fluid = 1-chi_solid;
- auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = compute_sub_iterations(solver_type,law,0.04,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+ auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,p1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,p2,bc_descriptor_list2,map1,map2);
+ auto [ur1_o1, Fjr1_o1, ur2_o1, Fjr2_o1, CR_ur_o1, CR_Fjr_o1] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,bc_descriptor_list2,map1,map2);
- //const auto& [sigma1_app, aL1_app, aT1_app, p1_app] = _apply_state_law(law,rho1_app,u1_app,E1_app,CG1_app,chi_fluid,lambda_f,mu_f,gamma_f,p_inf_f);
- //std::tie(sigma2_app, aL2_app, aT2_app, p2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+ auto [m1_app, rho1_app, u1_app, E1_app, CG1_app] = apply_fluxes(dt1/2., rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);
+ const auto& [sigma1_app, aL1_app, aT1_app, p1_app] = _apply_state_law(law, rho1_app, u1_app, E1_app, CG1_app, chi_fluid, lambda_f, mu_f, gamma_f, p_inf_f);
- //auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app] = compute_fluxes(solver_type,rho1_app, aL1_app, aT1_app, u1_app, sigma1_app, p1_app, bc_descriptor_list1,
- // rho2_app, aL2_app, aT2_app, u2_app, sigma2_app, p2_app, bc_descriptor_list2, map1, map2);
- //auto [ur1_app_o1, Fjr1_app_o1, ur2_app_o1, Fjr2_app_o1, CR_ur_app_o1, CR_Fjr_app_o1] = compute_fluxes(solver_type,rho1_app, aL1_app, aT1_app, u1_app, sigma1_app, bc_descriptor_list1,
- // rho2_app, aL2_app, aT2_app, u2_app, sigma2_app, bc_descriptor_list2, map1, map2);
+ auto [m2_app, rho2_app, u2_app, E2_app, CG2_app] = compute_sub_iterations(solver_type,law,0.03,dt1/2.,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+ const auto& [sigma2_app, aL2_app, aT2_app, p2_app] = _apply_state_law(law, rho2_app, u2_app, E2_app, CG2_app, chi_solid, lambda_s, mu_s, gamma_s, p_inf_s);
- //const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = order2_apply_fluxes(dt1, rho1, u1, E1, CG1, rho1_app, CG1_app, ur1_app, ur1_app_o1, Fjr1_app, Fjr1_app_o1);
- //auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,law,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,p2,bc_descriptor_list2,
- // CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+ auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app] = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,p1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,p2_app,bc_descriptor_list2,map1,map2);
+ auto [ur1_o1_app, Fjr1_o1_app, ur2_o1_app, Fjr2_o1_app, CR_ur_o1_app, CR_Fjr_o1_app] = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,bc_descriptor_list2,map1,map2);
- const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);
+ const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = order2_apply_fluxes(dt1, rho1, u1, E1, CG1, rho1_app, CG1_app, ur1_app, ur1_o1_app, Fjr1_app, Fjr1_o1_app);
+ auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = compute_sub_iterations_midpoint_method(solver_type,law,0.03,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) = mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
--
GitLab