diff --git a/src/incflo.H b/src/incflo.H index b9accecc..120c086e 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -163,7 +163,8 @@ public: void tracer_explicit_update_corrector(amrex::Vector const& tra_forces); void update_density (StepType step_type); - void update_tracer (StepType step_type, amrex::Vector& tra_forces); + void update_tracer (StepType step_type, amrex::Vector& tra_eta, + amrex::Vector& tra_forces); void update_velocity (StepType step_type, amrex::Vector& vel_eta, amrex::Vector& vel_forces); diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index 27df84ef..2e10a4cd 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -97,17 +97,22 @@ void incflo::ApplyCorrector() // in constructing the advection term // ********************************************************************************************** Vector vel_forces, tra_forces; - Vector vel_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)); - 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) { + for (int lev = 0; lev <= finest_level; ++lev) { + tra_forces.emplace_back(grids[lev], dmap[lev], m_ntrac, nghost_force(), + MFInfo(), Factory(lev)); + tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev)); + } //lev + } // if + // ************************************************************************************* // Compute the MAC-projected velocities at all levels // ************************************************************************************* @@ -133,6 +138,7 @@ void incflo::ApplyCorrector() compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 1); + compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); // Here we create divtau of the (n+1,*) state that was computed in the predictor; // we use this laps only if DiffusionType::Explicit @@ -143,14 +149,14 @@ void incflo::ApplyCorrector() } // ************************************************************************************* - // Update density first + // Update density // ************************************************************************************* update_density(StepType::Corrector); // ************************************************************************************* - // Update tracer next + // Update tracer // ************************************************************************************* - update_tracer(StepType::Corrector, tra_forces); + update_tracer(StepType::Corrector, tra_eta, tra_forces); // ************************************************************************************* // Update velocity diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index ebbcc734..b7306f20 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -94,10 +94,13 @@ void incflo::ApplyPredictor (bool incremental_projection) } } + // ************************************************************************************* + // Allocate space for half-time density + // ************************************************************************************* // Forcing terms Vector vel_forces, tra_forces; - Vector vel_eta; + Vector vel_eta, tra_eta; // ************************************************************************************* // Allocate space for the forcing terms @@ -111,6 +114,9 @@ void incflo::ApplyPredictor (bool incremental_projection) 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)); + } } // ************************************************************************************* @@ -123,23 +129,38 @@ void incflo::ApplyPredictor (bool incremental_projection) 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 // Note that for !advect_momentum, this actually computes divtau / rho // ************************************************************************************* - if (need_divtau() || use_tensor_correction ) { + if (need_divtau() || use_tensor_correction ) + { compute_divtau(get_divtau_old(),get_velocity_old_const(), get_density_old_const(),GetVecOfConstPtrs(vel_eta)); } + // ************************************************************************************* + // Compute explicit diffusive terms + // ************************************************************************************* + if (m_advect_tracer && need_divtau()) + { + compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(), + GetVecOfConstPtrs(tra_eta)); + } + // ********************************************************************************************** - // Compute the MAC-projected velocities at all levels + // Compute the forcing terms // ************************************************************************************* bool include_pressure_gradient = !(m_use_mac_phi_in_godunov); compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_old_const(), get_density_old_const(), get_tracer_old_const(), get_tracer_old_const(), include_pressure_gradient); + + // ********************************************************************************************** + // Compute the MAC-projected velocities at all levels + // ************************************************************************************* compute_MAC_projected_velocities(get_velocity_old_const(), get_density_old_const(), AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), GetVecOfPtrs(w_mac)), GetVecOfPtrs(vel_forces), m_cur_time); @@ -162,20 +183,18 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* update_density(StepType::Predictor); - // ************************************************************************************* + // ********************************************************************************************** // Update tracer - // ************************************************************************************* - update_tracer(StepType::Predictor, tra_forces); + // ********************************************************************************************** + update_tracer(StepType::Predictor, tra_eta, tra_forces); - // ************************************************************************************* + // ********************************************************************************************** // Update velocity - // ************************************************************************************* + // ********************************************************************************************** update_velocity(StepType::Predictor, vel_eta, vel_forces); // ********************************************************************************************** - // // Project velocity field, update pressure - // // ********************************************************************************************** ApplyProjection(get_density_nph_const(),new_time,m_dt,incremental_projection); @@ -191,9 +210,7 @@ void incflo::ApplyPredictor (bool incremental_projection) #ifdef AMREX_USE_EB // ********************************************************************************************** - // // Over-write velocity in cells with vfrac < 1e-4 - // // ********************************************************************************************** if (m_advection_type == "MOL") incflo_correct_small_cells(get_velocity_new(), diff --git a/src/incflo_update_tracer.cpp b/src/incflo_update_tracer.cpp index a3a3e8b0..3c602b34 100644 --- a/src/incflo_update_tracer.cpp +++ b/src/incflo_update_tracer.cpp @@ -2,32 +2,23 @@ using namespace amrex; -void incflo::update_tracer (StepType step_type, Vector& tra_forces) +void incflo::update_tracer (StepType step_type, Vector& tra_eta, Vector& tra_forces) { BL_PROFILE("incflo::update_tracer"); - Vector tra_eta; - Real new_time = m_cur_time + m_dt; if (m_advect_tracer) { // ************************************************************************************* - // Compute the tracer forcing terms ( forcing for (rho s) if conservative ) + // Compute the tracer forcing terms (forcing for (rho s), not for s) // ************************************************************************************* - for (int lev = 0; lev <= finest_level; ++lev) { - tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev)); - } - - compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + compute_tra_forces(GetVecOfPtrs(tra_forces), get_density_nph_const()); // ************************************************************************************* - // Compute explicit diffusive terms + // Compute explicit diffusive term (if corrector) // ************************************************************************************* - if (step_type == StepType::Predictor && need_divtau()) { - compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(), - GetVecOfConstPtrs(tra_eta)); - } else if (step_type == StepType::Corrector && m_diff_type == DiffusionType::Explicit) { + 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)); }