Merge branch 'development' into add_init_pressure_proj

cgilet · Jun 17, 2024 · f7fd6cb · f7fd6cb
2 parents 6f4c47d + bc553d1
commit f7fd6cb
Show file tree

Hide file tree

Showing 15 changed files with 947 additions and 973 deletions.
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -11,8 +11,12 @@ target_sources(incflo
    incflo_correct_small_cells.cpp
    incflo_redistribute.cpp
    incflo_explicit_update.cpp
-   incflo_tagging.cpp
    incflo_regrid.cpp
+   incflo_tagging.cpp
+   incflo_update_density.cpp
+   incflo_update_tracer.cpp
+   incflo_update_velocity.cpp
+   incflo_utils.cpp
    main.cpp
    )
 

diff --git a/src/Make.package b/src/Make.package
@@ -5,8 +5,12 @@ CEXE_sources += incflo_apply_corrector.cpp
 CEXE_sources += incflo_compute_dt.cpp
 CEXE_sources += incflo_compute_forces.cpp
 CEXE_sources += incflo_explicit_update.cpp
-CEXE_sources += incflo_tagging.cpp
 CEXE_sources += incflo_regrid.cpp
+CEXE_sources += incflo_tagging.cpp
+CEXE_sources += incflo_update_density.cpp
+CEXE_sources += incflo_update_tracer.cpp
+CEXE_sources += incflo_update_velocity.cpp
+CEXE_sources += incflo_utils.cpp
 CEXE_sources += main.cpp
 
 ifeq ($(USE_EB), TRUE)

diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp
@@ -107,21 +107,24 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 
     for (int lev = 0; lev <= finest_level; ++lev) {
         AMREX_D_TERM(
-           face_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));,
-           face_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));,
-           face_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev)););
+           face_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),u_mac[lev]->Factory());,
+           face_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),v_mac[lev]->Factory());,
+           face_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),w_mac[lev]->Factory()););
         AMREX_D_TERM(
-           flux_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));,
-           flux_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));,
-           flux_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev)););
+           flux_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),u_mac[lev]->Factory());,
+           flux_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),v_mac[lev]->Factory());,
+           flux_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),w_mac[lev]->Factory()););
 
