From 5d5ac55a235b270dfac18250ea8fad76e083b160 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 16 Jun 2024 11:22:11 -0700 Subject: [PATCH] fix corrector --- src/incflo_apply_corrector.cpp | 25 ++-- src/incflo_apply_predictor.cpp | 17 ++- src/incflo_update_tracer.cpp | 19 +++- src/incflo_update_velocity.cpp | 202 ++++++++++++++++++++------------- 4 files changed, 160 insertions(+), 103 deletions(-) diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index 2e10a4cd..f7627b3e 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -69,6 +69,8 @@ void incflo::ApplyCorrector() { BL_PROFILE("incflo::ApplyCorrector"); + constexpr Real m_half = Real(0.5); + // We use the new time value for things computed on the "*" state Real new_time = m_cur_time + m_dt; @@ -97,21 +99,19 @@ void incflo::ApplyCorrector() // in constructing the advection term // ********************************************************************************************** Vector vel_forces, tra_forces; - Vector vel_eta , tra_eta; - + Vector vel_eta, tra_eta; for (int lev = 0; lev <= finest_level; ++lev) { vel_forces.emplace_back(grids[lev], dmap[lev], AMREX_SPACEDIM, nghost_force(), MFInfo(), Factory(lev)); - vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); - } - - if (m_advect_tracer) { - for (int lev = 0; lev <= finest_level; ++lev) { + if (m_advect_tracer) { tra_forces.emplace_back(grids[lev], dmap[lev], m_ntrac, nghost_force(), MFInfo(), Factory(lev)); + } + vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); + if (m_advect_tracer) { tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev)); - } //lev - } // if + } + } // ************************************************************************************* // Compute the MAC-projected velocities at all levels @@ -135,10 +135,7 @@ void incflo::ApplyCorrector() // ************************************************************************************* // Compute viscosity / diffusive coefficients // ************************************************************************************* - compute_viscosity(GetVecOfPtrs(vel_eta), - get_density_new(), get_velocity_new(), - new_time, 1); - compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 1); // Here we create divtau of the (n+1,*) state that was computed in the predictor; // we use this laps only if DiffusionType::Explicit @@ -164,8 +161,8 @@ void incflo::ApplyCorrector() update_velocity(StepType::Corrector, vel_eta, vel_forces); // ********************************************************************************************** - // // Project velocity field, update pressure + // ********************************************************************************************** bool incremental_projection = false; ApplyProjection(get_density_nph_const(), new_time,m_dt,incremental_projection); diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index b7306f20..6ea2e23c 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -119,17 +119,12 @@ void incflo::ApplyPredictor (bool incremental_projection) } } - // ************************************************************************************* - // We now define the forcing terms to use in the Godunov prediction inside the predictor - // ************************************************************************************* - // ************************************************************************************* // Compute viscosity / diffusive coefficients // ************************************************************************************* compute_viscosity(GetVecOfPtrs(vel_eta), get_density_old(), get_velocity_old(), m_cur_time, 1); - compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); // ************************************************************************************* // Compute explicit viscous term @@ -142,12 +137,15 @@ void incflo::ApplyPredictor (bool incremental_projection) } // ************************************************************************************* - // Compute explicit diffusive terms + // Compute explicit diffusive term -- note this is used inside compute_convective_term // ************************************************************************************* - if (m_advect_tracer && need_divtau()) + if (m_advect_tracer) { - compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(), - GetVecOfConstPtrs(tra_eta)); + compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + if (need_divtau()) { + compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(), + GetVecOfConstPtrs(tra_eta)); + } } // ********************************************************************************************** @@ -170,6 +168,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // Compute the explicit advective terms R_u^(n+1/2), R_s^(n+1/2) and R_t^(n+1/2) // if (advection_type == "MOL" ) // Compute the explicit advective terms R_u^n , R_s^n and R_t^n + // Note that if advection_type != "MOL" then we call compute_tra_forces inside this routine // ************************************************************************************* compute_convective_term(get_conv_velocity_old(), get_conv_density_old(), get_conv_tracer_old(), get_velocity_old_const(), get_density_old_const(), get_tracer_old_const(), diff --git a/src/incflo_update_tracer.cpp b/src/incflo_update_tracer.cpp index 3c602b34..f59250ca 100644 --- a/src/incflo_update_tracer.cpp +++ b/src/incflo_update_tracer.cpp @@ -10,6 +10,10 @@ void incflo::update_tracer (StepType step_type, Vector& tra_eta, Vecto if (m_advect_tracer) { + // ************************************************************************************* + // Compute diffusion coefficient for tracer + // ************************************************************************************* + // ************************************************************************************* // Compute the tracer forcing terms (forcing for (rho s), not for s) // ************************************************************************************* @@ -18,11 +22,20 @@ void incflo::update_tracer (StepType step_type, Vector& tra_eta, Vecto // ************************************************************************************* // Compute explicit diffusive term (if corrector) // ************************************************************************************* - if (step_type == StepType::Corrector && m_diff_type == DiffusionType::Explicit) { - compute_laps(get_laps_new(), get_tracer_new_const(), get_density_new_const(), - GetVecOfConstPtrs(tra_eta)); + if (step_type == StepType::Corrector) + { + compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + if (m_diff_type == DiffusionType::Explicit) { + compute_laps(get_laps_new(), get_tracer_new_const(), get_density_new_const(), + GetVecOfConstPtrs(tra_eta)); + } } + // ************************************************************************************* + // Update the tracer next (note that dtdt already has rho in it) + // (rho trac)^new = (rho trac)^old + dt * ( + // div(rho trac u) + div (mu grad trac) + rho * f_t + // ************************************************************************************* if (step_type == StepType::Predictor) { tracer_explicit_update(tra_forces); } else if (step_type == StepType::Corrector) { diff --git a/src/incflo_update_velocity.cpp b/src/incflo_update_velocity.cpp index 48fbd183..7967c3f7 100644 --- a/src/incflo_update_velocity.cpp +++ b/src/incflo_update_velocity.cpp @@ -159,80 +159,40 @@ void incflo::update_velocity (StepType step_type, Vector& vel_eta, Vec #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { + for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { Box const& bx = mfi.tilebox(); Array4 const& vel = ld.velocity.array(mfi); - Array4 const& dvdt = ld.conv_velocity_o.const_array(mfi); + Array4 const& vel_o = ld.velocity_o.const_array(mfi); + Array4 const& dvdt = ld.conv_velocity.const_array(mfi); + Array4 const& dvdt_o = ld.conv_velocity_o.const_array(mfi); Array4 const& vel_f = vel_forces[lev].const_array(mfi); + Array4 const& rho_old = ld.density_o.const_array(mfi); Array4 const& rho_new = ld.density.const_array(mfi); Array4 const& rho_nph = ld.density_nph.const_array(mfi); - if (m_diff_type == DiffusionType::Implicit) { - - if (use_tensor_correction) - { - Array4 const& divtau_o = ld.divtau_o.const_array(mfi); - // Here divtau_o is the difference of tensor and scalar divtau_o! - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); - }); - } - } else { - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2));); - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2));); - }); - } - } - } - else if (m_diff_type == DiffusionType::Crank_Nicolson) + if (m_diff_type == DiffusionType::Explicit) { - Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + Array4 const& divtau = ld.divtau.const_array(mfi); + if (m_advect_momentum) { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+l_half*divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+l_half*divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+l_half*divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( + l_half*( dvdt_o(i,j,k,0)+ dvdt(i,j,k,0)) + + l_half*(divtau_o(i,j,k,0)+divtau(i,j,k,0)) + + rho_nph(i,j,k) * vel_f(i,j,k,0) );, + vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( + l_half*( dvdt_o(i,j,k,1)+ dvdt(i,j,k,1)) + + l_half*(divtau_o(i,j,k,1)+divtau(i,j,k,1)) + + rho_nph(i,j,k) * vel_f(i,j,k,1) );, + vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( + l_half*( dvdt_o(i,j,k,2)+ dvdt(i,j,k,2)) + + l_half*(divtau_o(i,j,k,2)+divtau(i,j,k,2)) + + rho_nph(i,j,k) * vel_f(i,j,k,2) );); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, vel(i,j,k,1) /= rho_new(i,j,k);, vel(i,j,k,2) /= rho_new(i,j,k);); @@ -240,24 +200,38 @@ void incflo::update_velocity (StepType step_type, Vector& vel_eta, Vec } else { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+l_half*divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+l_half*divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+l_half*divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( + l_half*( dvdt_o(i,j,k,0)+ dvdt(i,j,k,0)) + + l_half*(divtau_o(i,j,k,0)+divtau(i,j,k,0)) + + vel_f(i,j,k,0) );, + vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( + l_half*( dvdt_o(i,j,k,1)+ dvdt(i,j,k,1)) + + l_half*(divtau_o(i,j,k,1)+divtau(i,j,k,1)) + + vel_f(i,j,k,1) );, + vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( + l_half*( dvdt_o(i,j,k,2)+ dvdt(i,j,k,2)) + + l_half*(divtau_o(i,j,k,2)+divtau(i,j,k,2)) + + vel_f(i,j,k,2) );); }); } } - else if (m_diff_type == DiffusionType::Explicit) + else if (m_diff_type == DiffusionType::Crank_Nicolson) { Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + if (m_advect_momentum) { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( + l_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + +l_half*(divtau_o(i,j,k,0) ) + rho_nph(i,j,k)*vel_f(i,j,k,0) );, + vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( + l_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) + +l_half*(divtau_o(i,j,k,1) ) + rho_nph(i,j,k)*vel_f(i,j,k,1) );, + vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( + l_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + +l_half*(divtau_o(i,j,k,2) ) + rho_nph(i,j,k)*vel_f(i,j,k,2) );); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, vel(i,j,k,1) /= rho_new(i,j,k);, vel(i,j,k,2) /= rho_new(i,j,k);); @@ -265,15 +239,89 @@ void incflo::update_velocity (StepType step_type, Vector& vel_eta, Vec } else { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( + l_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + + l_half* divtau_o(i,j,k,0) + vel_f(i,j,k,0) );, + vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( + l_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) + + l_half* divtau_o(i,j,k,1) + vel_f(i,j,k,1) );, + vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( + l_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + + l_half* divtau_o(i,j,k,2) + vel_f(i,j,k,2) );); }); } } - } // mfi - } // lev + else if (m_diff_type == DiffusionType::Implicit) + { + if (use_tensor_correction) + { + Array4 const& divtau = ld.divtau.const_array(mfi); + if (m_advect_momentum) + { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( + l_half*(dvdt_o(i,j,k,0) + dvdt(i,j,k,0)) + + rho_nph(i,j,k) * vel_f(i,j,k,0) + divtau(i,j,k,0));, + vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( + l_half*(dvdt_o(i,j,k,1) + dvdt(i,j,k,1)) + + rho_nph(i,j,k) * vel_f(i,j,k,1) + divtau(i,j,k,1));, + vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( + l_half*(dvdt_o(i,j,k,2) + dvdt(i,j,k,2)) + + rho_nph(i,j,k) * vel_f(i,j,k,2) + divtau(i,j,k,2));); + if (i == 20 and j == 36) amrex::Print() <<" MOM " << vel(i,j,k,0) << std::endl; + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( + l_half*( dvdt_o(i,j,k,0) + dvdt(i,j,k,0)) + + vel_f(i,j,k,0) + divtau(i,j,k,0));, + vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( + l_half*(dvdt_o(i,j,k,1) + dvdt(i,j,k,1)) + + vel_f(i,j,k,1) + divtau(i,j,k,1) );, + vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( + l_half*( dvdt_o(i,j,k,2) + dvdt(i,j,k,2)) + + vel_f(i,j,k,2) + divtau(i,j,k,2) );); + }); + } + } else { + if (m_advect_momentum) + { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( + l_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + + rho_nph(i,j,k) * vel_f(i,j,k,0) );, + vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( + l_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) + + rho_nph(i,j,k) * vel_f(i,j,k,1) );, + vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( + l_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + + rho_nph(i,j,k) * vel_f(i,j,k,2) );); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( + l_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + vel_f(i,j,k,0) );, + vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( + l_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) + vel_f(i,j,k,1) );, + vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( + l_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + vel_f(i,j,k,2) );); + }); + } + } + } + } + } // lev } // Corrector // *************************************************************************************