diff --git a/src/scheme/LocalDtAcousticSolver.cpp b/src/scheme/LocalDtAcousticSolver.cpp index 256cd28506ced0e5a09049a561bfb03d148f3548..f1b43366318042fad3da83a8928352c5078863b5 100644 --- a/src/scheme/LocalDtAcousticSolver.cpp +++ b/src/scheme/LocalDtAcousticSolver.cpp @@ -767,37 +767,31 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2); auto [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, rho2, c2, u2, p2, bc_descriptor_list2,map1,map2); - auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1); + const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1); double dt2 = dt1/q; const double& gamma = 1.4; double sum_dt = 0; - std::shared_ptr<const DiscreteFunctionVariant> new_p2; + do{ - while(sum_dt != dt1){ + const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>(); + const DiscreteVectorFunction& u_d = new_u2->get<DiscreteVectorFunction>(); + const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>(); + const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d); + const DiscreteScalarFunction& p_d = (gamma-1) * rho_d * eps; + const DiscreteScalarFunction& c_d = sqrt(gamma * p_d / rho_d); - if(sum_dt + dt2 > dt1){ - dt2 = dt1 - sum_dt; - } - - std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2); - - DiscreteScalarFunction rho_d = new_rho2->get<DiscreteScalarFunction>(); - DiscreteVectorFunction u_d = new_u2->get<DiscreteVectorFunction>(); - DiscreteScalarFunction E_d = new_E2->get<DiscreteScalarFunction>(); - DiscreteScalarFunction eps = E_d - 0.5 * dot(u_d,u_d); - DiscreteScalarFunction p_d = (gamma-1) * rho_d * eps; - DiscreteScalarFunction c_d = sqrt(gamma * p_d / rho_d); - - new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d); - std::shared_ptr<const DiscreteFunctionVariant> new_c2 = std::make_shared<DiscreteFunctionVariant>(c_d); + const std::shared_ptr<const DiscreteFunctionVariant>& new_p2 = std::make_shared<const DiscreteFunctionVariant>(p_d); + const std::shared_ptr<const DiscreteFunctionVariant>& new_c2 = std::make_shared<const DiscreteFunctionVariant>(c_d); auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2); + std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2); + sum_dt += dt2; - } - + }while(sum_dt < dt1); + std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2); return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};