-        divu[lev].define(vel[lev]->boxArray(),dmap[lev],1,4,MFInfo(),Factory(lev));
-        if (m_advect_momentum)
+        divu[lev].define(vel[lev]->boxArray(),dmap[lev],1,4,MFInfo(),vel[lev]->Factory());
+
+        if (m_advect_momentum) {
             rhovel[lev].define(vel[lev]->boxArray(),dmap[lev],AMREX_SPACEDIM,
-                               vel[lev]->nGrow(),MFInfo(),Factory(lev));
-        if (m_advect_tracer && m_ntrac > 0 && any_conserv_trac)
-            rhotrac[lev].define(vel[lev]->boxArray(),dmap[lev],tracer[lev]->nComp(),
-                                tracer[lev]->nGrow(),MFInfo(),Factory(lev));
+                               vel[lev]->nGrow(),MFInfo(),vel[lev]->Factory());
+        }
+        if (m_advect_tracer && m_ntrac > 0 && any_conserv_trac) {
+            rhotrac[lev].define(tracer[lev]->boxArray(),dmap[lev],tracer[lev]->nComp(),
+                                tracer[lev]->nGrow(),MFInfo(),tracer[lev]->Factory());
+        }
 
         AMREX_D_TERM(faces[lev][0] = &face_x[lev];,
                      faces[lev][1] = &face_y[lev];,

diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H
@@ -28,7 +28,6 @@ public:
 
     void compute_laps (amrex::Vector<amrex::MultiFab*> const& laps,
                        amrex::Vector<amrex::MultiFab const*> const& a_scalar,
-                       amrex::Vector<amrex::MultiFab const*> const& /*density*/,
                        amrex::Vector<amrex::MultiFab const*> const& eta);
 
     void compute_divtau (amrex::Vector<amrex::MultiFab*> const& a_divtau,

diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp
@@ -474,7 +474,6 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
 
 void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
                                       Vector<MultiFab const*> const& a_scalar,
-                                      Vector<MultiFab const*> const& /*a_density*/,
                                       Vector<MultiFab const*> const& a_eta)
 {
     BL_PROFILE("DiffusionScalarOp::compute_laps");

diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp
@@ -57,10 +57,9 @@ incflo::compute_divtau(Vector<MultiFab      *> const& divtau,
 void
 incflo::compute_laps(Vector<MultiFab      *> const& laps,
                      Vector<MultiFab const*> const& scalar,
-                     Vector<MultiFab const*> const& density,
                      Vector<MultiFab const*> const& eta)
 {
-    get_diffusion_scalar_op()->compute_laps(laps, scalar, density, eta);
+    get_diffusion_scalar_op()->compute_laps(laps, scalar, eta);
 }
 
 void

diff --git a/src/incflo.H b/src/incflo.H
@@ -21,6 +21,10 @@
 #include <DiffusionTensorOp.H>
 #include <DiffusionScalarOp.H>
 
+enum struct StepType {
+    Predictor, Corrector
+};
+
 class incflo : public amrex::AmrCore
 {
 public:
@@ -131,8 +135,8 @@ public:
     //
     ///////////////////////////////////////////////////////////////////////////
 
-    void ApplyPredictor(bool incremental_projection = false);
-    void ApplyCorrector();
+    void ApplyPredictor (bool incremental_projection = false);
+    void ApplyCorrector ();
     void compute_convective_term (amrex::Vector<amrex::MultiFab*> const& conv_u,
                                   amrex::Vector<amrex::MultiFab*> const& conv_r,
                                   amrex::Vector<amrex::MultiFab*> const& conv_t,
@@ -158,6 +162,12 @@ public:
     void tracer_explicit_update(amrex::Vector<amrex::MultiFab> const& tra_forces);
     void tracer_explicit_update_corrector(amrex::Vector<amrex::MultiFab> const& tra_forces);
 
+    void update_density  (StepType step_type);
+    void update_tracer   (StepType step_type, amrex::Vector<amrex::MultiFab>& tra_eta,
+                                              amrex::Vector<amrex::MultiFab>& tra_forces);
+    void update_velocity (StepType step_type, amrex::Vector<amrex::MultiFab>& vel_eta,
+                                              amrex::Vector<amrex::MultiFab>& vel_forces);
+
     ///////////////////////////////////////////////////////////////////////////
     //
     // derive
@@ -201,11 +211,8 @@ public:
     [[nodiscard]] amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>
     average_scalar_eta_to_faces (int lev, int comp, amrex::MultiFab const& cc_eta) const;
 
-    // FIXME? compute_laps never used density in any meaningful way, but
-    // then removing density would break anything downstream....
     void compute_laps (amrex::Vector<amrex::MultiFab      *> const& laps,
                        amrex::Vector<amrex::MultiFab const*> const& scalar,
-                       amrex::Vector<amrex::MultiFab const*> const& density,
                        amrex::Vector<amrex::MultiFab const*> const& eta);
 
     void diffuse_scalar (amrex::Vector<amrex::MultiFab      *> const& scalar,
@@ -595,6 +602,7 @@ private:
         amrex::MultiFab density;
         amrex::MultiFab density_eb;
         amrex::MultiFab density_o;
+        amrex::MultiFab density_nph;
         amrex::MultiFab tracer;
         amrex::MultiFab tracer_eb;
         amrex::MultiFab tracer_o;
@@ -751,6 +759,7 @@ private:
     amrex::Vector<amrex::MultiFab*> get_tracer_eb () noexcept;
     amrex::Vector<amrex::MultiFab*> get_density_old () noexcept;
     amrex::Vector<amrex::MultiFab*> get_density_new () noexcept;
+    amrex::Vector<amrex::MultiFab*> get_density_nph () noexcept;
     amrex::Vector<amrex::MultiFab*> get_tracer_old () noexcept;
     amrex::Vector<amrex::MultiFab*> get_tracer_new () noexcept;
     amrex::Vector<amrex::MultiFab*> get_mac_phi () noexcept;
@@ -774,6 +783,7 @@ private:
     [[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_tracer_eb () const noexcept;
     [[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_density_old_const () const noexcept;
     [[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_density_new_const () const noexcept;
+    [[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_density_nph_const () const noexcept;
     [[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_tracer_old_const () const noexcept;
     [[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_tracer_new_const () const noexcept;
     [[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_vel_forces_const () const noexcept;