From a0d0d326974c26afdd77c729d5d4b4484d96e6b7 Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Fri, 5 Jan 2024 16:09:21 -0700 Subject: [PATCH] MD's xmas present: EB Godunov (#731) Co-authored-by: Jon Rood Co-authored-by: Bruce Perry Co-authored-by: Bruce Perry <53018946+baperry2@users.noreply.github.com> --- Docs/sphinx/Algorithms.rst | 5 + .../{eb-symmetry.inp => eb-symmetry-mol.inp} | 0 Exec/RegTests/EB-C1/eb-symmetry-plm.inp | 71 ++ Exec/RegTests/EB-C1/ebmms-2.inp | 71 ++ .../{eb-c8-rere.inp => eb-c8-rere-mol.inp} | 0 Exec/RegTests/EB-C8/eb-c8-rere-plm.inp | 90 ++ Exec/RegTests/MassCons/masscons-mol-eb.inp | 84 ++ Exec/RegTests/MassCons/masscons-plm-eb.inp | 84 ++ Exec/RegTests/Sod/shock-cylinder.inp | 101 ++ Source/EB.H | 6 +- Source/Godunov.H | 922 ++++++++++++------ Source/Godunov.cpp | 819 +++++++++++++++- Source/Hydro.H | 784 +++++++++++++++ Source/Hydro.cpp | 412 +++++++- Source/PLM.H | 250 ++++- Source/PeleC.cpp | 4 - Source/main.cpp | 3 +- Tests/CMakeLists.txt | 18 +- 18 files changed, 3365 insertions(+), 359 deletions(-) rename Exec/RegTests/EB-C1/{eb-symmetry.inp => eb-symmetry-mol.inp} (100%) create mode 100644 Exec/RegTests/EB-C1/eb-symmetry-plm.inp create mode 100644 Exec/RegTests/EB-C1/ebmms-2.inp rename Exec/RegTests/EB-C8/{eb-c8-rere.inp => eb-c8-rere-mol.inp} (100%) create mode 100644 Exec/RegTests/EB-C8/eb-c8-rere-plm.inp create mode 100644 Exec/RegTests/MassCons/masscons-mol-eb.inp create mode 100644 Exec/RegTests/MassCons/masscons-plm-eb.inp create mode 100644 Exec/RegTests/Sod/shock-cylinder.inp diff --git a/Docs/sphinx/Algorithms.rst b/Docs/sphinx/Algorithms.rst index e3f3c724f..e4988e69f 100644 --- a/Docs/sphinx/Algorithms.rst +++ b/Docs/sphinx/Algorithms.rst @@ -502,6 +502,11 @@ than PPM simulations. Three dimensional energy spectrum at :math:`t = 5\tau`. Solid red: PPM at :math:`N=128^3`; dashed green: MOL at :math:`N=128^3`; dot-dashed blue: PPM at :math:`N=512^3`; dotted orange: MOL at :math:`N=512^3`; dashed black: spectral code. +Advection schemes and EB +~~~~~~~~~~~~~~~~~~~~~~~~ + +The two main advection schemes (PLM and MOL) described above support the use of complex geometry in the domain (as implemented using the :ref:`Embedded Boundary (EB) formulation `). While the PLM advection scheme is preferred numerically (less dissipative scheme), the PLM implementation with EB is a recent addition and has not had the same extensive testing as the MOL scheme with EB. Both PLM and MOL may also have issues when the EB intersects the boundary at an angle (this remains to be solved). + Diffusion --------- diff --git a/Exec/RegTests/EB-C1/eb-symmetry.inp b/Exec/RegTests/EB-C1/eb-symmetry-mol.inp similarity index 100% rename from Exec/RegTests/EB-C1/eb-symmetry.inp rename to Exec/RegTests/EB-C1/eb-symmetry-mol.inp diff --git a/Exec/RegTests/EB-C1/eb-symmetry-plm.inp b/Exec/RegTests/EB-C1/eb-symmetry-plm.inp new file mode 100644 index 000000000..27aa9bb52 --- /dev/null +++ b/Exec/RegTests/EB-C1/eb-symmetry-plm.inp @@ -0,0 +1,71 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 100 +stop_time = 0.000005 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 0 +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = -1.0 -1.0 -1.0 +geometry.prob_hi = 1.0 1.0 1.0 +# use with single level +amr.n_cell = 16 16 16 + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +pelec.lo_bc = "Hard" "Hard" "Hard" +pelec.hi_bc = "Hard" "Hard" "Hard" + +# WHICH PHYSICS +pelec.do_hydro = 1 +pelec.diffuse_vel = 1 +pelec.diffuse_temp = 1 +pelec.do_mol = 0 +pelec.do_react = 0 +pelec.do_mms = 1 +pelec.masa_solution_name = ad_cns_3d_les_sph + +# TIME STEP CONTROL +pelec.cfl = 0.1 # cfl number for hyperbolic system +pelec.init_shrink = 0.3 # scale back initial timestep +pelec.change_max = 1.1 # max time step growth +pelec.dt_cutoff = 5.e-20 # level 0 timestep below which we halt + +# DIAGNOSTICS & VERBOSITY +pelec.sum_interval = 1 # timesteps between computing mass +pelec.v = 1 # verbosity in PeleC.cpp +amr.v = 1 # verbosity in Amr.cpp +amr.data_log = datlog mmslog +#amr.grid_log = grdlog # name of grid logging file + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed +#amr.max_level = 1 # maximum level number allowed +#amr.ref_ratio = 2 2 2 2 # refinement ratio +#amr.regrid_int = 2 2 2 2 # how often to regrid +amr.blocking_factor = 8 # block factor in grid generation +amr.max_grid_size = 32 +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.check_file = chk # root name of checkpoint file +amr.check_int = 1000 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_file = plt # root name of plotfile +amr.plot_int = 1000 # number of timesteps between plotfiles +amr.plot_vars = density Temp +amr.derive_plot_vars = x_velocity y_velocity z_velocity magvel magvort pressure rhommserror ummserror vmmserror wmmserror pmmserror + +# PROBLEM PARAMETERS +prob.rho_x_fact = 0.0 + +# EB +eb2.geom_type = sphere +eb2.sphere_radius = 0.25 +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_has_fluid_inside = 0 +pelec.eb_isothermal = 0 +ebd.boundary_grad_stencil_type=0 diff --git a/Exec/RegTests/EB-C1/ebmms-2.inp b/Exec/RegTests/EB-C1/ebmms-2.inp new file mode 100644 index 000000000..5eba1cbe8 --- /dev/null +++ b/Exec/RegTests/EB-C1/ebmms-2.inp @@ -0,0 +1,71 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 100 +stop_time = 0.000005 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 0 +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = -1.0 -1.0 -1.0 +geometry.prob_hi = 1.0 1.0 1.0 +# use with single level +amr.n_cell = 16 16 16 + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +pelec.lo_bc = "Hard" "Hard" "Hard" +pelec.hi_bc = "Hard" "Hard" "Hard" + +# WHICH PHYSICS +pelec.do_hydro = 1 +pelec.diffuse_vel = 1 +pelec.diffuse_temp = 1 +pelec.do_mol = 0 +pelec.do_react = 0 +pelec.do_mms = 1 +pelec.masa_solution_name = ad_cns_3d_les_sph + +# TIME STEP CONTROL +pelec.cfl = 0.1 # cfl number for hyperbolic system +pelec.init_shrink = 0.3 # scale back initial timestep +pelec.change_max = 1.1 # max time step growth +pelec.dt_cutoff = 5.e-20 # level 0 timestep below which we halt + +# DIAGNOSTICS & VERBOSITY +pelec.sum_interval = 1 # timesteps between computing mass +pelec.v = 1 # verbosity in PeleC.cpp +amr.v = 1 # verbosity in Amr.cpp +amr.data_log = datlog mmslog +#amr.grid_log = grdlog # name of grid logging file + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed +#amr.max_level = 1 # maximum level number allowed +#amr.ref_ratio = 2 2 2 2 # refinement ratio +#amr.regrid_int = 2 2 2 2 # how often to regrid +amr.blocking_factor = 8 # block factor in grid generation +amr.max_grid_size = 64 +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.check_file = chk # root name of checkpoint file +amr.check_int = 1000 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_file = plt # root name of plotfile +amr.plot_int = 1000 # number of timesteps between plotfiles +amr.plot_vars = density Temp +amr.derive_plot_vars = x_velocity y_velocity z_velocity magvel magvort pressure rhommserror ummserror vmmserror wmmserror pmmserror + +# PROBLEM PARAMETERS + +# EB +eb2.geom_type = sphere +eb2.sphere_radius = 0.25 +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_has_fluid_inside = 0 +pelec.eb_boundary_T = 300. +pelec.eb_isothermal = 1 +ebd.boundary_grad_stencil_type=0 diff --git a/Exec/RegTests/EB-C8/eb-c8-rere.inp b/Exec/RegTests/EB-C8/eb-c8-rere-mol.inp similarity index 100% rename from Exec/RegTests/EB-C8/eb-c8-rere.inp rename to Exec/RegTests/EB-C8/eb-c8-rere-mol.inp diff --git a/Exec/RegTests/EB-C8/eb-c8-rere-plm.inp b/Exec/RegTests/EB-C8/eb-c8-rere-plm.inp new file mode 100644 index 000000000..0f0e00570 --- /dev/null +++ b/Exec/RegTests/EB-C8/eb-c8-rere-plm.inp @@ -0,0 +1,90 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 10000 +stop_time = 1.959e-6 #final time is 0.2*L*sqrt(rhoL/pL) + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 1 +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = 0.0 -0.13 0.0 +geometry.prob_hi = 0.87 0.653 0.0435 +amr.n_cell = 160 144 8 + +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< +# Interior, UserBC, Symmetry, SlipWall, NoSlipWall +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< +pelec.lo_bc = "FOExtrap" "NoSlipWall" "Interior" +pelec.hi_bc = "FOExtrap" "NoSlipWall" "Interior" + +# WHICH PHYSICS +pelec.do_hydro = 1 +pelec.do_mol = 0 +pelec.diffuse_vel = 0 +pelec.diffuse_temp = 0 +pelec.diffuse_spec = 0 +pelec.do_react = 0 +pelec.diffuse_enth = 0 + +# TIME STEP CONTROL +pelec.dt_cutoff = 5.e-20 # level 0 timestep below which we halt +pelec.cfl = 0.8 # cfl number for hyperbolic system +pelec.init_shrink = 0.8 # scale back initial timestep +pelec.change_max = 1.05 # scale back initial timestep + +# DIAGNOSTICS & VERBOSITY +pelec.sum_interval = 1 # timesteps between computing mass +pelec.v = 1 # verbosity in PeleC.cpp +amr.v = 1 # verbosity in Amr.cpp +amr.data_log = datlog + +# REFINEMENT / REGRIDDING +amr.max_level = 1 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 2 2 2 2 # how often to regrid +amr.n_error_buf = 0 0 0 0 # number of buffer cells in error est +amr.blocking_factor = 8 # block factor in grid generation +amr.max_grid_size = 64 +amr.grid_eff = 0.99 + +# CHECKPOINT FILES +amr.checkpoint_files_output = 0 +amr.check_file = chk # root name of checkpoint file +amr.check_int = -1 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_files_output = 1 +amr.plot_file = plt # root name of plotfile +amr.plot_int = 100 # number of timesteps between plotfiles +amr.derive_plot_vars = ALL + +# PROBLEM PARAMETERS +prob.p_l = 1e7 +prob.rho_l = 9.6e-4 +prob.p_r = 1e6 +prob.rho_r = 1.2e-4 +prob.angle = 30.0 +prob.left_gas = N2 +prob.right_gas = HE + +# Problem setup +pelec.eb_boundary_T = 300. +pelec.eb_isothermal = 0 + +# TAGGING +tagging.denerr = 1e20 +tagging.dengrad = 4e-5 +tagging.max_denerr_lev = 3 +tagging.max_dengrad_lev = 3 +tagging.tempgrad = 50.0 +tagging.max_tempgrad_lev = 3 +tagging.max_vfracerr_lev = 0 +tagging.eb_refine_type = "adaptive" +tagging.max_eb_refine_lev = 0 +tagging.min_eb_refine_lev = 0 + +eb2.geom_type = "rotated_box" +eb2.box_lo = -10.0 -0.172 -10.0 +eb2.box_hi = 10.0 0.172 10.0 +eb2.box_rotation = 30.0 +eb2.box_rotation_axe = 2 +eb2.box_has_fluid_inside = 1 +ebd.boundary_grad_stencil_type = 0 \ No newline at end of file diff --git a/Exec/RegTests/MassCons/masscons-mol-eb.inp b/Exec/RegTests/MassCons/masscons-mol-eb.inp new file mode 100644 index 000000000..4deb24c41 --- /dev/null +++ b/Exec/RegTests/MassCons/masscons-mol-eb.inp @@ -0,0 +1,84 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 30 +stop_time = 0.2 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 0 +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = -0.5 -0.5 -0.5 +geometry.prob_hi = 0.5 0.5 0.5 +amr.n_cell = 8 8 8 + +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< +# Interior, UserBC, Symmetry, SlipWall, NoSlipWall +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< + +pelec.lo_bc = "SlipWall" "NoSlipWall" "Symmetry" +pelec.hi_bc = "Hard" "Hard" "Hard" +prob.wall_type = 1 0 1 + +# WHICH PHYSICS +pelec.mol_iorder = 2 +pelec.do_hydro = 1 +pelec.do_mol = 1 +pelec.diffuse_vel = 1 +pelec.diffuse_temp = 1 +pelec.diffuse_spec = 1 +pelec.do_react = 0 +pelec.diffuse_enth = 1 +pelec.add_ext_src = 0 +pelec.external_forcing = 0.0 0.0 0.0 + +transport.const_viscosity = 1 +transport.const_conductivity = 2.7271624e+04 + +# TIME STEP CONTROL +pelec.cfl = 0.1 # cfl number for hyperbolic system +pelec.init_shrink = 1.0 # scale back initial timestep +pelec.change_max = 1.05 # scale back initial timestep +pelec.dt_cutoff = 5.e-20 # level 0 timestep below which we halt + +# DIAGNOSTICS & VERBOSITY +pelec.sum_interval = 1 # timesteps between computing mass +pelec.v = 1 # verbosity in PeleC.cpp +amr.v = 1 # verbosity in Amr.cpp +amr.data_log = datlog + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed +#amr.ref_ratio = 2 2 2 2 # refinement ratio +#amr.regrid_int = 2 2 2 2 # how often to regrid +amr.blocking_factor = 8 # block factor in grid generation +amr.max_grid_size = 64 +amr.n_error_buf = 12 8 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.checkpoint_files_output = 0 +amr.check_file = chk # root name of checkpoint file +amr.check_int = 500 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_files_output = 0 +amr.plot_file = plt # root name of plotfile +amr.plot_int = -1 # number of timesteps between plotfiles +amr.plot_vars = density Temp +amr.derive_plot_vars = x_velocity y_velocity z_velocity magvel magvort pressure +pelec.plot_rhoy = 0 +pelec.plot_massfrac = 1 + +# PROBLEM PARAMETERS +prob.T_mean = 750.0 +prob.u0 = 10000.0 +prob.v0 = 8000.0 +prob.w0 = 5000.0 + +# Problem setup +eb2.geom_type = sphere +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.2 +eb2.sphere_has_fluid_inside = 0 +pelec.eb_isothermal = 0 + +#amrex.fpe_trap_invalid = 1 +#amrex.fpe_trap_zero = 1 +#amrex.fpe_trap_overflow = 1 diff --git a/Exec/RegTests/MassCons/masscons-plm-eb.inp b/Exec/RegTests/MassCons/masscons-plm-eb.inp new file mode 100644 index 000000000..4c2d78f3b --- /dev/null +++ b/Exec/RegTests/MassCons/masscons-plm-eb.inp @@ -0,0 +1,84 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 30 +stop_time = 0.2 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 0 +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = -0.5 -0.5 -0.5 +geometry.prob_hi = 0.5 0.5 0.5 +amr.n_cell = 16 16 16 + +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< +# Interior, UserBC, Symmetry, SlipWall, NoSlipWall +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< + +pelec.lo_bc = "SlipWall" "NoSlipWall" "Symmetry" +pelec.hi_bc = "UserBC" "UserBC" "UserBC" +prob.wall_type = 1 0 1 + +# WHICH PHYSICS +pelec.ppm_type = 0 +pelec.do_hydro = 1 +pelec.do_mol = 0 +pelec.diffuse_vel = 1 +pelec.diffuse_temp = 1 +pelec.diffuse_spec = 1 +pelec.do_react = 0 +pelec.diffuse_enth = 1 +pelec.add_ext_src = 0 +pelec.external_forcing = 0.0 0.0 0.0 + +transport.const_viscosity = 1 +transport.const_conductivity = 2.7271624e+04 + +# TIME STEP CONTROL +pelec.cfl = 0.9 # cfl number for hyperbolic system +pelec.init_shrink = 1.0 # scale back initial timestep +pelec.change_max = 1.05 # scale back initial timestep +pelec.dt_cutoff = 5.e-20 # level 0 timestep below which we halt + +# DIAGNOSTICS & VERBOSITY +pelec.sum_interval = 1 # timesteps between computing mass +pelec.v = 1 # verbosity in PeleC.cpp +amr.v = 1 # verbosity in Amr.cpp +amr.data_log = datlog + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed +#amr.ref_ratio = 2 2 2 2 # refinement ratio +#amr.regrid_int = 2 2 2 2 # how often to regrid +amr.blocking_factor = 8 # block factor in grid generation +amr.max_grid_size = 64 +amr.n_error_buf = 12 8 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.checkpoint_files_output = 0 +amr.check_file = chk # root name of checkpoint file +amr.check_int = 500 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_files_output = 0 +amr.plot_file = plt # root name of plotfile +amr.plot_int = -1 # number of timesteps between plotfiles +amr.plot_vars = density Temp +amr.derive_plot_vars = x_velocity y_velocity z_velocity magvel magvort pressure +pelec.plot_rhoy = 0 +pelec.plot_massfrac = 1 + +# PROBLEM PARAMETERS +prob.T_mean = 750.0 +prob.u0 = 10000.0 +prob.v0 = 8000.0 +prob.w0 = 5000.0 + +# Problem setup +eb2.geom_type = sphere +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.2 +eb2.sphere_has_fluid_inside = 0 +pelec.eb_isothermal = 0 + +#amrex.fpe_trap_invalid = 1 +#amrex.fpe_trap_zero = 1 +#amrex.fpe_trap_overflow = 1 diff --git a/Exec/RegTests/Sod/shock-cylinder.inp b/Exec/RegTests/Sod/shock-cylinder.inp new file mode 100644 index 000000000..7eea690cb --- /dev/null +++ b/Exec/RegTests/Sod/shock-cylinder.inp @@ -0,0 +1,101 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 10000 +stop_time = 0.2 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 0 +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = 0 0 0 +geometry.prob_hi = 1.0 1.0 1.0 +#amr.n_cell = 32 8 8 +amr.n_cell = 16 16 16 + +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< +# Interior, UserBC, Symmetry, SlipWall, NoSlipWall +# >>>>>>>>>>>>> BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<< +pelec.lo_bc = "Hard" "SlipWall" "SlipWall" +pelec.hi_bc = "Hard" "SlipWall" "SlipWall" + +# WHICH PHYSICS +pelec.do_hydro = 1 +pelec.diffuse_vel = 0 +pelec.diffuse_temp = 0 +pelec.diffuse_spec = 0 +pelec.do_react = 0 +pelec.do_mol = 0 + +# TIME STEP CONTROL +pelec.cfl = 0.9 # cfl number for hyperbolic system +pelec.init_shrink = 0.1 # scale back initial timestep +pelec.change_max = 1.05 # scale back initial timestep +pelec.dt_cutoff = 5.e-20 # level 0 timestep below which we halt + +# DIAGNOSTICS & VERBOSITY +pelec.sum_interval = 1 # timesteps between computing mass +pelec.v = 1 # verbosity in PeleC.cpp +amr.v = 1 # verbosity in Amr.cpp +amr.data_log = datlog + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed +# amr.ref_ratio = 2 2 2 2 # refinement ratio +# amr.regrid_int = 2 2 2 2 # how often to regrid +amr.blocking_factor = 8 # block factor in grid generation +amr.max_grid_size = 8 +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.checkpoint_files_output = 0 +amr.check_file = chk # root name of checkpoint file +amr.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_files_output = 1 +amr.plot_file = plt # root name of plotfile +amr.plot_int = 100 # number of timesteps between plotfiles +amr.plot_vars = density Temp +amr.derive_plot_vars = ALL + +# PROBLEM PARAMETERS +prob.rho_l = 3.6737153092795505 +prob.u_l = 3.8260444105770732 +prob.p_l = 22.613625000000006 +prob.rho_r = 1.0 +prob.u_r = 0.0 +prob.p_r = 2.5 +prob.idir = 1 +prob.frac = 0.2 + +# TAGGING +tagging.denerr = 3 +tagging.dengrad = 0.01 +tagging.max_denerr_lev = 3 +tagging.max_dengrad_lev = 3 +tagging.presserr = 3 +tagging.pressgrad = 0.01 +tagging.max_presserr_lev = 3 +tagging.max_pressgrad_lev = 3 + +# EB +ebd.boundary_grad_stencil_type = 0 +pelec.eb_boundary_T = 300. +pelec.eb_isothermal = 1 + +# for 3D: +eb2.geom_type = "cylinder" +eb2.cylinder_direction = 2 +eb2.cylinder_center = 0.5 0.5 0.5 +eb2.cylinder_radius = 0.1 +eb2.cylinder_height = 1000.0 +eb2.cylinder_has_fluid_inside = 0 + +# for 2D (need a sphere): +# eb2.geom_type = "sphere" +# eb2.sphere_radius = 0.1 +# eb2.sphere_center = 0.5 0.5 0.5 +# eb2.sphere_has_fluid_inside = 0 + +# fabarray.mfiter_tile_size = 1024000 64 64 +# amrex.fpe_trap_invalid = 1 +# amrex.fpe_trap_zero = 1 +# amrex.fpe_trap_overflow = 1 diff --git a/Source/EB.H b/Source/EB.H index c7a9768bc..4c3d3ee45 100644 --- a/Source/EB.H +++ b/Source/EB.H @@ -186,15 +186,15 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE bool is_cut_neighborhood( - const amrex::IntVect iv, + const amrex::IntVect& iv, amrex::Array4 const& flags, - const int ng = 1) + const amrex::IntVect& ng = amrex::IntVect(1)) { // Check if there's a cut cell in an ng neighborhood around iv bool has_cut_cell = flags(iv).isSingleValued(); for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { const auto dimvec = amrex::IntVect::TheDimensionVector(idir); - for (int n = 1; n <= ng; n++) { + for (int n = 1; n <= ng[idir]; n++) { has_cut_cell |= flags(iv - n * dimvec).isSingleValued() || flags(iv + n * dimvec).isSingleValued(); } diff --git a/Source/Godunov.H b/Source/Godunov.H index 523b4e457..3e0592709 100644 --- a/Source/Godunov.H +++ b/Source/Godunov.H @@ -7,6 +7,7 @@ #include "Constants.H" #include "IndexDefines.H" #include "Riemann.H" +#include "EB.H" AMREX_GPU_DEVICE AMREX_FORCE_INLINE @@ -28,31 +29,47 @@ AMREX_FORCE_INLINE void pc_transdo_passive( const amrex::IntVect iv, - const amrex::IntVect ivp0, - const amrex::IntVect ivp1, + const amrex::IntVect ivpn, + const amrex::IntVect ivpt, const int n, const int u_offset, const int q_offset, const amrex::Real cdtdx, const amrex::Real flxrho, - const amrex::Real rrr, - const amrex::Real rrl, - amrex::Array4 const& qyp, - amrex::Array4 const& qym, + const bool divide_by_rho, + amrex::Array4 const& qnormp, + amrex::Array4 const& qnormm, amrex::Array4 const& flxx, amrex::Array4 const& qp, - amrex::Array4 const& qm) + amrex::Array4 const& qm, + const bool no_cov_face, + const bool lo_face_not_covered, + const bool hi_face_not_covered) { const int uc = u_offset + n; const int qc = q_offset + n; - const amrex::Real compn = cdtdx * (flxx(ivp0, uc) - flxx(iv, uc)); - amrex::Real rrnew = rrr - flxrho; - amrex::Real compo = rrr * qyp(iv, qc) - compn; - qp(iv, qc) = compo / rrnew; - - rrnew = rrl - flxrho; - compo = rrl * qym(ivp1, qc) - compn; - qm(ivp1, qc) = compo / rrnew; + const amrex::Real compn = + no_cov_face ? cdtdx * (flxx(ivpt, uc) - flxx(iv, uc)) : amrex::Real(0.0); + + amrex::Real rr, rrnew, compo; + + if (lo_face_not_covered) { + rr = qnormp(iv, QRHO); + rrnew = divide_by_rho ? rr - flxrho : 1.0; + compo = rr * qnormp(iv, qc) - compn; + qp(iv, qc) = compo / rrnew; + } else { + qp(iv, qc) = qnormp(iv, qc); + } + + if (hi_face_not_covered) { + rr = qnormm(ivpn, QRHO); + rrnew = divide_by_rho ? rr - flxrho : 1.0; + compo = rr * qnormm(ivpn, qc) - compn; + qm(ivpn, qc) = compo / rrnew; + } else { + qm(ivpn, qc) = qnormm(ivpn, qc); + } } AMREX_GPU_DEVICE @@ -60,39 +77,58 @@ AMREX_FORCE_INLINE void pc_transdd_passive( const amrex::IntVect iv, - const amrex::IntVect ivp0, - const amrex::IntVect ivp1, - const amrex::IntVect ivp2, + const amrex::IntVect ivpn, + const amrex::IntVect ivpt0, + const amrex::IntVect ivpt1, const int n, const int u_offset, const int q_offset, - const amrex::Real cdtdx, - const amrex::Real cdtdy, + const amrex::Real cdtdx0, + const amrex::Real cdtdx1, const amrex::Real hdt, const amrex::Real rrr, const amrex::Real rrl, const amrex::Real rrnewr, const amrex::Real rrnewl, amrex::Array4 const& srcq, - amrex::Array4 const& qzp, - amrex::Array4 const& qzm, + amrex::Array4 const& qnormp, + amrex::Array4 const& qnormm, amrex::Array4 const& flxx, amrex::Array4 const& flxy, amrex::Array4 const& qp, - amrex::Array4 const& qm) + amrex::Array4 const& qm, + const bool no_cov_face, + const bool lo_face_not_covered, + const bool hi_face_not_covered) { const int uc = u_offset + n; const int qc = q_offset + n; const amrex::Real srcpass = srcq(iv, qc); - const amrex::Real compn = cdtdx * (flxx(ivp0, uc) - flxx(iv, uc)) + - cdtdy * (flxy(ivp1, uc) - flxy(iv, uc)); + + const amrex::Real compn = no_cov_face + ? cdtdx0 * (flxx(ivpt0, uc) - flxx(iv, uc)) + + cdtdx1 * (flxy(ivpt1, uc) - flxy(iv, uc)) + : amrex::Real(0.0); + + amrex::Real compo; + // qp - amrex::Real compo = rrr * qzp(iv, qc) - compn; - qp(iv, qc) = compo / rrnewr + hdt * srcpass; + if (lo_face_not_covered) { + compo = rrr * qnormp(iv, qc) - compn; + qp(iv, qc) = compo / rrnewr + hdt * srcpass; + } else { + qp(iv, qc) = qnormp(iv, qc); + } + // qm - compo = rrl * qzm(ivp2, qc) - compn; - qm(ivp2, qc) = compo / rrnewl + hdt * srcpass; + + if (hi_face_not_covered) { + compo = rrl * qnormm(ivpn, qc) - compn; + qm(ivpn, qc) = compo / rrnewl + hdt * srcpass; + } else { + qm(ivpn, qc) = qnormm(ivpn, qc); + } } AMREX_GPU_DEVICE @@ -100,34 +136,50 @@ AMREX_FORCE_INLINE void pc_transd_passive( const amrex::IntVect iv, - const amrex::IntVect ivpt, const amrex::IntVect ivpn, + const amrex::IntVect ivpt, const int n, const int u_offset, const int q_offset, const amrex::Real cdtdx, const amrex::Real hdt, const amrex::Real flxrho, - const amrex::Real rrr, - const amrex::Real rrl, + const bool divide_by_rho, amrex::Array4 const& srcQ, - amrex::Array4 const& qyp, - amrex::Array4 const& qym, + amrex::Array4 const& qnormp, + amrex::Array4 const& qnormm, amrex::Array4 const& flxx, amrex::Array4 const& qp, - amrex::Array4 const& qm) + amrex::Array4 const& qm, + const bool no_cov_face, + const bool lo_face_not_covered, + const bool hi_face_not_covered) { const int uc = u_offset + n; const int qc = q_offset + n; const amrex::Real srcpass = srcQ(iv, qc); - const amrex::Real compn = cdtdx * (flxx(ivpt, uc) - flxx(iv, uc)); - amrex::Real rrnew = rrr - flxrho; - amrex::Real compo = rrr * qyp(iv, qc) - compn; - qp(iv, qc) = compo / rrnew + hdt * srcpass; - - rrnew = rrl - flxrho; - compo = rrl * qym(ivpn, qc) - compn; - qm(ivpn, qc) = compo / rrnew + hdt * srcpass; + const amrex::Real compn = + no_cov_face ? cdtdx * (flxx(ivpt, uc) - flxx(iv, uc)) : 0.0; + + amrex::Real rr, rrnew, compo; + + if (lo_face_not_covered) { + rr = qnormp(iv, QRHO); + rrnew = divide_by_rho ? rr - flxrho : 1.0; + compo = rr * qnormp(iv, qc) - compn; + qp(iv, qc) = compo / rrnew + hdt * srcpass; + } else { + qp(iv, qc) = qnormp(iv, qc); + } + + if (hi_face_not_covered) { + rr = qnormm(ivpn, QRHO); + rrnew = divide_by_rho ? rr - flxrho : 1.0; + compo = rr * qnormm(ivpn, qc) - compn; + qm(ivpn, qc) = compo / rrnew + hdt * srcpass; + } else { + qm(ivpn, qc) = qnormm(ivpn, qc); + } } AMREX_GPU_DEVICE @@ -207,6 +259,44 @@ flatten( return 1.0 - amrex::max(chi2 * z2, chi * z); } +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +flatten_eb( + const int i, + const int j, + const int k, + const int dir, + amrex::Array4 const& flags, + amrex::Array4 const& q) +{ + + const amrex::IntVect iv{AMREX_D_DECL(i, j, k)}; + const amrex::IntVect dvec(amrex::IntVect::TheDimensionVector(dir)); + + const amrex::IntVect ivm3(iv - 3 * dvec); + const amrex::IntVect ivp3(iv + 3 * dvec); + +#if (AMREX_SPACEDIM == 2) + if ( + (!q.contains(ivp3[0], ivp3[1], k)) || (!q.contains(ivm3[0], ivm3[1], k))) { + return 1.0; + } +#elif (AMREX_SPACEDIM == 3) + if ( + !q.contains(ivp3[0], ivp3[1], ivp3[2]) || + !q.contains(ivm3[0], ivm3[1], ivm3[2])) { + return 1.0; + } +#endif + + const amrex::IntVect ng(3 * dvec); + if (is_cut_neighborhood(iv, flags, ng)) { + return 1.0; + } + return flatten(AMREX_D_DECL(i, j, k), dir, q); +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -373,125 +463,185 @@ pc_transdo( const int other_dir, amrex::Array4 const& qm, amrex::Array4 const& qp, - amrex::Array4 const& qym, - amrex::Array4 const& qyp, + amrex::Array4 const& qnormm, + amrex::Array4 const& qnormp, amrex::Array4 const& flxx, amrex::Array4 const& qa, amrex::Array4 const& qint, - const amrex::Real cdtdx) + const amrex::Real cdtdx, + amrex::Array4 const& norm_area = {}, + amrex::Array4 const& trans_area = {}) { const amrex::IntVect iv{AMREX_D_DECL(i, j, k)}; - const amrex::IntVect ivp0(iv + amrex::IntVect::TheDimensionVector(other_dir)); - const amrex::IntVect ivp1(iv + amrex::IntVect::TheDimensionVector(dir)); + const amrex::IntVect ivpn(iv + amrex::IntVect::TheDimensionVector(dir)); + const amrex::IntVect ivpt(iv + amrex::IntVect::TheDimensionVector(other_dir)); const int qvidx = (other_dir == 0) ? GDU : ((other_dir == 1) ? GDV : GDW); - const amrex::Real flxrho = cdtdx * (flxx(ivp0, URHO) - flxx(iv, URHO)); - const amrex::Real flxu = cdtdx * (flxx(ivp0, UMX) - flxx(iv, UMX)); - const amrex::Real flxv = cdtdx * (flxx(ivp0, UMY) - flxx(iv, UMY)); - const amrex::Real flxw = cdtdx * (flxx(ivp0, UMZ) - flxx(iv, UMZ)); - const amrex::Real flxe = cdtdx * (flxx(ivp0, UEDEN) - flxx(iv, UEDEN)); - const amrex::Real rrr = qyp(iv, QRHO); - const amrex::Real rrl = qym(ivp1, QRHO); + // clang-format off + // ivpn is an offset by 1 in the normal direction + // ivpt is an offset by 1 in the transverse direction + // If (dir == 0 and other_dir == 1), ivpt = (i,j+1,k) and ivpn = (i+1,j,k) and qvidx = GDV + // If (dir == 0 and other_dir == 2), ivpt = (i,j,k+1) and ivpn = (i+1,j,k) and qvidx = GDW + // If (dir == 1 and other_dir == 0), ivpt = (i+1,j,k) and ivpn = (i,j+1,k) and qvidx = GDU + // If (dir == 1 and other_dir == 2), ivpt = (i,j,k+1) and ivpn = (i,j+1,k) and qvidx = GDW + // If (dir == 2 and other_dir == 0), ivpt = (i+1,j,k) and ivpn = (i,j,k+1) and qvidx = GDU + // If (dir == 2 and other_dir == 1), ivpt = (i,j+1,k) and ivpn = (i,j,k+1) and qvidx = GDV + // clang-format on + + bool no_cov_face = true; + bool lo_face_not_covered = true; + bool hi_face_not_covered = true; + + if (trans_area) { + no_cov_face = (trans_area(ivpt) > 0.0) && (trans_area(iv) > 0.0); + } + if (norm_area) { + lo_face_not_covered = (norm_area(iv) > amrex::Real(0.0)); + hi_face_not_covered = (norm_area(ivpn) > amrex::Real(0.0)); + } + + const amrex::Real flxrho = no_cov_face + ? cdtdx * (flxx(ivpt, URHO) - flxx(iv, URHO)) + : amrex::Real(0.0); + const amrex::Real flxu = + no_cov_face ? cdtdx * (flxx(ivpt, UMX) - flxx(iv, UMX)) : amrex::Real(0.0); + const amrex::Real flxv = + no_cov_face ? cdtdx * (flxx(ivpt, UMY) - flxx(iv, UMY)) : amrex::Real(0.0); + const amrex::Real flxw = + no_cov_face ? cdtdx * (flxx(ivpt, UMZ) - flxx(iv, UMZ)) : amrex::Real(0.0); + const amrex::Real flxe = no_cov_face + ? cdtdx * (flxx(ivpt, UEDEN) - flxx(iv, UEDEN)) + : amrex::Real(0.0); const amrex::Real c = qa(iv, QGAMC); // Update passive variables #if NUM_ADV > 0 for (int n = 0; n < NUM_ADV; n++) { pc_transdo_passive( - iv, ivp0, ivp1, n, UFA, QFA, cdtdx, flxrho, rrr, rrl, qyp, qym, flxx, qp, - qm); + iv, ivpn, ivpt, n, UFA, QFA, cdtdx, flxrho, true, qnormp, qnormm, flxx, + qp, qm, no_cov_face, lo_face_not_covered, hi_face_not_covered); } #endif for (int n = 0; n < NUM_SPECIES; n++) { pc_transdo_passive( - iv, ivp0, ivp1, n, UFS, QFS, cdtdx, flxrho, rrr, rrl, qyp, qym, flxx, qp, - qm); + iv, ivpn, ivpt, n, UFS, QFS, cdtdx, flxrho, true, qnormp, qnormm, flxx, + qp, qm, no_cov_face, lo_face_not_covered, hi_face_not_covered); } #if NUM_AUX > 0 for (int n = 0; n < NUM_AUX; n++) { pc_transdo_passive( - iv, ivp0, ivp1, n, UFX, QFX, cdtdx, flxrho, rrr, rrl, qyp, qym, flxx, qp, - qm); + iv, ivpn, ivpt, n, UFX, QFX, cdtdx, flxrho, true, qnormp, qnormm, flxx, + qp, qm, no_cov_face, lo_face_not_covered, hi_face_not_covered); } #endif #if NUM_LIN > 0 for (int n = 0; n < NUM_LIN; n++) { pc_transdo_passive( - iv, ivp0, ivp1, n, ULIN, QLIN, cdtdx, 0., 1., 1., qyp, qym, flxx, qp, qm); + iv, ivpn, ivpt, n, ULIN, QLIN, cdtdx, 0., false, qnormp, qnormm, flxx, qp, + qm, no_cov_face, lo_face_not_covered, hi_face_not_covered); } #endif // Update hydro vars - const amrex::Real pggp = qint(ivp0, GDPRES); - const amrex::Real pggm = qint(iv, GDPRES); - const amrex::Real ugp = qint(ivp0, qvidx); - const amrex::Real ugm = qint(iv, qvidx); - - const amrex::Real dAup = pggp * ugp - pggm * ugm; - const amrex::Real pav = 0.5 * (pggp + pggm); - const amrex::Real dAu = ugp - ugm; + amrex::Real dAup = 0.0, pav = 0.0, dAu = 0.0; + if (no_cov_face) { + const amrex::Real pggp = qint(ivpt, GDPRES); + const amrex::Real pggm = qint(iv, GDPRES); + const amrex::Real ugp = qint(ivpt, qvidx); + const amrex::Real ugm = qint(iv, qvidx); + + dAup = pggp * ugp - pggm * ugm; + pav = 0.5 * (pggp + pggm); + dAu = ugp - ugm; + } // QP - // Convert to conservative - amrex::Real rur = qyp(iv, QU); - amrex::Real rvr = qyp(iv, QV); - amrex::Real rwr = qyp(iv, QW); - const amrex::Real ekinr = 0.5 * rrr * (rur * rur + rvr * rvr + rwr * rwr); - rur *= rrr; - rvr *= rrr; - rwr *= rrr; - - const amrex::Real rer = qyp(iv, QREINT) + ekinr; - // Add transverse predictor - const amrex::Real rrnewr = rrr - flxrho; - const amrex::Real runewr = rur - flxu; - const amrex::Real rvnewr = rvr - flxv; - const amrex::Real rwnewr = rwr - flxw; - const amrex::Real renewr = rer - flxe; - - // Convert back to primitive - qp(iv, QRHO) = rrnewr; - qp(iv, QU) = runewr / rrnewr; - qp(iv, QV) = rvnewr / rrnewr; - qp(iv, QW) = rwnewr / rrnewr; - const amrex::Real rhoekinr = - 0.5 * (runewr * runewr + rvnewr * rvnewr + rwnewr * rwnewr) / rrnewr; - amrex::Real pnewr = qyp(iv, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.0)); - qp(iv, QPRES) = - amrex::max(pnewr, std::numeric_limits::min()); - qp(iv, QREINT) = renewr - rhoekinr; + if (lo_face_not_covered) { + // Convert to conservative + const amrex::Real rrr = qnormp(iv, QRHO); + amrex::Real rur = qnormp(iv, QU); + amrex::Real rvr = qnormp(iv, QV); + amrex::Real rwr = qnormp(iv, QW); + const amrex::Real ekinr = 0.5 * rrr * (rur * rur + rvr * rvr + rwr * rwr); + rur *= rrr; + rvr *= rrr; + rwr *= rrr; + + const amrex::Real rer = qnormp(iv, QREINT) + ekinr; + // Add transverse predictor + amrex::Real rrnewr = rrr - flxrho; + amrex::Real runewr = rur - flxu; + amrex::Real rvnewr = rvr - flxv; + amrex::Real rwnewr = rwr - flxw; + amrex::Real renewr = rer - flxe; + + // Convert back to primitive + qp(iv, QRHO) = rrnewr; + qp(iv, QU) = runewr / rrnewr; + qp(iv, QV) = rvnewr / rrnewr; + qp(iv, QW) = rwnewr / rrnewr; + + const amrex::Real rhoekinr = + 0.5 * (runewr * runewr + rvnewr * rvnewr + rwnewr * rwnewr) / rrnewr; + + qp(iv, QREINT) = renewr - rhoekinr; + + amrex::Real pnewr = + qnormp(iv, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.0)); + ; + qp(iv, QPRES) = amrex::max(pnewr, constants::very_small_num()); + } else { + qp(iv, QRHO) = qnormp(iv, QRHO); + qp(iv, QU) = qnormp(iv, QU); + qp(iv, QV) = qnormp(iv, QV); + qp(iv, QREINT) = qnormp(iv, QREINT); + qp(iv, QPRES) = qnormp(iv, QPRES); + } + // QM - // Conversion to Conservative - amrex::Real rul = qym(ivp1, QU); - amrex::Real rvl = qym(ivp1, QV); - amrex::Real rwl = qym(ivp1, QW); - const amrex::Real ekinl = 0.5 * rrl * (rul * rul + rvl * rvl + rwl * rwl); - rul *= rrl; - rvl *= rrl; - rwl *= rrl; - const amrex::Real rel = qym(ivp1, QREINT) + ekinl; - - // Transverse fluxes - const amrex::Real rrnewl = rrl - flxrho; - const amrex::Real runewl = rul - flxu; - const amrex::Real rvnewl = rvl - flxv; - const amrex::Real rwnewl = rwl - flxw; - const amrex::Real renewl = rel - flxe; - - qm(ivp1, QRHO) = rrnewl; - qm(ivp1, QU) = runewl / rrnewl; - qm(ivp1, QV) = rvnewl / rrnewl; - qm(ivp1, QW) = rwnewl / rrnewl; - const amrex::Real rhoekinl = - 0.5 * (runewl * runewl + rvnewl * rvnewl + rwnewl * rwnewl) / rrnewl; - - qm(ivp1, QREINT) = renewl - rhoekinl; - const amrex::Real pnewl = - qym(ivp1, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.0)); - qm(ivp1, QPRES) = - amrex::max(pnewl, std::numeric_limits::min()); + if (hi_face_not_covered) { + // Conversion to Conservative + const amrex::Real rrl = qnormm(ivpn, QRHO); + amrex::Real rul = qnormm(ivpn, QU); + amrex::Real rvl = qnormm(ivpn, QV); + amrex::Real rwl = qnormm(ivpn, QW); + const amrex::Real ekinl = 0.5 * rrl * (rul * rul + rvl * rvl + rwl * rwl); + rul *= rrl; + rvl *= rrl; + rwl *= rrl; + const amrex::Real rel = qnormm(ivpn, QREINT) + ekinl; + + // Transverse fluxes + amrex::Real rrnewl = rrl - flxrho; + amrex::Real runewl = rul - flxu; + amrex::Real rvnewl = rvl - flxv; + amrex::Real rwnewl = rwl - flxw; + amrex::Real renewl = rel - flxe; + + qm(ivpn, QRHO) = rrnewl; + qm(ivpn, QU) = runewl / rrnewl; + qm(ivpn, QV) = rvnewl / rrnewl; + qm(ivpn, QW) = rwnewl / rrnewl; + const amrex::Real rhoekinl = + 0.5 * (runewl * runewl + rvnewl * rvnewl + rwnewl * rwnewl) / rrnewl; + + qm(ivpn, QREINT) = renewl - rhoekinl; + + amrex::Real pnewl = + qnormm(ivpn, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.0)); + qm(ivpn, QPRES) = + amrex::max(pnewl, constants::very_small_num()); + + } else { + qm(ivpn, QRHO) = qnormm(ivpn, QRHO); + qm(ivpn, QU) = qnormm(ivpn, QU); + qm(ivpn, QV) = qnormm(ivpn, QV); + qm(ivpn, QW) = qnormm(ivpn, QW); + qm(ivpn, QREINT) = qnormm(ivpn, QREINT); + qm(ivpn, QPRES) = qnormm(ivpn, QPRES); + } } // dir corrected from other two dirs @@ -503,8 +653,8 @@ pc_transdd( const int dir, amrex::Array4 const& qm, amrex::Array4 const& qp, - amrex::Array4 const& qzm, - amrex::Array4 const& qzp, + amrex::Array4 const& qnormm, + amrex::Array4 const& qnormp, amrex::Array4 const& flxx, amrex::Array4 const& flxy, amrex::Array4 const& qx, @@ -513,7 +663,10 @@ pc_transdd( amrex::Array4 const& srcq, const amrex::Real hdt, const amrex::Real cdtdx0, - const amrex::Real cdtdx1) + const amrex::Real cdtdx1, + amrex::Array4 const& norm_area = {}, + amrex::Array4 const& trans0_area = {}, + amrex::Array4 const& trans1_area = {}) { const amrex::IntVect iv{AMREX_D_DECL(i, j, k)}; const int qvidx0 = (dir == 0) ? GDV : GDU; @@ -525,20 +678,53 @@ pc_transdd( {bdim[0] * 0 + bdim[1] * 1 + bdim[2] * 2, bdim[0] * 1 + bdim[1] * 0 + bdim[2] * 0, bdim[0] * 2 + bdim[1] * 2 + bdim[2] * 1}}; - const amrex::IntVect ivp0(iv + amrex::IntVect::TheDimensionVector(l_idx[1])); - const amrex::IntVect ivp1(iv + amrex::IntVect::TheDimensionVector(l_idx[2])); - const amrex::IntVect ivp2(iv + amrex::IntVect::TheDimensionVector(dir)); - - const amrex::Real flxrho = cdtdx0 * (flxx(ivp0, URHO) - flxx(iv, URHO)) + - cdtdx1 * (flxy(ivp1, URHO) - flxy(iv, URHO)); - const amrex::Real flxu = cdtdx0 * (flxx(ivp0, UMX) - flxx(iv, UMX)) + - cdtdx1 * (flxy(ivp1, UMX) - flxy(iv, UMX)); - const amrex::Real flxv = cdtdx0 * (flxx(ivp0, UMY) - flxx(iv, UMY)) + - cdtdx1 * (flxy(ivp1, UMY) - flxy(iv, UMY)); - const amrex::Real flxw = cdtdx0 * (flxx(ivp0, UMZ) - flxx(iv, UMZ)) + - cdtdx1 * (flxy(ivp1, UMZ) - flxy(iv, UMZ)); - const amrex::Real flxe = cdtdx0 * (flxx(ivp0, UEDEN) - flxx(iv, UEDEN)) + - cdtdx1 * (flxy(ivp1, UEDEN) - flxy(iv, UEDEN)); + const amrex::IntVect ivpt0(iv + amrex::IntVect::TheDimensionVector(l_idx[1])); + const amrex::IntVect ivpt1(iv + amrex::IntVect::TheDimensionVector(l_idx[2])); + const amrex::IntVect ivpn(iv + amrex::IntVect::TheDimensionVector(dir)); + + // clang-format off + // ivpn is an offset by 1 in the normal direction + // ivpt0 is an offset by 1 in the first transverse direction + // ivpt1 is an offset by 1 in the second transverse direction + // If (dir == 0), ivpt0 = (i,j+1,k) and ivpt1 = (i,j,k+1) and ivpn = (i+1,j,k) and qvidx0 = GDV and qvidx1 = GVW + // If (dir == 1), ivpt0 = (i+1,j,k) and ivpt1 = (i,j,k+1) and ivpn = (i,j+1,k) and qvidx0 = GDU and qvidx1 = GVW + // If (dir == 2), ivpt0 = (i+1,j,k) and ivpt1 = (i,j+1,k) and ivpn = (i,j,k+1) and qvidx0 = GDU and qvidx1 = GVV + // clang-format on + + bool no_cov_face = true; + bool lo_face_not_covered = true; + bool hi_face_not_covered = true; +#ifdef AMREX_USE_EB + if (trans0_area && trans1_area) { + no_cov_face = (trans0_area(ivpt0) > 0.0) && (trans1_area(ivpt1) > 0.0) && + (trans0_area(iv) > 0.0) && (trans1_area(iv) > 0.0); + } + if (norm_area) { + lo_face_not_covered = (norm_area(iv) > amrex::Real(0.0)); + hi_face_not_covered = (norm_area(ivpn) > amrex::Real(0.0)); + } +#endif + + const amrex::Real flxrho = no_cov_face + ? cdtdx0 * (flxx(ivpt0, URHO) - flxx(iv, URHO)) + + cdtdx1 * (flxy(ivpt1, URHO) - flxy(iv, URHO)) + : amrex::Real(0.0); + const amrex::Real flxu = no_cov_face + ? cdtdx0 * (flxx(ivpt0, UMX) - flxx(iv, UMX)) + + cdtdx1 * (flxy(ivpt1, UMX) - flxy(iv, UMX)) + : amrex::Real(0.0); + const amrex::Real flxv = no_cov_face + ? cdtdx0 * (flxx(ivpt0, UMY) - flxx(iv, UMY)) + + cdtdx1 * (flxy(ivpt1, UMY) - flxy(iv, UMY)) + : amrex::Real(0.0); + const amrex::Real flxw = no_cov_face + ? cdtdx0 * (flxx(ivpt0, UMZ) - flxx(iv, UMZ)) + + cdtdx1 * (flxy(ivpt1, UMZ) - flxy(iv, UMZ)) + : amrex::Real(0.0); + const amrex::Real flxe = no_cov_face + ? cdtdx0 * (flxx(ivpt0, UEDEN) - flxx(iv, UEDEN)) + + cdtdx1 * (flxy(ivpt1, UEDEN) - flxy(iv, UEDEN)) + : amrex::Real(0.0); const amrex::Real c = qa(iv, QGAMC); const amrex::Real srcrho = srcq(iv, QRHO); const amrex::Real srcu = srcq(iv, QU); @@ -548,103 +734,134 @@ pc_transdd( const amrex::Real srcp = srcq(iv, QPRES); // Update passive variables - const amrex::Real rrr = qzp(iv, QRHO); - const amrex::Real rrl = qzm(ivp2, QRHO); + const amrex::Real rrr = qnormp(iv, QRHO); + const amrex::Real rrl = qnormm(ivpn, QRHO); const amrex::Real rrnewl = rrl - flxrho; const amrex::Real rrnewr = rrr - flxrho; #if NUM_ADV > 0 for (int n = 0; n < NUM_ADV; n++) { pc_transdd_passive( - iv, ivp0, ivp1, ivp2, n, UFA, QFA, cdtdx0, cdtdx1, hdt, rrr, rrl, rrnewr, - rrnewl, srcq, qzp, qzm, flxx, flxy, qp, qm); + iv, ivpn, ivpt0, ivpt1, n, UFA, QFA, cdtdx0, cdtdx1, hdt, rrr, rrl, + rrnewr, rrnewl, srcq, qnormp, qnormm, flxx, flxy, qp, qm, no_cov_face, + lo_face_not_covered, hi_face_not_covered); } #endif for (int n = 0; n < NUM_SPECIES; n++) { pc_transdd_passive( - iv, ivp0, ivp1, ivp2, n, UFS, QFS, cdtdx0, cdtdx1, hdt, rrr, rrl, rrnewr, - rrnewl, srcq, qzp, qzm, flxx, flxy, qp, qm); + iv, ivpn, ivpt0, ivpt1, n, UFS, QFS, cdtdx0, cdtdx1, hdt, rrr, rrl, + rrnewr, rrnewl, srcq, qnormp, qnormm, flxx, flxy, qp, qm, no_cov_face, + lo_face_not_covered, hi_face_not_covered); } #if NUM_AUX > 0 for (int n = 0; n < NUM_AUX; n++) { pc_transdd_passive( - iv, ivp0, ivp1, ivp2, n, UFX, QFX, cdtdx0, cdtdx1, hdt, rrr, rrl, rrnewr, - rrnewl, srcq, qzp, qzm, flxx, flxy, qp, qm); + iv, ivpn, ivpt0, ivpt1, n, UFX, QFX, cdtdx0, cdtdx1, hdt, rrr, rrl, + rrnewr, rrnewl, srcq, qnormp, qnormm, flxx, flxy, qp, qm, no_cov_face, + lo_face_not_covered, hi_face_not_covered); } #endif #if NUM_LIN > 0 for (int n = 0; n < NUM_LIN; n++) { pc_transdd_passive( - iv, ivp0, ivp1, ivp2, n, ULIN, QLIN, cdtdx0, cdtdx1, hdt, 1., 1., 1., 1., - srcq, qzp, qzm, flxx, flxy, qp, qm); + iv, ivpn, ivpt0, ivpt1, n, ULIN, QLIN, cdtdx0, cdtdx1, hdt, 1., 1., 1., + 1., srcq, qnormp, qnormm, flxx, flxy, qp, qm, no_cov_face, + lo_face_not_covered, hi_face_not_covered); } #endif // Update hydro vars - const amrex::Real pggpx = qx(ivp0, GDPRES); - const amrex::Real pggmx = qx(iv, GDPRES); - const amrex::Real ugpx = qx(ivp0, qvidx0); - const amrex::Real ugmx = qx(iv, qvidx0); - - const amrex::Real dAupx = pggpx * ugpx - pggmx * ugmx; - const amrex::Real pavx = 0.5 * (pggpx + pggmx); - const amrex::Real dAux = ugpx - ugmx; - - const amrex::Real pggpy = qy(ivp1, GDPRES); - const amrex::Real pggmy = qy(iv, GDPRES); - const amrex::Real ugpy = qy(ivp1, qvidx1); - const amrex::Real ugmy = qy(iv, qvidx1); - - const amrex::Real dAupy = pggpy * ugpy - pggmy * ugmy; - const amrex::Real pavy = 0.5 * (pggpy + pggmy); - const amrex::Real dAuy = ugpy - ugmy; + // Update hydro vars + amrex::Real dAupx = 0.0, pavx = 0.0, dAux = 0.0; + amrex::Real dAupy = 0.0, pavy = 0.0, dAuy = 0.0; + if (no_cov_face) { + const amrex::Real pggpx = qx(ivpt0, GDPRES); + const amrex::Real pggmx = qx(iv, GDPRES); + const amrex::Real ugpx = qx(ivpt0, qvidx0); + const amrex::Real ugmx = qx(iv, qvidx0); + + dAupx = pggpx * ugpx - pggmx * ugmx; + pavx = 0.5 * (pggpx + pggmx); + dAux = ugpx - ugmx; + + const amrex::Real pggpy = qy(ivpt1, GDPRES); + const amrex::Real pggmy = qy(iv, GDPRES); + const amrex::Real ugpy = qy(ivpt1, qvidx1); + const amrex::Real ugmy = qy(iv, qvidx1); + + dAupy = pggpy * ugpy - pggmy * ugmy; + pavy = 0.5 * (pggpy + pggmy); + dAuy = ugpy - ugmy; + } const amrex::Real pxnew = cdtdx0 * (dAupx + pavx * dAux * (c - 1.0)); const amrex::Real pynew = cdtdx1 * (dAupy + pavy * dAuy * (c - 1.0)); // qp state - const amrex::Real rur = rrr * qzp(iv, QU); - const amrex::Real rvr = rrr * qzp(iv, QV); - const amrex::Real rwr = rrr * qzp(iv, QW); - const amrex::Real ekinr = 0.5 * (rur * rur + rvr * rvr + rwr * rwr) / rrr; - const amrex::Real rer = qzp(iv, QREINT) + ekinr; - - const amrex::Real runewr = rur - flxu; - const amrex::Real rvnewr = rvr - flxv; - const amrex::Real rwnewr = rwr - flxw; - const amrex::Real renewr = rer - flxe; - - qp(iv, QRHO) = rrnewr + hdt * srcrho; - qp(iv, QU) = runewr / rrnewr + hdt * srcu; - qp(iv, QV) = rvnewr / rrnewr + hdt * srcv; - qp(iv, QW) = rwnewr / rrnewr + hdt * srcw; - const amrex::Real rhoekinr = - 0.5 * (runewr * runewr + rvnewr * rvnewr + rwnewr * rwnewr) / rrnewr; - qp(iv, QREINT) = renewr - rhoekinr + hdt * srce; - amrex::Real temppres = qzp(iv, QPRES) - pxnew - pynew + hdt * srcp; - qp(iv, QPRES) = - amrex::max(temppres, std::numeric_limits::min()); + if (lo_face_not_covered) { + const amrex::Real rur = rrr * qnormp(iv, QU); + const amrex::Real rvr = rrr * qnormp(iv, QV); + const amrex::Real rwr = rrr * qnormp(iv, QW); + const amrex::Real ekinr = 0.5 * (rur * rur + rvr * rvr + rwr * rwr) / rrr; + const amrex::Real rer = qnormp(iv, QREINT) + ekinr; + + amrex::Real runewr = rur - flxu; + amrex::Real rvnewr = rvr - flxv; + amrex::Real rwnewr = rwr - flxw; + amrex::Real renewr = rer - flxe; + + qp(iv, QRHO) = rrnewr + hdt * srcrho; + qp(iv, QU) = runewr / rrnewr + hdt * srcu; + qp(iv, QV) = rvnewr / rrnewr + hdt * srcv; + qp(iv, QW) = rwnewr / rrnewr + hdt * srcw; + + const amrex::Real rhoekinr = + 0.5 * (runewr * runewr + rvnewr * rvnewr + rwnewr * rwnewr) / rrnewr; + qp(iv, QREINT) = renewr - rhoekinr + hdt * srce; + + qp(iv, QPRES) = qnormp(iv, QPRES) - pxnew - pynew + hdt * srcp; + qp(iv, QPRES) = amrex::max( + qp(iv, QPRES), std::numeric_limits::min()); + } else { + qp(iv, QRHO) = qnormp(iv, QRHO); + qp(iv, QU) = qnormp(iv, QU); + qp(iv, QV) = qnormp(iv, QV); + qp(iv, QREINT) = qnormp(iv, QREINT); + qp(iv, QPRES) = qnormp(iv, QPRES); + } // qm state - const amrex::Real rul = rrl * qzm(ivp2, QU); - const amrex::Real rvl = rrl * qzm(ivp2, QV); - const amrex::Real rwl = rrl * qzm(ivp2, QW); - const amrex::Real ekinl = 0.5 * (rul * rul + rvl * rvl + rwl * rwl) / rrl; - const amrex::Real rel = qzm(ivp2, QREINT) + ekinl; - - const amrex::Real runewl = rul - flxu; - const amrex::Real rvnewl = rvl - flxv; - const amrex::Real rwnewl = rwl - flxw; - const amrex::Real renewl = rel - flxe; - - qm(ivp2, QRHO) = rrnewl + hdt * srcrho; - qm(ivp2, QU) = runewl / rrnewl + hdt * srcu; - qm(ivp2, QV) = rvnewl / rrnewl + hdt * srcv; - qm(ivp2, QW) = rwnewl / rrnewl + hdt * srcw; - const amrex::Real rhoekinl = - 0.5 * (runewl * runewl + rvnewl * rvnewl + rwnewl * rwnewl) / rrnewl; - qm(ivp2, QREINT) = renewl - rhoekinl + hdt * srce; - temppres = qzm(ivp2, QPRES) - pxnew - pynew + hdt * srcp; - qm(ivp2, QPRES) = - amrex::max(temppres, std::numeric_limits::min()); + if (hi_face_not_covered) { + const amrex::Real rul = rrl * qnormm(ivpn, QU); + const amrex::Real rvl = rrl * qnormm(ivpn, QV); + const amrex::Real rwl = rrl * qnormm(ivpn, QW); + const amrex::Real ekinl = 0.5 * (rul * rul + rvl * rvl + rwl * rwl) / rrl; + const amrex::Real rel = qnormm(ivpn, QREINT) + ekinl; + + amrex::Real runewl = rul - flxu; + amrex::Real rvnewl = rvl - flxv; + amrex::Real rwnewl = rwl - flxw; + amrex::Real renewl = rel - flxe; + + qm(ivpn, QRHO) = rrnewl + hdt * srcrho; + qm(ivpn, QU) = runewl / rrnewl + hdt * srcu; + qm(ivpn, QV) = rvnewl / rrnewl + hdt * srcv; + qm(ivpn, QW) = rwnewl / rrnewl + hdt * srcw; + + const amrex::Real rhoekinl = + 0.5 * (runewl * runewl + rvnewl * rvnewl + rwnewl * rwnewl) / rrnewl; + qm(ivpn, QREINT) = renewl - rhoekinl + hdt * srce; + + qm(ivpn, QPRES) = qnormm(ivpn, QPRES) - pxnew - pynew + hdt * srcp; + qm(ivpn, QPRES) = amrex::max( + qm(ivpn, QPRES), std::numeric_limits::min()); + + } else { + qm(ivpn, QRHO) = qnormm(ivpn, QRHO); + qm(ivpn, QU) = qnormm(ivpn, QU); + qm(ivpn, QV) = qnormm(ivpn, QV); + qm(ivpn, QW) = qnormm(ivpn, QW); + qm(ivpn, QREINT) = qnormm(ivpn, QREINT); + qm(ivpn, QPRES) = qnormm(ivpn, QPRES); + } } // 2D version of transdd and transdo @@ -657,14 +874,16 @@ pc_transd( const int dir, amrex::Array4 const& qm, amrex::Array4 const& qp, - amrex::Array4 const& qym, - amrex::Array4 const& qyp, + amrex::Array4 const& qnormm, + amrex::Array4 const& qnormp, amrex::Array4 const& flxx, amrex::Array4 const& srcQ, amrex::Array4 const& qa, amrex::Array4 const& qint, const amrex::Real hdt, - const amrex::Real cdtdx) + const amrex::Real cdtdx, + amrex::Array4 const& norm_area = {}, + amrex::Array4 const& trans_area = {}) { const amrex::IntVect iv{AMREX_D_DECL(i, j, k)}; // ivpt is an offset by 1 in the transverse direction @@ -677,107 +896,160 @@ pc_transd( iv + amrex::IntVect::TheDimensionVector(dir == 0 ? 0 : 1)); const int qvidx = (dir == 0) ? GDV : GDU; - const amrex::Real flxrho = cdtdx * (flxx(ivpt, URHO) - flxx(iv, URHO)); - const amrex::Real flxu = cdtdx * (flxx(ivpt, UMX) - flxx(iv, UMX)); - const amrex::Real flxv = cdtdx * (flxx(ivpt, UMY) - flxx(iv, UMY)); - const amrex::Real flxe = cdtdx * (flxx(ivpt, UEDEN) - flxx(iv, UEDEN)); + bool no_cov_face = true; + bool lo_face_not_covered = true; + bool hi_face_not_covered = true; +#ifdef AMREX_USE_EB + if (trans_area) { + no_cov_face = (trans_area(ivpt) > 0.0) && (trans_area(iv) > 0.0); + } + if (norm_area) { + lo_face_not_covered = (norm_area(iv) > amrex::Real(0.0)); + hi_face_not_covered = (norm_area(ivpn) > amrex::Real(0.0)); + } +#endif + + const amrex::Real flxrho = no_cov_face + ? cdtdx * (flxx(ivpt, URHO) - flxx(iv, URHO)) + : amrex::Real(0.0); + const amrex::Real flxu = + no_cov_face ? cdtdx * (flxx(ivpt, UMX) - flxx(iv, UMX)) : amrex::Real(0.0); + const amrex::Real flxv = + no_cov_face ? cdtdx * (flxx(ivpt, UMY) - flxx(iv, UMY)) : amrex::Real(0.0); + const amrex::Real flxe = no_cov_face + ? cdtdx * (flxx(ivpt, UEDEN) - flxx(iv, UEDEN)) + : amrex::Real(0.0); const amrex::Real srcr = srcQ(iv, QRHO); const amrex::Real srce = srcQ(iv, QREINT); const amrex::Real srcp = srcQ(iv, QPRES); - const amrex::Real rrr = qyp(iv, QRHO); - const amrex::Real rrl = qym(ivpn, QRHO); const amrex::Real c = qa(iv, QGAMC); // Update passive variables #if AMREX_SPACEDIM == 2 pc_transd_passive( - iv, ivpt, ivpn, 0, UMZ, QW, cdtdx, hdt, flxrho, rrr, rrl, srcQ, qyp, qym, - flxx, qp, qm); + iv, ivpn, ivpt, 0, UMZ, QW, cdtdx, hdt, flxrho, true, srcQ, qnormp, qnormm, + flxx, qp, qm, no_cov_face, lo_face_not_covered, hi_face_not_covered); #endif for (int n = 0; n < NUM_ADV; n++) { pc_transd_passive( - iv, ivpt, ivpn, n, UFA, QFA, cdtdx, hdt, flxrho, rrr, rrl, srcQ, qyp, qym, - flxx, qp, qm); + iv, ivpn, ivpt, n, UFA, QFA, cdtdx, hdt, flxrho, true, srcQ, qnormp, + qnormm, flxx, qp, qm, no_cov_face, lo_face_not_covered, + hi_face_not_covered); } for (int n = 0; n < NUM_SPECIES; n++) { pc_transd_passive( - iv, ivpt, ivpn, n, UFS, QFS, cdtdx, hdt, flxrho, rrr, rrl, srcQ, qyp, qym, - flxx, qp, qm); + iv, ivpn, ivpt, n, UFS, QFS, cdtdx, hdt, flxrho, true, srcQ, qnormp, + qnormm, flxx, qp, qm, no_cov_face, lo_face_not_covered, + hi_face_not_covered); } for (int n = 0; n < NUM_AUX; n++) { pc_transd_passive( - iv, ivpt, ivpn, n, UFX, QFX, cdtdx, hdt, flxrho, rrr, rrl, srcQ, qyp, qym, - flxx, qp, qm); + iv, ivpn, ivpt, n, UFX, QFX, cdtdx, hdt, flxrho, true, srcQ, qnormp, + qnormm, flxx, qp, qm, no_cov_face, lo_face_not_covered, + hi_face_not_covered); } for (int n = 0; n < NUM_LIN; n++) { pc_transd_passive( - iv, ivpt, ivpn, n, ULIN, QLIN, cdtdx, hdt, 0., 1., 1., srcQ, qyp, qym, - flxx, qp, qm); + iv, ivpn, ivpt, n, ULIN, QLIN, cdtdx, hdt, 0., false, srcQ, qnormp, + qnormm, flxx, qp, qm, no_cov_face, lo_face_not_covered, + hi_face_not_covered); } - const amrex::Real pggp = qint(ivpt, GDPRES); - const amrex::Real pggm = qint(iv, GDPRES); - const amrex::Real ugp = qint(ivpt, qvidx); - const amrex::Real ugm = qint(iv, qvidx); - const amrex::Real dAup = pggp * ugp - pggm * ugm; - const amrex::Real pav = 0.5 * (pggp + pggm); - const amrex::Real dAu = ugp - ugm; + amrex::Real dAup = 0.0, pav = 0.0, dAu = 0.0; + if (no_cov_face) { + const amrex::Real pggp = qint(ivpt, GDPRES); + const amrex::Real pggm = qint(iv, GDPRES); + const amrex::Real ugp = qint(ivpt, qvidx); + const amrex::Real ugm = qint(iv, qvidx); + + dAup = pggp * ugp - pggm * ugm; + pav = 0.5 * (pggp + pggm); + dAu = ugp - ugm; + } // QP - // Convert to conservative - amrex::Real rur = qyp(iv, QU); - amrex::Real rvr = qyp(iv, QV); - const amrex::Real ekinr = 0.5 * rrr * (rur * rur + rvr * rvr); - rur *= rrr; - rvr *= rrr; + if (lo_face_not_covered) { + // Convert to conservative + const amrex::Real rrr = qnormp(iv, QRHO); + amrex::Real rur = qnormp(iv, QU); + amrex::Real rvr = qnormp(iv, QV); + const amrex::Real ekinr = 0.5 * rrr * (rur * rur + rvr * rvr); + rur *= rrr; + rvr *= rrr; + + const amrex::Real rer = qnormp(iv, QREINT) + ekinr; + // Add transverse predictor + amrex::Real rrnewr = rrr - flxrho; + amrex::Real runewr = rur - flxu; + amrex::Real rvnewr = rvr - flxv; + amrex::Real renewr = rer - flxe; + + // Convert back to primitive + qp(iv, QRHO) = rrnewr + hdt * srcr; + qp(iv, QU) = runewr / rrnewr + hdt * srcQ(iv, QU); + qp(iv, QV) = rvnewr / rrnewr + hdt * srcQ(iv, QV); + + const amrex::Real rhoekinr = + 0.5 * (runewr * runewr + rvnewr * rvnewr) / rrnewr; + qp(iv, QREINT) = renewr - rhoekinr + hdt * srce; + + amrex::Real pnewr = + qnormp(iv, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.)); + pnewr += hdt * srcp; + qp(iv, QPRES) = pnewr; + + qp(iv, QPRES) = amrex::max( + qp(iv, QPRES), std::numeric_limits::min()); - const amrex::Real rer = qyp(iv, QREINT) + ekinr; - // Add transverse predictor - const amrex::Real rrnewr = rrr - flxrho; - const amrex::Real runewr = rur - flxu; - const amrex::Real rvnewr = rvr - flxv; - const amrex::Real renewr = rer - flxe; - - // Convert back to primitive - qp(iv, QRHO) = rrnewr + hdt * srcr; - qp(iv, QU) = runewr / rrnewr + hdt * srcQ(iv, QU); - qp(iv, QV) = rvnewr / rrnewr + hdt * srcQ(iv, QV); - const amrex::Real rhoekinr = - 0.5 * (runewr * runewr + rvnewr * rvnewr) / rrnewr; - amrex::Real pnewr = qyp(iv, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.)); - pnewr += hdt * srcp; - qp(iv, QPRES) = - amrex::max(pnewr, std::numeric_limits::min()); - qp(iv, QREINT) = renewr - rhoekinr + hdt * srce; + } else { + qp(iv, QRHO) = qnormp(iv, QRHO); + qp(iv, QU) = qnormp(iv, QU); + qp(iv, QV) = qnormp(iv, QV); + qp(iv, QREINT) = qnormp(iv, QREINT); + qp(iv, QPRES) = qnormp(iv, QPRES); + } // QM - // Conversion to Conservative - amrex::Real rul = qym(ivpn, QU); - amrex::Real rvl = qym(ivpn, QV); - const amrex::Real ekinl = 0.5 * rrl * (rul * rul + rvl * rvl); - rul *= rrl; - rvl *= rrl; - const amrex::Real rel = qym(ivpn, QREINT) + ekinl; + if (hi_face_not_covered) { + // Conversion to Conservative + const amrex::Real rrl = qnormm(ivpn, QRHO); + amrex::Real rul = qnormm(ivpn, QU); + amrex::Real rvl = qnormm(ivpn, QV); + const amrex::Real ekinl = 0.5 * rrl * (rul * rul + rvl * rvl); + rul *= rrl; + rvl *= rrl; + const amrex::Real rel = qnormm(ivpn, QREINT) + ekinl; + + // Transverse fluxes + amrex::Real rrnewl = rrl - flxrho; + amrex::Real runewl = rul - flxu; + amrex::Real rvnewl = rvl - flxv; + amrex::Real renewl = rel - flxe; + + qm(ivpn, QRHO) = rrnewl + hdt * srcr; + qm(ivpn, QU) = runewl / rrnewl + hdt * srcQ(iv, QU); + qm(ivpn, QV) = rvnewl / rrnewl + hdt * srcQ(iv, QV); + + const amrex::Real rhoekinl = + 0.5 * (runewl * runewl + rvnewl * rvnewl) / rrnewl; + qm(ivpn, QREINT) = renewl - rhoekinl + hdt * srce; + + amrex::Real pnewl = + qnormm(ivpn, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.)); + qm(ivpn, QPRES) = pnewl + hdt * srcp; + + qm(ivpn, QPRES) = amrex::max( + qm(ivpn, QPRES), std::numeric_limits::min()); - // Transverse fluxes - const amrex::Real rrnewl = rrl - flxrho; - const amrex::Real runewl = rul - flxu; - const amrex::Real rvnewl = rvl - flxv; - const amrex::Real renewl = rel - flxe; - - qm(ivpn, QRHO) = rrnewl + hdt * srcr; - qm(ivpn, QU) = runewl / rrnewl + hdt * srcQ(iv, QU); - qm(ivpn, QV) = rvnewl / rrnewl + hdt * srcQ(iv, QV); - const amrex::Real rhoekinl = - 0.5 * (runewl * runewl + rvnewl * rvnewl) / rrnewl; - - qm(ivpn, QREINT) = renewl - rhoekinl + hdt * srce; - const amrex::Real pnewl = - qym(ivpn, QPRES) - cdtdx * (dAup + pav * dAu * (c - 1.)); - qm(ivpn, QPRES) = amrex::max( - pnewl + hdt * srcp, std::numeric_limits::min()); + } else { + qm(ivpn, QRHO) = qnormm(ivpn, QRHO); + qm(ivpn, QU) = qnormm(ivpn, QU); + qm(ivpn, QV) = qnormm(ivpn, QV); + qm(ivpn, QREINT) = qnormm(ivpn, QREINT); + qm(ivpn, QPRES) = qnormm(ivpn, QPRES); + } } // Use interface states from Riemann solver for pdivu. @@ -825,7 +1097,11 @@ pc_artif_visc( amrex::Array4 const& u, amrex::Real const dx, amrex::Real const difmag, - const int dir) + const int dir, + const int domlo, + const int domhi, + const int bclo, + const int bchi) { const amrex::IntVect iv{AMREX_D_DECL(i, j, k)}; const amrex::IntVect ivm(iv - amrex::IntVect::TheDimensionVector(dir)); @@ -838,12 +1114,12 @@ pc_artif_visc( bdim[0] * 2 + bdim[1] * 2 + bdim[2] * 1}}; const amrex::IntVect ivpj(iv + amrex::IntVect::TheDimensionVector(l_idx[1])); -#if AMREX_SPACEDIM == 3 - const amrex::IntVect ivpk(iv + amrex::IntVect::TheDimensionVector(l_idx[2])); - const amrex::IntVect ivpp( - iv + amrex::IntVect::TheDimensionVector(l_idx[1]) + - amrex::IntVect::TheDimensionVector(l_idx[2])); -#endif + AMREX_D_TERM(, , + const amrex::IntVect ivpk( + iv + amrex::IntVect::TheDimensionVector(l_idx[2])); + const amrex::IntVect ivpp( + iv + amrex::IntVect::TheDimensionVector(l_idx[1]) + + amrex::IntVect::TheDimensionVector(l_idx[2]));); const amrex::Real div = difmag * amrex::min( @@ -851,9 +1127,27 @@ pc_artif_visc( divu(iv), 0.5 * (divu(iv) + divu(ivpj)), 0.25 * (divu(iv) + divu(ivpj) + divu(ivpk) + divu(ivpp)))); + bool at_bndry = ((dir == 0) && (i == domlo) && (bclo == NoSlipWall)) || + ((dir == 0) && (i == domhi + 1) && (bchi == NoSlipWall)) || + ((dir == 1) && (j == domlo) && (bclo == NoSlipWall)) || + ((dir == 1) && (j == domhi + 1) && (bchi == NoSlipWall)); +#if (AMREX_SPACEDIM == 3) + at_bndry = at_bndry || ((dir == 2) && (k == domlo) && (bclo == NoSlipWall)) || + ((dir == 2) && (k == domhi + 1) && (bchi == NoSlipWall)); +#endif + for (int n = 0; n < NVAR; ++n) { if (n != UTEMP) { - flx(iv, n) += dx * div * (u(iv, n) - u(ivm, n)); +#if (AMREX_SPACEDIM == 2) + bool is_tang_vel = (dir == 0 && n == UMY) || (dir == 1 && n == UMX); +#elif (AMREX_SPACEDIM == 3) + bool is_tang_vel = (dir == 0 && (n == UMY || n == UMZ)) || + (dir == 1 && (n == UMX || n == UMZ)) || + (dir == 2 && (n == UMX || n == UMY)); +#endif + if (!(at_bndry && is_tang_vel)) { + flx(iv, n) += dx * div * (u(iv, n) - u(ivm, n)); + } } } flx(iv, UTEMP) = 0.0; @@ -883,6 +1177,25 @@ void pc_umeth_3D( const bool use_hybrid_weno, const int weno_scheme); +void pc_umeth_eb_3D( + amrex::Box const& bx_to_fill, + const int* bclo, + const int* bchi, + const int* domlo, + const int* domhi, + amrex::Array4 const& q, + amrex::Array4 const& qaux, + amrex::Array4 const& srcQ, + const amrex::GpuArray, 3>& flx, + const amrex::GpuArray, 3>& qec, + const amrex::GpuArray, 3>& ap, + amrex::Array4 const& flag_arr, + const amrex::GpuArray del, + const amrex::Real dt, + const int ppm_type, + const bool use_flattening, + const int plm_iorder); + #elif AMREX_SPACEDIM == 2 void pc_umeth_2D( @@ -906,6 +1219,25 @@ void pc_umeth_2D( const bool use_flattening, const bool use_hybrid_weno, const int weno_scheme); + +void pc_umeth_eb_2D( + amrex::Box const& bx_to_fill, + const int* bclo, + const int* bchi, + const int* domlo, + const int* domhi, + amrex::Array4 const& q, + amrex::Array4 const& qaux, + amrex::Array4 const& srcQ, + const amrex::GpuArray, 2>& flx, + const amrex::GpuArray, 2>& qec, + const amrex::GpuArray, 2>& ap, + amrex::Array4 const& flag_arr, + const amrex::GpuArray del, + const amrex::Real dt, + const int ppm_type, + const bool use_flattening, + const int plm_iorder); #endif #endif diff --git a/Source/Godunov.cpp b/Source/Godunov.cpp index 0c72b25f0..d121b36c1 100644 --- a/Source/Godunov.cpp +++ b/Source/Godunov.cpp @@ -11,6 +11,7 @@ pc_low_order_boundary( const int domlo, const int domhi, const int plm_iorder, + const bool use_flattening, const int idir, const amrex::Real dx, const amrex::Real dt, @@ -28,8 +29,18 @@ pc_low_order_boundary( auto const& qbparr = qbp.array(); amrex::ParallelFor(bdbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real slope[QVAR]; + + amrex::Real flat = 1.0; + // Calculate flattening in-place + if (use_flattening) { + for (int dir_flat = 0; dir_flat < AMREX_SPACEDIM; dir_flat++) { + flat = amrex::min( + flat, flatten(AMREX_D_DECL(i, j, k), dir_flat, q)); + } + } + for (int n = 0; n < QVAR; ++n) { - slope[n] = plm_slope(AMREX_D_DECL(i, j, k), n, idir, q, plm_iorder); + slope[n] = plm_slope(AMREX_D_DECL(i, j, k), n, idir, q, flat, plm_iorder); } pc_plm_d( AMREX_D_DECL(i, j, k), idir, qbmarr, qbparr, slope, q, qa(i, j, k, QC), @@ -123,16 +134,25 @@ pc_umeth_3D( auto const& qzparr = qzp.array(); // Put the PLM and slopes in the same kernel launch to avoid unnecessary - // launch overhead Pelec_Slope_* are SIMD as well as PeleC_plm_* which loop - // over the same box + // launch if (ppm_type == 0) { amrex::ParallelFor( bxg2, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real slope[QVAR]; + + amrex::Real flat = 1.0; + // Calculate flattening in-place + if (use_flattening) { + for (int dir_flat = 0; dir_flat < AMREX_SPACEDIM; dir_flat++) { + flat = amrex::min(flat, flatten(i, j, k, dir_flat, q)); + } + } + // X slopes and interp int idir = 0; for (int n = 0; n < QVAR; ++n) { - slope[n] = plm_slope(AMREX_D_DECL(i, j, k), n, 0, q, plm_iorder); + slope[n] = + plm_slope(AMREX_D_DECL(i, j, k), n, 0, q, flat, plm_iorder); } pc_plm_d( AMREX_D_DECL(i, j, k), idir, qxmarr, qxparr, slope, q, @@ -141,7 +161,8 @@ pc_umeth_3D( // Y slopes and interp idir = 1; for (int n = 0; n < QVAR; n++) { - slope[n] = plm_slope(AMREX_D_DECL(i, j, k), n, 1, q, plm_iorder); + slope[n] = + plm_slope(AMREX_D_DECL(i, j, k), n, 1, q, flat, plm_iorder); } pc_plm_d( AMREX_D_DECL(i, j, k), idir, qymarr, qyparr, slope, q, @@ -150,7 +171,8 @@ pc_umeth_3D( // Z slopes and interp idir = 2; for (int n = 0; n < QVAR; ++n) { - slope[n] = plm_slope(AMREX_D_DECL(i, j, k), n, 2, q, plm_iorder); + slope[n] = + plm_slope(AMREX_D_DECL(i, j, k), n, 2, q, flat, plm_iorder); } pc_plm_d( AMREX_D_DECL(i, j, k), idir, qzmarr, qzparr, slope, q, @@ -428,7 +450,7 @@ pc_umeth_3D( if (bfbx.ok()) { pc_low_order_boundary( bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, - idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); } } if ( @@ -440,7 +462,7 @@ pc_umeth_3D( if (bfbx.ok()) { pc_low_order_boundary( bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, - idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); } } } @@ -453,6 +475,539 @@ pc_umeth_3D( }); } +void +pc_umeth_eb_3D( + amrex::Box const& bx_to_fill, + const int* bclo, + const int* bchi, + const int* domlo, + const int* domhi, + amrex::Array4 const& q, + amrex::Array4 const& qaux, + amrex::Array4 const& srcQ, + const amrex::GpuArray, 3>& flx, + const amrex::GpuArray, 3>& qec, + const amrex::GpuArray, 3>& ap, + amrex::Array4 const& flag_arr, + const amrex::GpuArray del, + const amrex::Real dt, + const int ppm_type, + const bool use_flattening, + const int plm_iorder) +{ + int cdir; + + amrex::Real const dx = del[0]; + amrex::Real const dy = del[1]; + amrex::Real const dz = del[2]; + + amrex::Real const hdt = amrex::Real(0.5) * dt; + amrex::Real const hdtdx = hdt / dx; + amrex::Real const hdtdy = hdt / dy; + amrex::Real const hdtdz = hdt / dz; + amrex::Real const cdtdx = amrex::Real(1.0 / 3.0) * dt / dx; + amrex::Real const cdtdy = amrex::Real(1.0 / 3.0) * dt / dy; + amrex::Real const cdtdz = amrex::Real(1.0 / 3.0) * dt / dz; + + const int bclx = bclo[0]; + const int bcly = bclo[1]; + const int bclz = bclo[2]; + const int bchx = bchi[0]; + const int bchy = bchi[1]; + const int bchz = bchi[2]; + const int dlx = domlo[0]; + const int dly = domlo[1]; + const int dlz = domlo[2]; + const int dhx = domhi[0]; + const int dhy = domhi[1]; + const int dhz = domhi[2]; + + const amrex::Box& bxg1 = grow(bx_to_fill, 1); + const amrex::Box& bxg2 = grow(bx_to_fill, 2); + + // X data + cdir = 0; + const amrex::Box& xmbx = growHi(bxg2, cdir, 1); + amrex::FArrayBox qxm(xmbx, QVAR, amrex::The_Async_Arena()); + auto const& qxmarr = qxm.array(); + amrex::FArrayBox qxp(bxg2, QVAR, amrex::The_Async_Arena()); + auto const& qxparr = qxp.array(); + + // Y data + cdir = 1; + const amrex::Box& ymbx = growHi(bxg2, cdir, 1); + amrex::FArrayBox qym(ymbx, QVAR, amrex::The_Async_Arena()); + auto const& qymarr = qym.array(); + amrex::FArrayBox qyp(bxg2, QVAR, amrex::The_Async_Arena()); + auto const& qyparr = qyp.array(); + + // Z data + cdir = 2; + const amrex::Box& zmbx = growHi(bxg2, cdir, 1); + amrex::FArrayBox qzm(zmbx, QVAR, amrex::The_Async_Arena()); + auto const& qzmarr = qzm.array(); + amrex::FArrayBox qzp(bxg2, QVAR, amrex::The_Async_Arena()); + auto const& qzparr = qzp.array(); + // + // Put the PLM and slopes in the same kernel launch to avoid unnecessary + // launch overhead Note that we compute the qm and qp values on bxg2 -5,-5,-5 + // + if (ppm_type == 0) { + amrex::ParallelFor( + bxg2, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (flag_arr(i, j, k).isCovered()) { + return; + } + + amrex::Real slope[QVAR]; + + amrex::Real flat = 1.0; + // Calculate flattening in-place + if (use_flattening) { + for (int dir_flat = 0; dir_flat < AMREX_SPACEDIM; dir_flat++) { + flat = amrex::min( + flat, flatten_eb(i, j, k, dir_flat, flag_arr, q)); + } + } + + // X slopes and interp + int idir = 0; + for (int n = 0; n < QVAR; ++n) { + slope[n] = plm_slope_eb( + AMREX_D_DECL(i, j, k), n, 0, flag_arr, q, flat, plm_iorder); + } + // cells -5,-5,-5 + pc_plm_d_eb( + i, j, k, idir, qxmarr, qxparr, slope, q, qaux(i, j, k, QC), dx, dt, + ap[idir]); + + // Y slopes and interp + idir = 1; + for (int n = 0; n < QVAR; n++) { + slope[n] = plm_slope_eb( + AMREX_D_DECL(i, j, k), n, 1, flag_arr, q, flat, plm_iorder); + } + // cells -5,-5,-5 + pc_plm_d_eb( + i, j, k, idir, qymarr, qyparr, slope, q, qaux(i, j, k, QC), dy, dt, + ap[idir]); + + // Z slopes and interp + idir = 2; + for (int n = 0; n < QVAR; ++n) { + slope[n] = plm_slope_eb( + AMREX_D_DECL(i, j, k), n, 2, flag_arr, q, flat, plm_iorder); + } + // cells -5,-5,-5 + pc_plm_d_eb( + i, j, k, idir, qzmarr, qzparr, slope, q, qaux(i, j, k, QC), dz, dt, + ap[idir]); + }); + } else { + amrex::Error("ppm_type must be 0 (PLM) when using EB"); + } + + // These are the first flux estimates as per the corner-transport-upwind + // method + + // ************************************************************************************* + // X initial fluxes + // ************************************************************************************* + cdir = 0; + const amrex::Box& xflxbx = surroundingNodes(grow(bxg2, cdir, -1), cdir); + + amrex::FArrayBox fx(xflxbx, NVAR, amrex::The_Async_Arena()); + auto const& fxarr = fx.array(); + + amrex::FArrayBox qgdx(xflxbx, NGDNV, amrex::The_Async_Arena()); + auto const& gdtempx = qgdx.array(); + + // -4,-5,-5 + amrex::ParallelFor( + xflxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclx, bchx, dlx, dhx, qxmarr, qxparr, fxarr, gdtempx, qaux, + cdir); + } + }); + + // ************************************************************************************* + // Y initial fluxes + // ************************************************************************************* + cdir = 1; + const amrex::Box& yflxbx = surroundingNodes(grow(bxg2, cdir, -1), cdir); + + amrex::FArrayBox fy(yflxbx, NVAR, amrex::The_Async_Arena()); + auto const& fyarr = fy.array(); + + amrex::FArrayBox qgdy(yflxbx, NGDNV, amrex::The_Async_Arena()); + auto const& gdtempy = qgdy.array(); + + // -5,-4,-5 + amrex::ParallelFor( + yflxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bcly, bchy, dly, dhy, qymarr, qyparr, fyarr, gdtempy, qaux, + cdir); + } + }); + + // ************************************************************************************* + // Z initial fluxes + // ************************************************************************************* + cdir = 2; + const amrex::Box& zflxbx = surroundingNodes(grow(bxg2, cdir, -1), cdir); + + amrex::FArrayBox fz(zflxbx, NVAR, amrex::The_Async_Arena()); + auto const& fzarr = fz.array(); + + amrex::FArrayBox qgdz(zflxbx, NGDNV, amrex::The_Async_Arena()); + auto const& gdtempz = qgdz.array(); + + // -5,-5,-4 + amrex::ParallelFor( + zflxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclz, bchz, dlz, dhz, qzmarr, qzparr, fzarr, gdtempz, qaux, + cdir); + } + }); + + // ************************************************************************************* + // X interface corrections + // ************************************************************************************* + cdir = 0; + amrex::FArrayBox qxym(xmbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qxyp(xmbx, QVAR, amrex::The_Async_Arena()); + auto const& qmxy = qxym.array(); + auto const& qpxy = qxyp.array(); + + amrex::FArrayBox qxzm(xmbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qxzp(xmbx, QVAR, amrex::The_Async_Arena()); + auto const& qmxz = qxzm.array(); + auto const& qpxz = qxzp.array(); + + // + amrex::Box xybx(bxg2); + xybx.grow(1, -1); + amrex::ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // fyarr was made on y-faces to -5,-4,-5 + // this loop is over cells to -5,-4,-5 + if (!flag_arr(i, j, k).isCovered()) { + // X|Y + pc_transdo( + i, j, k, cdir, 1, qmxy, qpxy, qxmarr, qxparr, fyarr, qaux, gdtempy, + cdtdy, ap[cdir], ap[1]); + } + }); + + amrex::Box xzbx(bxg2); + xzbx.grow(2, -1); + amrex::ParallelFor(xzbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // fzarr was made on z-faces to -5,-5,-4 + // this loop is over cells to -5,-5,-4 + if (!flag_arr(i, j, k).isCovered()) { + // X|Z + pc_transdo( + i, j, k, cdir, 2, qmxz, qpxz, qxmarr, qxparr, fzarr, qaux, gdtempz, + cdtdz, ap[cdir], ap[2]); + } + }); + + const amrex::Box& txfxbx = surroundingNodes(bxg2, cdir); + amrex::FArrayBox fluxxy(txfxbx, NVAR, amrex::The_Async_Arena()); + amrex::FArrayBox fluxxz(txfxbx, NVAR, amrex::The_Async_Arena()); + amrex::FArrayBox gdvxyfab(txfxbx, NGDNV, amrex::The_Async_Arena()); + amrex::FArrayBox gdvxzfab(txfxbx, NGDNV, amrex::The_Async_Arena()); + + auto const& flxy = fluxxy.array(); + auto const& flxz = fluxxz.array(); + auto const& qxy = gdvxyfab.array(); + auto const& qxz = gdvxzfab.array(); + + // Riemann problem X|Y + amrex::Box xycmpbx(surroundingNodes(bxg1, 0).grow(2, 1)); + // this loop is over x-faces to -4,-4,-5 + amrex::ParallelFor( + xycmpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // X|Y + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclx, bchx, dlx, dhx, qmxy, qpxy, flxy, qxy, qaux, cdir); + } + }); + + // Riemann problem X|Z + amrex::Box xzcmpbx(surroundingNodes(bxg1, 0).grow(1, 1)); + // this loop is over x-faces to -4,-5,-4 + amrex::ParallelFor( + xzcmpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // X|Z + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclx, bchx, dlx, dhx, qmxz, qpxz, flxz, qxz, qaux, cdir); + } + }); + + // ************************************************************************************* + // Y interface corrections + // ************************************************************************************* + + cdir = 1; + amrex::FArrayBox qyxm(ymbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qyxp(ymbx, QVAR, amrex::The_Async_Arena()); + auto const& qmyx = qyxm.array(); + auto const& qpyx = qyxp.array(); + + amrex::FArrayBox qyzm(ymbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qyzp(ymbx, QVAR, amrex::The_Async_Arena()); + auto const& qmyz = qyzm.array(); + auto const& qpyz = qyzp.array(); + + amrex::Box yxbx(bxg2); + yxbx.grow(0, -1); + amrex::ParallelFor(yxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Y|X + if (!flag_arr(i, j, k).isCovered()) { + pc_transdo( + i, j, k, cdir, 0, qmyx, qpyx, qymarr, qyparr, fxarr, qaux, gdtempx, + cdtdx, ap[cdir], ap[0]); + } + }); + + amrex::Box yzbx(bxg2); + yzbx.grow(2, -1); + amrex::ParallelFor(yzbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Y|Z + if (!flag_arr(i, j, k).isCovered()) { + pc_transdo( + i, j, k, cdir, 2, qmyz, qpyz, qymarr, qyparr, fzarr, qaux, gdtempz, + cdtdz, ap[cdir], ap[2]); + } + }); + + // Riemann problem Y|X Y|Z + const amrex::Box& tyfxbx = surroundingNodes(bxg2, cdir); + amrex::FArrayBox fluxyx(tyfxbx, NVAR, amrex::The_Async_Arena()); + amrex::FArrayBox fluxyz(tyfxbx, NVAR, amrex::The_Async_Arena()); + amrex::FArrayBox gdvyxfab(tyfxbx, NGDNV, amrex::The_Async_Arena()); + amrex::FArrayBox gdvyzfab(tyfxbx, NGDNV, amrex::The_Async_Arena()); + + auto const& flyx = fluxyx.array(); + auto const& flyz = fluxyz.array(); + auto const& qyx = gdvyxfab.array(); + auto const& qyz = gdvyzfab.array(); + + // Riemann problem Y|X + amrex::Box yxcmpbx(surroundingNodes(bxg1, 1).grow(2, 1)); + amrex::ParallelFor( + yxcmpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Y|X + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bcly, bchy, dly, dhy, qmyx, qpyx, flyx, qyx, qaux, cdir); + } + }); + + // Riemann problem Y|Z + amrex::Box yzcmpbx(surroundingNodes(bxg1, 1).grow(0, 1)); + amrex::ParallelFor( + yzcmpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Y|Z + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bcly, bchy, dly, dhy, qmyz, qpyz, flyz, qyz, qaux, cdir); + } + }); + + // ************************************************************************************* + // Z interface corrections + // ************************************************************************************* + cdir = 2; + + amrex::FArrayBox qzxm(zmbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qzxp(zmbx, QVAR, amrex::The_Async_Arena()); + auto const& qmzx = qzxm.array(); + auto const& qpzx = qzxp.array(); + + amrex::FArrayBox qzym(zmbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qzyp(zmbx, QVAR, amrex::The_Async_Arena()); + auto const& qmzy = qzym.array(); + auto const& qpzy = qzyp.array(); + + amrex::Box zxbx(bxg2); + zxbx.grow(0, -1); + amrex::ParallelFor(zxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Z|X + if (!flag_arr(i, j, k).isCovered()) { + pc_transdo( + i, j, k, cdir, 0, qmzx, qpzx, qzmarr, qzparr, fxarr, qaux, gdtempx, + cdtdx, ap[cdir], ap[0]); + } + }); + + amrex::Box zybx(bxg2); + zybx.grow(1, -1); + amrex::ParallelFor(zybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Z|Y + if (!flag_arr(i, j, k).isCovered()) { + pc_transdo( + i, j, k, cdir, 1, qmzy, qpzy, qzmarr, qzparr, fyarr, qaux, gdtempy, + cdtdy, ap[cdir], ap[1]); + } + }); + + // Riemann problem Z|X Z|Y + const amrex::Box& tzfxbx = surroundingNodes(bxg2, cdir); + amrex::FArrayBox fluxzx(tzfxbx, NVAR, amrex::The_Async_Arena()); + amrex::FArrayBox fluxzy(tzfxbx, NVAR, amrex::The_Async_Arena()); + amrex::FArrayBox gdvzxfab(tzfxbx, NGDNV, amrex::The_Async_Arena()); + amrex::FArrayBox gdvzyfab(tzfxbx, NGDNV, amrex::The_Async_Arena()); + + auto const& flzx = fluxzx.array(); + auto const& flzy = fluxzy.array(); + auto const& qzx = gdvzxfab.array(); + auto const& qzy = gdvzyfab.array(); + + // Riemann problem Z|X + amrex::Box zxcmpbx(surroundingNodes(bxg1, 2).grow(1, 1)); + amrex::ParallelFor( + zxcmpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Z|X + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclz, bchz, dlz, dhz, qmzx, qpzx, flzx, qzx, qaux, cdir); + } + }); + + // Riemann problem Z|Y + amrex::Box zycmpbx(surroundingNodes(bxg1, 2).grow(0, 1)); + amrex::ParallelFor( + zycmpbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Z|Y + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclz, bchz, dlz, dhz, qmzy, qpzy, flzy, qzy, qaux, cdir); + } + }); + + amrex::FArrayBox qmfab(bxg2, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qpfab(bxg1, QVAR, amrex::The_Async_Arena()); + auto const& qm = qmfab.array(); + auto const& qp = qpfab.array(); + + // ************************************************************************************* + // Final X steps + // ************************************************************************************* + cdir = 0; + + // this loop is over cells to -4,-4,-4 + amrex::ParallelFor(bxg1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (!flag_arr(i, j, k).isCovered()) { + pc_transdd( + i, j, k, cdir, qm, qp, qxmarr, qxparr, flyz, flzy, qyz, qzy, qaux, srcQ, + hdt, hdtdy, hdtdz, ap[0], ap[1], ap[2]); + } + }); + + // This box must be grown by one in transverse direction so we can + // do tangential interpolation when taking divergence later + amrex::Box xfbx(surroundingNodes(bx_to_fill, 0)); + xfbx.grow(1, 1); + xfbx.grow(2, 1); + amrex::ParallelFor(xfbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclx, bchx, dlx, dhx, qm, qp, flx[cdir], qec[cdir], qaux, + cdir); + } + }); + + // ************************************************************************************* + // Final Y steps + // ************************************************************************************* + cdir = 1; + // this loop is over cells to -4,-4,-4 + amrex::ParallelFor(bxg1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (!flag_arr(i, j, k).isCovered()) { + pc_transdd( + i, j, k, cdir, qm, qp, qymarr, qyparr, flxz, flzx, qxz, qzx, qaux, srcQ, + hdt, hdtdx, hdtdz, ap[1], ap[0], ap[2]); + } + }); + + // This box must be grown by one in transverse direction so we can + // do tangential interpolation when taking divergence later + amrex::Box yfbx(surroundingNodes(bx_to_fill, 1)); + yfbx.grow(0, 1); + yfbx.grow(2, 1); + amrex::ParallelFor(yfbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bcly, bchy, dly, dhy, qm, qp, flx[cdir], qec[cdir], qaux, + cdir); + } + }); + + // ************************************************************************************* + // Final Z steps + // ************************************************************************************* + cdir = 2; + amrex::ParallelFor(bxg1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (!flag_arr(i, j, k).isCovered()) { + pc_transdd( + i, j, k, cdir, qm, qp, qzmarr, qzparr, flxy, flyx, qxy, qyx, qaux, srcQ, + hdt, hdtdx, hdtdy, ap[2], ap[0], ap[1]); + } + }); + + // This box must be grown by one in transverse direction so we can + // do tangential interpolation when taking divergence later + amrex::Box zfbx(surroundingNodes(bx_to_fill, 2)); + zfbx.grow(0, 1); + zfbx.grow(1, 1); + amrex::ParallelFor(zfbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclz, bchz, dlz, dhz, qm, qp, flx[cdir], qec[cdir], qaux, + cdir); + } + }); + + // Fix bcnormal boundaries - always use PLM and don't do N+1/2 predictor + // because the user specifies conditions at N + // bcnormal used for "Hard" (inflow) and "UserBC" (user_bc) + for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { + if ( + bclo[idir] == PCPhysBCType::inflow || + bclo[idir] == PCPhysBCType::user_bc) { + // Box for fluxes at this boundary + amrex::Box bfbx = surroundingNodes(bx_to_fill, idir); + bfbx.setBig(idir, domlo[idir]); + if (bfbx.ok()) { + pc_low_order_boundary( + bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + } + } + if ( + bchi[idir] == PCPhysBCType::inflow || + bchi[idir] == PCPhysBCType::user_bc) { + // Box for fluxes at this boundary + amrex::Box bfbx = surroundingNodes(bx_to_fill, idir); + bfbx.setSmall(idir, domhi[idir] + 1); + if (bfbx.ok()) { + pc_low_order_boundary( + bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + } + } + } +} + #elif AMREX_SPACEDIM == 2 void @@ -518,16 +1073,28 @@ pc_umeth_2D( amrex::ParallelFor( bxg2, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real slope[QVAR]; + + amrex::Real flat = 1.0; + // Calculate flattening in-place + if (use_flattening) { + for (int dir_flat = 0; dir_flat < AMREX_SPACEDIM; dir_flat++) { + flat = amrex::min( + flat, flatten(AMREX_D_DECL(i, j, k), dir_flat, q)); + } + } + // X slopes and interp for (int n = 0; n < QVAR; ++n) - slope[n] = plm_slope(AMREX_D_DECL(i, j, k), n, 0, q, plm_iorder); + slope[n] = + plm_slope(AMREX_D_DECL(i, j, k), n, 0, q, flat, plm_iorder); pc_plm_d( AMREX_D_DECL(i, j, k), 0, qxmarr, qxparr, slope, q, qaux(i, j, k, QC), dx, dt); // Y slopes and interp for (int n = 0; n < QVAR; n++) - slope[n] = plm_slope(AMREX_D_DECL(i, j, k), n, 1, q, plm_iorder); + slope[n] = + plm_slope(AMREX_D_DECL(i, j, k), n, 1, q, flat, plm_iorder); pc_plm_d( AMREX_D_DECL(i, j, k), 1, qymarr, qyparr, slope, q, qaux(i, j, k, QC), dy, dt); @@ -632,7 +1199,7 @@ pc_umeth_2D( if (bfbx.ok()) { pc_low_order_boundary( bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, - idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); } } if ( @@ -644,7 +1211,7 @@ pc_umeth_2D( if (bfbx.ok()) { pc_low_order_boundary( bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, - idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); } } } @@ -654,4 +1221,232 @@ pc_umeth_2D( pc_pdivu(i, j, k, pdivu, qec[0], qec[1], a[0], a[1], vol); }); } + +void +pc_umeth_eb_2D( + amrex::Box const& bx_to_fill, + const int* bclo, + const int* bchi, + const int* domlo, + const int* domhi, + amrex::Array4 const& q, + amrex::Array4 const& qaux, + amrex::Array4 const& srcQ, + const amrex::GpuArray, 2>& flx, + const amrex::GpuArray, 2>& qec, + const amrex::GpuArray, 2>& ap, + amrex::Array4 const& flag_arr, + const amrex::GpuArray del, + const amrex::Real dt, + const int ppm_type, + const bool use_flattening, + const int plm_iorder) +{ + BL_PROFILE("Godunov_umeth_2D_eb()"); + + amrex::Real const dx = del[0]; + amrex::Real const dy = del[1]; + amrex::Real const hdt = 0.5 * dt; + amrex::Real const hdtdy = 0.5 * dt / dy; + amrex::Real const hdtdx = 0.5 * dt / dx; + + const int bclx = bclo[0]; + const int bcly = bclo[1]; + const int bchx = bchi[0]; + const int bchy = bchi[1]; + const int dlx = domlo[0]; + const int dly = domlo[1]; + const int dhx = domhi[0]; + const int dhy = domhi[1]; + + const amrex::Box& bxg1 = grow(bx_to_fill, 1); + const amrex::Box& bxg2 = grow(bx_to_fill, 2); + + // X data + int cdir = 0; + const amrex::Box& xmbx = growHi(bxg2, cdir, 1); + const amrex::Box& xflxbx = surroundingNodes(bxg1, cdir); + amrex::FArrayBox qxm(xmbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qxp(bxg2, QVAR, amrex::The_Async_Arena()); + auto const& qxmarr = qxm.array(); + auto const& qxparr = qxp.array(); + + // Y data + cdir = 1; + const amrex::Box& ymbx = growHi(bxg2, cdir, 1); + const amrex::Box& yflxbx = surroundingNodes(bxg1, cdir); + amrex::FArrayBox qym(ymbx, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qyp(bxg2, QVAR, amrex::The_Async_Arena()); + auto const& qymarr = qym.array(); + auto const& qyparr = qyp.array(); + + if (ppm_type == 0) { + amrex::ParallelFor( + bxg2, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (flag_arr(i, j, k).isCovered()) { + return; + } + + amrex::Real slope[QVAR]; + + amrex::Real flat = 1.0; + // Calculate flattening in-place + if (use_flattening) { + for (int dir_flat = 0; dir_flat < AMREX_SPACEDIM; dir_flat++) { + flat = amrex::min( + flat, flatten_eb(i, j, k, dir_flat, flag_arr, q)); + } + } + + // + // X slopes and interp + // + for (int n = 0; n < QVAR; ++n) { + slope[n] = plm_slope_eb( + AMREX_D_DECL(i, j, k), n, 0, flag_arr, q, flat, plm_iorder); + } + pc_plm_d_eb( + i, j, k, 0, qxmarr, qxparr, slope, q, qaux(i, j, k, QC), dx, dt, + ap[0]); + + // + // Y slopes and interp + // + for (int n = 0; n < QVAR; n++) { + slope[n] = plm_slope_eb( + AMREX_D_DECL(i, j, k), n, 1, flag_arr, q, flat, plm_iorder); + } + pc_plm_d_eb( + i, j, k, 1, qymarr, qyparr, slope, q, qaux(i, j, k, QC), dy, dt, + ap[1]); + }); + } else { + amrex::Error("ppm_type must be 0 (PLM) when using EB"); + } + + // ******************************************************************************* + // These are the first flux estimates as per the corner-transport-upwind + // method X initial fluxes + // ******************************************************************************* + cdir = 0; + amrex::FArrayBox fx(xflxbx, NVAR, amrex::The_Async_Arena()); + auto const& fxarr = fx.array(); + amrex::FArrayBox qgdx(bxg2, NGDNV, amrex::The_Async_Arena()); + auto const& gdtemp = qgdx.array(); + amrex::ParallelFor( + xflxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclx, bchx, dlx, dhx, qxmarr, qxparr, fxarr, gdtemp, qaux, + cdir); + } + }); + + // ******************************************************************************* + // These are the first flux estimates as per the corner-transport-upwind + // method Y initial fluxes + // ******************************************************************************* + cdir = 1; + amrex::FArrayBox fy(yflxbx, NVAR, amrex::The_Async_Arena()); + auto const& fyarr = fy.array(); + amrex::ParallelFor( + yflxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bcly, bchy, dly, dhy, qymarr, qyparr, fyarr, qec[cdir], qaux, + cdir); + } + }); + + // ******************************************************************************* + // X interface corrections + // ******************************************************************************* + cdir = 0; + const amrex::Box& tybx = bxg1; + amrex::FArrayBox qm(bxg2, QVAR, amrex::The_Async_Arena()); + amrex::FArrayBox qp(bxg1, QVAR, amrex::The_Async_Arena()); + auto const& qmarr = qm.array(); + auto const& qparr = qp.array(); + + amrex::ParallelFor(tybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (!flag_arr(i, j, k).isCovered()) { + pc_transd( + AMREX_D_DECL(i, j, k), cdir, qmarr, qparr, qxmarr, qxparr, fyarr, srcQ, + qaux, qec[1], hdt, hdtdy, ap[cdir], ap[1]); + } + }); + + // This box must be grown by one in transverse direction so we can + // do tangential interpolation when taking divergence later + const amrex::Box& xfxbx = + surroundingNodes(grow(bx_to_fill, 1, 1 - cdir), cdir); + + // Final Riemann problem X + amrex::ParallelFor(xfxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bclx, bchx, dlx, dhx, qmarr, qparr, flx[cdir], qec[cdir], qaux, + cdir); + } + }); + + // ******************************************************************************* + // Y interface corrections + // ******************************************************************************* + cdir = 1; + const amrex::Box& txbx = bxg1; + + amrex::ParallelFor(txbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (!flag_arr(i, j, k).isCovered()) { + pc_transd( + AMREX_D_DECL(i, j, k), cdir, qmarr, qparr, qymarr, qyparr, fxarr, srcQ, + qaux, gdtemp, hdt, hdtdx, ap[cdir], ap[0]); + } + }); + + // This box must be grown by one in transverse direction so we can + // do tangential interpolation when taking divergence later + const amrex::Box& yfxbx = + surroundingNodes(grow(bx_to_fill, 1, 1 - cdir), cdir); + + // Final Riemann problem Y + amrex::ParallelFor(yfxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (ap[cdir](i, j, k) > 0.) { + pc_cmpflx( + i, j, k, bcly, bchy, dly, dhy, qmarr, qparr, flx[cdir], qec[cdir], qaux, + cdir); + } + }); + + // Fix bcnormal boundaries - always use PLM and don't do N+1/2 predictor + // because the user specifies conditions at N + // bcnormal used for "Hard" (inflow) and "UserBC" (user_bc) + for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { + if ( + bclo[idir] == PCPhysBCType::inflow || + bclo[idir] == PCPhysBCType::user_bc) { + // Box for fluxes at this boundary + amrex::Box bfbx = surroundingNodes(bx_to_fill, idir); + bfbx.setBig(idir, domlo[idir]); + if (bfbx.ok()) { + pc_low_order_boundary( + bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + } + } + if ( + bchi[idir] == PCPhysBCType::inflow || + bchi[idir] == PCPhysBCType::user_bc) { + // Box for fluxes at this boundary + amrex::Box bfbx = surroundingNodes(bx_to_fill, idir); + bfbx.setSmall(idir, domhi[idir] + 1); + if (bfbx.ok()) { + pc_low_order_boundary( + bfbx, bclo[idir], bchi[idir], domlo[idir], domhi[idir], plm_iorder, + use_flattening, idir, del[idir], dt, q, qaux, flx[idir], qec[idir]); + } + } + } +} + #endif diff --git a/Source/Hydro.H b/Source/Hydro.H index 7ce4ea196..c84d951b7 100644 --- a/Source/Hydro.H +++ b/Source/Hydro.H @@ -6,6 +6,7 @@ #include "PelePhysics.H" #include "Utilities.H" #include "Godunov.H" +#include "AMReX_EB_Redistribution.H" AMREX_GPU_DEVICE AMREX_FORCE_INLINE @@ -63,6 +64,407 @@ pc_srctoprim( #endif } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +pc_compute_hyp_wallflux( + int i, + int j, + int k, + amrex::Array4 const& q, + amrex::Array4 const& qa, + const amrex::Real axm, + const amrex::Real axp, + const amrex::Real aym, + const amrex::Real ayp, +#if (AMREX_SPACEDIM == 3) + const amrex::Real azm, + const amrex::Real azp, +#endif + amrex::GpuArray& fluxw) noexcept +{ + + // Here we have 2 options. Option 1: call the Riemann solver with dummy data. + // Option 2 is set the known flux at the wall, as we do in MOL. We are going + // with option 2. + // clang-format off + // Option 1 + // AMREX_D_TERM(amrex::Real u = q(i, j, k, QU);, amrex::Real v = q(i, j, k, QV); + // , amrex::Real w = q(i, j, k, QW);); + // const amrex::Real apnorm = std::sqrt(AMREX_D_TERM( + // (axm - axp) * (axm - axp), +(aym - ayp) * (aym - ayp), + // +(azm - azp) * (azm - azp))); + // const amrex::Real apnorminv = 1. / apnorm; + // amrex::Real un = AMREX_D_TERM( + // u * (axm - axp) * apnorminv, +v * (aym - ayp) * apnorminv, + // +w * (azm - azp) * apnorminv); + + // amrex::Real momfluxn, momfluxt1, momfluxt2; + // amrex::Real dummy_gu, dummy_gv, dummy_gv2, dummy_gp, dummy_gg; + // amrex::Real rho = q(i, j, k, QRHO); + // amrex::Real p = q(i, j, k, QPRES); + // amrex::Real cav = qa(i, j, k, QC); + // amrex::Real ustar; + // amrex::Real sp[NUM_SPECIES] = {0.0}; + // for (int n = 0; n < NUM_SPECIES; ++n) { + // sp[n] = q(i, j, k, QFS + n); + // } + + // const int bc_test_val = 1; + // riemann( + // rho, un, 0., 0., p, sp, rho, un, 0., 0., p, sp, bc_test_val, cav, ustar, + // fluxw[URHO], &fluxw[UFS], momfluxn, momfluxt1, momfluxt2, fluxw[UEDEN], + // fluxw[UEINT], dummy_gu, dummy_gv, dummy_gv2, dummy_gp, dummy_gg); + + // // enforce that no spurious flux was created by the Riemann solver + // for (int n = 0; n < NVAR; n++) { + // fluxw[n] = 0.; + // } + // AMREX_D_TERM(fluxw[UMX] = (axm - axp) * momfluxn; + // , fluxw[UMY] = (aym - ayp) * momfluxn; + // , fluxw[UMZ] = (azm - azp) * momfluxn;); + // clang-format on + + // Option 2 + amrex::ignore_unused(qa); + amrex::Real p = q(i, j, k, QPRES); + AMREX_D_TERM(fluxw[UMX] = (axm - axp) * p;, fluxw[UMY] = (aym - ayp) * p; + , fluxw[UMZ] = (azm - azp) * p;); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +eb_compute_div( + int i, + int j, + int k, + int n, + amrex::IntVect const& blo, + amrex::IntVect const& bhi, + amrex::Array4 const& q, + amrex::Array4 const& qaux, + amrex::Array4 const& divu, + AMREX_D_DECL( + amrex::Array4 const& u, + amrex::Array4 const& v, + amrex::Array4 const& w), + AMREX_D_DECL( + amrex::Array4 const& fx, + amrex::Array4 const& fy, + amrex::Array4 const& fz), + amrex::Array4 const& flag, + amrex::Array4 const& vfrc, + amrex::Array4 const& redistwgt, + AMREX_D_DECL( + amrex::Array4 const& apx, + amrex::Array4 const& apy, + amrex::Array4 const& apz), + AMREX_D_DECL( + amrex::Array4 const& fcx, + amrex::Array4 const& fcy, + amrex::Array4 const& fcz), + amrex::GpuArray const& dxinv, + const int eb_weights_type) +{ + AMREX_D_TERM(bool x_high = (i == bhi[0]);, bool y_high = (j == bhi[1]); + , bool z_high = (k == bhi[2])); + bool valid_cell = AMREX_D_TERM( + (blo[0] <= i) && (i <= bhi[0]), &&(blo[1] <= j) && (j <= bhi[1]), + &&(blo[2] <= k) && (k <= bhi[2])); + +#if (AMREX_SPACEDIM == 2) + amrex::Real cell_vol_inv = dxinv[0] * dxinv[1]; + if (flag(i, j, k).isCovered()) { + divu(i, j, k, n) = 0.0; + if (valid_cell) { + fx(i, j, k, n) = 0.; + fy(i, j, k, n) = 0.; + if (x_high) { + fx(i + 1, j, k, n) = 0.; + } + if (y_high) { + fy(i, j + 1, k, n) = 0.; + } + } + } else if (flag(i, j, k).isRegular()) { + divu(i, j, k, n) = cell_vol_inv * (u(i + 1, j, k, n) - u(i, j, k, n)) + + cell_vol_inv * (v(i, j + 1, k, n) - v(i, j, k, n)); + if (valid_cell) { + fx(i, j, k, n) = u(i, j, k, n); + fy(i, j, k, n) = v(i, j, k, n); + if (x_high) { + fx(i + 1, j, k, n) = u(i + 1, j, k, n); + } + if (y_high) { + fy(i, j + 1, k, n) = v(i, j + 1, k, n); + } + } + } else { + amrex::Real fxm = u(i, j, k, n); + if (apx(i, j, k) != 0.0 && apx(i, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i, j, k, 0))); + amrex::Real fracy = flag(i, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i, j, k, 0)) + : 0.0; + fxm = (1.0 - fracy) * fxm + fracy * u(i, jj, k, n); + } + if (valid_cell) { + fx(i, j, k, n) = fxm; + } + + amrex::Real fxp = u(i + 1, j, k, n); + if (apx(i + 1, j, k) != 0.0 && apx(i + 1, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i + 1, j, k, 0))); + amrex::Real fracy = flag(i + 1, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i + 1, j, k, 0)) + : 0.0; + fxp = (1.0 - fracy) * fxp + fracy * u(i + 1, jj, k, n); + } + if (valid_cell && x_high) { + fx(i + 1, j, k, n) = fxp; + } + + amrex::Real fym = v(i, j, k, n); + if (apy(i, j, k) != 0.0 && apy(i, j, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j, k, 0))); + amrex::Real fracx = flag(i, j, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j, k, 0)) + : 0.0; + fym = (1.0 - fracx) * fym + fracx * v(ii, j, k, n); + } + if (valid_cell) { + fy(i, j, k, n) = fym; + } + + amrex::Real fyp = v(i, j + 1, k, n); + if (apy(i, j + 1, k) != 0.0 && apy(i, j + 1, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j + 1, k, 0))); + amrex::Real fracx = flag(i, j + 1, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j + 1, k, 0)) + : 0.0; + fyp = (1.0 - fracx) * fyp + fracx * v(ii, j + 1, k, n); + } + if (valid_cell && y_high) { + fy(i, j + 1, k, n) = fyp; + } + + divu(i, j, k, n) = (1.0 / vfrc(i, j, k)) * cell_vol_inv * + (apx(i + 1, j, k) * fxp - apx(i, j, k) * fxm + + apy(i, j + 1, k) * fyp - apy(i, j, k) * fym); + + amrex::GpuArray flux_hyp_wall = {0.0}; + pc_compute_hyp_wallflux( + i, j, k, q, qaux, apx(i, j, k), apx(i + 1, j, k), apy(i, j, k), + apy(i, j + 1, k), flux_hyp_wall); + + // Here we assume dx == dy == dz + divu(i, j, k, n) += flux_hyp_wall[n] * dxinv[0] / vfrc(i, j, k); + } + +#else // 3-d starts here + + amrex::Real cell_vol_inv = dxinv[0] * dxinv[1] * dxinv[2]; + if (flag(i, j, k).isCovered()) { + divu(i, j, k, n) = 0.0; + if (valid_cell) { + fx(i, j, k, n) = 0.; + fy(i, j, k, n) = 0.; + fz(i, j, k, n) = 0.; + if (x_high) { + fx(i + 1, j, k, n) = 0.; + } + if (y_high) { + fy(i, j + 1, k, n) = 0.; + } + if (z_high) { + fz(i, j, k + 1, n) = 0.; + } + } + } else if (flag(i, j, k).isRegular()) { + divu(i, j, k, n) = cell_vol_inv * (u(i + 1, j, k, n) - u(i, j, k, n)) + + cell_vol_inv * (v(i, j + 1, k, n) - v(i, j, k, n)) + + cell_vol_inv * (w(i, j, k + 1, n) - w(i, j, k, n)); + if (valid_cell) { + fx(i, j, k, n) = u(i, j, k, n); + fy(i, j, k, n) = v(i, j, k, n); + fz(i, j, k, n) = w(i, j, k, n); + if (x_high) { + fx(i + 1, j, k, n) = u(i + 1, j, k, n); + } + if (y_high) { + fy(i, j + 1, k, n) = v(i, j + 1, k, n); + } + if (z_high) { + fz(i, j, k + 1, n) = w(i, j, k + 1, n); + } + } + } else { + amrex::Real fxm = u(i, j, k, n); + if (apx(i, j, k) != 0.0 && apx(i, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i, j, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcx(i, j, k, 1))); + amrex::Real fracy = flag(i, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i, j, k, 0)) + : 0.0; + amrex::Real fracz = flag(i, j, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcx(i, j, k, 1)) + : 0.0; + fxm = (1.0 - fracy) * (1.0 - fracz) * fxm + + fracy * (1.0 - fracz) * u(i, jj, k, n) + + fracz * (1.0 - fracy) * u(i, j, kk, n) + + fracy * fracz * u(i, jj, kk, n); + } + if (valid_cell) { + fx(i, j, k, n) = fxm; + } + + amrex::Real fxp = u(i + 1, j, k, n); + if (apx(i + 1, j, k) != 0.0 && apx(i + 1, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i + 1, j, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcx(i + 1, j, k, 1))); + amrex::Real fracy = flag(i + 1, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i + 1, j, k, 0)) + : 0.0; + amrex::Real fracz = flag(i + 1, j, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcx(i + 1, j, k, 1)) + : 0.0; + fxp = (1.0 - fracy) * (1.0 - fracz) * fxp + + fracy * (1.0 - fracz) * u(i + 1, jj, k, n) + + fracz * (1.0 - fracy) * u(i + 1, j, kk, n) + + fracy * fracz * u(i + 1, jj, kk, n); + } + if (valid_cell && x_high) { + fx(i + 1, j, k, n) = fxp; + } + + amrex::Real fym = v(i, j, k, n); + if (apy(i, j, k) != 0.0 && apy(i, j, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcy(i, j, k, 1))); + amrex::Real fracx = flag(i, j, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j, k, 0)) + : 0.0; + amrex::Real fracz = flag(i, j, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcy(i, j, k, 1)) + : 0.0; + fym = (1.0 - fracx) * (1.0 - fracz) * fym + + fracx * (1.0 - fracz) * v(ii, j, k, n) + + fracz * (1.0 - fracx) * v(i, j, kk, n) + + fracx * fracz * v(ii, j, kk, n); + } + if (valid_cell) { + fy(i, j, k, n) = fym; + } + + amrex::Real fyp = v(i, j + 1, k, n); + if (apy(i, j + 1, k) != 0.0 && apy(i, j + 1, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j + 1, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcy(i, j + 1, k, 1))); + amrex::Real fracx = flag(i, j + 1, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j + 1, k, 0)) + : 0.0; + amrex::Real fracz = flag(i, j + 1, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcy(i, j + 1, k, 1)) + : 0.0; + fyp = (1.0 - fracx) * (1.0 - fracz) * fyp + + fracx * (1.0 - fracz) * v(ii, j + 1, k, n) + + fracz * (1.0 - fracx) * v(i, j + 1, kk, n) + + fracx * fracz * v(ii, j + 1, kk, n); + } + if (valid_cell && y_high) { + fy(i, j + 1, k, n) = fyp; + } + + amrex::Real fzm = w(i, j, k, n); + if (apz(i, j, k) != 0.0 && apz(i, j, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k, 0))); + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k, 1))); + amrex::Real fracx = flag(i, j, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcz(i, j, k, 0)) + : 0.0; + amrex::Real fracy = flag(i, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcz(i, j, k, 1)) + : 0.0; + + fzm = (1.0 - fracx) * (1.0 - fracy) * fzm + + fracx * (1.0 - fracy) * w(ii, j, k, n) + + fracy * (1.0 - fracx) * w(i, jj, k, n) + + fracx * fracy * w(ii, jj, k, n); + } + if (valid_cell) { + fz(i, j, k, n) = fzm; + } + + amrex::Real fzp = w(i, j, k + 1, n); + if (apz(i, j, k + 1) != 0.0 && apz(i, j, k + 1) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k + 1, 0))); + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k + 1, 1))); + amrex::Real fracx = flag(i, j, k + 1).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcz(i, j, k + 1, 0)) + : 0.0; + amrex::Real fracy = flag(i, j, k + 1).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcz(i, j, k + 1, 1)) + : 0.0; + fzp = (1.0 - fracx) * (1.0 - fracy) * fzp + + fracx * (1.0 - fracy) * w(ii, j, k + 1, n) + + fracy * (1.0 - fracx) * w(i, jj, k + 1, n) + + fracx * fracy * w(ii, jj, k + 1, n); + } + if (valid_cell && z_high) { + fz(i, j, k + 1, n) = fzp; + } + + amrex::GpuArray flux_hyp_wall = {0.0}; + pc_compute_hyp_wallflux( + i, j, k, q, qaux, apx(i, j, k), apx(i + 1, j, k), apy(i, j, k), + apy(i, j + 1, k), apz(i, j, k), apz(i, j, k + 1), flux_hyp_wall); + + // With EB we assume dx == dy == dz + // NOTE: we have already made the fluxes extensive so we define by dx*dy*dz + // here ... + divu(i, j, k, n) = + cell_vol_inv / vfrc(i, j, k) * + (apx(i + 1, j, k) * fxp - apx(i, j, k) * fxm + apy(i, j + 1, k) * fyp - + apy(i, j, k) * fym + apz(i, j, k + 1) * fzp - apz(i, j, k) * fzm); + divu(i, j, k, n) += flux_hyp_wall[n] * dxinv[0] / vfrc(i, j, k); + } +#endif + + // The operations following this assume we have returned the negative of the + // divergence of fluxes. + divu(i, j, k, n) *= -1.0; + + // Go ahead and make the redistwgt array here since we'll need it in + // flux_redistribute + if (eb_weights_type == 0) { + redistwgt(i, j, k) = 1.0; + } else if (eb_weights_type == 1) { + redistwgt(i, j, k) = + q(i, j, k, QRHO) * + (q(i, j, k, QREINT) + 0.5 * (AMREX_D_TERM( + q(i, j, k, QU) * q(i, j, k, QU), + +q(i, j, k, QV) * q(i, j, k, QV), + +q(i, j, k, QW) * q(i, j, k, QW)))); + } else if (eb_weights_type == 2) { + redistwgt(i, j, k) = q(i, j, k, QRHO); + } else if (eb_weights_type == 3) { + redistwgt(i, j, k) = vfrc(i, j, k); + } +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -109,6 +511,319 @@ pc_divu( divu(i, j, k) = AMREX_D_TERM(ux, +vy, +wz); } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +eb_divu( + int i, + int j, + int k, + amrex::Array4 const& q, + amrex::Array4 const& divu, + amrex::Array4 const& vfrac, + amrex::GpuArray const& dxInv) +{ +#if (AMREX_SPACEDIM == 2) + const amrex::Real surrounding_vol = vfrac(i, j, k) + vfrac(i - 1, j, k) + + vfrac(i, j - 1, k) + + vfrac(i - 1, j - 1, k); + const amrex::Real cutoff = 4.0 * (1.0 - 1.e-12); +#elif (AMREX_SPACEDIM == 3) + const amrex::Real surrounding_vol = + vfrac(i, j, k) + vfrac(i - 1, j, k) + vfrac(i, j - 1, k) + + vfrac(i - 1, j - 1, k) + vfrac(i, j, k - 1) + vfrac(i - 1, j, k - 1) + + vfrac(i, j - 1, k - 1) + vfrac(i - 1, j - 1, k - 1); + const amrex::Real cutoff = 8.0 * (1.0 - 1.e-12); +#endif + + divu(i, j, k) = 0.0; + + if (surrounding_vol >= cutoff) { + +#if (AMREX_SPACEDIM == 2) + const amrex::Real ux = 0.5 * + (q(i, j, k, QU) - q(i - 1, j, k, QU) + + q(i, j - 1, k, QU) - q(i - 1, j - 1, k, QU)) * + dxInv[0]; + const amrex::Real vy = 0.5 * + (q(i, j, k, QV) - q(i, j - 1, k, QV) + + q(i - 1, j, k, QV) - q(i - 1, j - 1, k, QV)) * + dxInv[1]; +#elif (AMREX_SPACEDIM == 3) + const amrex::Real ux = + 0.25 * + (q(i, j, k, QU) - q(i - 1, j, k, QU) + q(i, j, k - 1, QU) - + q(i - 1, j, k - 1, QU) + q(i, j - 1, k, QU) - q(i - 1, j - 1, k, QU) + + q(i, j - 1, k - 1, QU) - q(i - 1, j - 1, k - 1, QU)) * + dxInv[0]; + + const amrex::Real vy = + 0.25 * + (q(i, j, k, QV) - q(i, j - 1, k, QV) + q(i, j, k - 1, QV) - + q(i, j - 1, k - 1, QV) + q(i - 1, j, k, QV) - q(i - 1, j - 1, k, QV) + + q(i - 1, j, k - 1, QV) - q(i - 1, j - 1, k - 1, QV)) * + dxInv[1]; + + const amrex::Real wz = + 0.25 * + (q(i, j, k, QW) - q(i, j, k - 1, QW) + q(i - 1, j, k, QW) - + q(i - 1, j, k - 1, QW) + q(i, j - 1, k, QW) - q(i, j - 1, k - 1, QW) + + q(i - 1, j - 1, k, QW) - q(i - 1, j - 1, k - 1, QW)) * + dxInv[2]; +#endif + divu(i, j, k) = AMREX_D_TERM(ux, +vy, +wz); + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +pc_eb_pdivu( + int i, + int j, + int k, + amrex::Array4 const& q, + AMREX_D_DECL( + amrex::Array4 const& q1, + amrex::Array4 const& q2, + amrex::Array4 const& q3), + amrex::Array4 const& divu, + amrex::Array4 const& flag, + amrex::Array4 const& vfrc, + AMREX_D_DECL( + amrex::Array4 const& apx, + amrex::Array4 const& apy, + amrex::Array4 const& apz), + AMREX_D_DECL( + amrex::Array4 const& fcx, + amrex::Array4 const& fcy, + amrex::Array4 const& fcz), + amrex::GpuArray const& dxinv) +{ + amrex::Real pdivu; + +#if (AMREX_SPACEDIM == 2) + if (flag(i, j, k).isRegular()) { + pdivu = 0.5 * ((q1(i + 1, j, k, GDPRES) + q1(i, j, k, GDPRES)) * + (q1(i + 1, j, k, GDU) - q1(i, j, k, GDU)) * dxinv[0] + + (q2(i, j + 1, k, GDPRES) + q2(i, j, k, GDPRES)) * + (q2(i, j + 1, k, GDV) - q2(i, j, k, GDV)) * dxinv[1]); + divu(i, j, k, UEINT) -= pdivu; + } else if (!flag(i, j, k).isCovered()) { + amrex::Real fxm, fxp, fym, fyp; + if (apx(i, j, k) > 0.0) { + fxm = q1(i, j, k, GDU); + if (apx(i, j, k) != 0.0 && apx(i, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i, j, k, 0))); + amrex::Real fracy = flag(i, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i, j, k, 0)) + : 0.0; + fxm = (1.0 - fracy) * fxm + fracy * q1(i, jj, k, GDU); + } + } else { + fxm = amrex::Real(0.0); + } + + if (apx(i + 1, j, k) > 0.0) { + fxp = q1(i + 1, j, k, GDU); + if (apx(i + 1, j, k) != 0.0 && apx(i + 1, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i + 1, j, k, 0))); + amrex::Real fracy = flag(i + 1, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i + 1, j, k, 0)) + : 0.0; + fxp = (1.0 - fracy) * fxp + fracy * q1(i + 1, jj, k, GDU); + } + } else { + fxp = amrex::Real(0.0); + } + + if (apy(i, j, k) > 0.0) { + fym = q2(i, j, k, GDV); + if (apy(i, j, k) != 0.0 && apy(i, j, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j, k, 0))); + amrex::Real fracx = flag(i, j, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j, k, 0)) + : 0.0; + fym = (1.0 - fracx) * fym + fracx * q2(ii, j, k, GDV); + } + } else { + fym = amrex::Real(0.0); + } + + if (apy(i, j + 1, k) > 0.0) { + fyp = q2(i, j + 1, k, GDV); + if (apy(i, j + 1, k) != 0.0 && apy(i, j + 1, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j + 1, k, 0))); + amrex::Real fracx = flag(i, j + 1, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j + 1, k, 0)) + : 0.0; + fyp = (1.0 - fracx) * fyp + fracx * q2(ii, j + 1, k, GDV); + } + } else { + fyp = amrex::Real(0.0); + } + + divu(i, j, k, UEINT) -= + (1.0 / vfrc(i, j, k)) * q(i, j, k, QPRES) * + (dxinv[0] * (apx(i + 1, j, k) * fxp - apx(i, j, k) * fxm) + + dxinv[1] * (apy(i, j + 1, k) * fyp - apy(i, j, k) * fym)); + } + +#else // 3-d starts here + + if (flag(i, j, k).isRegular()) { + pdivu = 0.5 * ((q1(i + 1, j, k, GDPRES) + q1(i, j, k, GDPRES)) * + (q1(i + 1, j, k, GDU) - q1(i, j, k, GDU)) * dxinv[0] + + (q2(i, j + 1, k, GDPRES) + q2(i, j, k, GDPRES)) * + (q2(i, j + 1, k, GDV) - q2(i, j, k, GDV)) * dxinv[1] + + (q3(i, j, k + 1, GDPRES) + q3(i, j, k, GDPRES)) * + (q3(i, j, k + 1, GDW) - q3(i, j, k, GDW)) * dxinv[2]); + + divu(i, j, k, UEINT) -= pdivu; + } else if (!flag(i, j, k).isCovered()) { + amrex::Real fxm, fxp, fym, fyp, fzm, fzp; + if (apx(i, j, k) > 0.0) { + fxm = q1(i, j, k, GDU); + if (apx(i, j, k) != 0.0 && apx(i, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i, j, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcx(i, j, k, 1))); + amrex::Real fracy = flag(i, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i, j, k, 0)) + : 0.0; + amrex::Real fracz = flag(i, j, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcx(i, j, k, 1)) + : 0.0; + fxm = (1.0 - fracy) * (1.0 - fracz) * fxm + + fracy * (1.0 - fracz) * q1(i, jj, k, GDU) + + fracz * (1.0 - fracy) * q1(i, j, kk, GDU) + + fracy * fracz * q1(i, jj, kk, GDU); + } + } else { + fxm = amrex::Real(0.0); + } + + if (apx(i + 1, j, k) > 0.0) { + fxp = q1(i + 1, j, k, GDU); + if (apx(i + 1, j, k) != 0.0 && apx(i + 1, j, k) != 1.0) { + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcx(i + 1, j, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcx(i + 1, j, k, 1))); + amrex::Real fracy = flag(i + 1, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcx(i + 1, j, k, 0)) + : 0.0; + amrex::Real fracz = flag(i + 1, j, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcx(i + 1, j, k, 1)) + : 0.0; + fxp = (1.0 - fracy) * (1.0 - fracz) * fxp + + fracy * (1.0 - fracz) * q1(i + 1, jj, k, GDU) + + fracz * (1.0 - fracy) * q1(i + 1, j, kk, GDU) + + fracy * fracz * q1(i + 1, jj, kk, GDU); + } + } else { + fxp = amrex::Real(0.0); + } + + if (apy(i, j, k) > 0.0) { + fym = q2(i, j, k, GDV); + if (apy(i, j, k) != 0.0 && apy(i, j, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcy(i, j, k, 1))); + amrex::Real fracx = flag(i, j, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j, k, 0)) + : 0.0; + amrex::Real fracz = flag(i, j, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcy(i, j, k, 1)) + : 0.0; + fym = (1.0 - fracx) * (1.0 - fracz) * fym + + fracx * (1.0 - fracz) * q2(ii, j, k, GDV) + + fracz * (1.0 - fracx) * q2(i, j, kk, GDV) + + fracx * fracz * q2(ii, j, kk, GDV); + } + } else { + fym = amrex::Real(0.0); + } + + if (apy(i, j + 1, k) > 0.0) { + fyp = q2(i, j + 1, k, GDV); + if (apy(i, j + 1, k) != 0.0 && apy(i, j + 1, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcy(i, j + 1, k, 0))); + int kk = + k + static_cast(amrex::Math::copysign(1.0, fcy(i, j + 1, k, 1))); + amrex::Real fracx = flag(i, j + 1, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcy(i, j + 1, k, 0)) + : 0.0; + amrex::Real fracz = flag(i, j + 1, k).isConnected(0, 0, kk - k) + ? amrex::Math::abs(fcy(i, j + 1, k, 1)) + : 0.0; + fyp = (1.0 - fracx) * (1.0 - fracz) * fyp + + fracx * (1.0 - fracz) * q2(ii, j + 1, k, GDV) + + fracz * (1.0 - fracx) * q2(i, j + 1, kk, GDV) + + fracx * fracz * q2(ii, j + 1, kk, GDV); + } + } else { + fyp = amrex::Real(0.0); + } + + if (apz(i, j, k) > 0.0) { + fzm = q3(i, j, k, GDW); + if (apz(i, j, k) != 0.0 && apz(i, j, k) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k, 0))); + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k, 1))); + amrex::Real fracx = flag(i, j, k).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcz(i, j, k, 0)) + : 0.0; + amrex::Real fracy = flag(i, j, k).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcz(i, j, k, 1)) + : 0.0; + + fzm = (1.0 - fracx) * (1.0 - fracy) * fzm + + fracx * (1.0 - fracy) * q3(ii, j, k, GDW) + + fracy * (1.0 - fracx) * q3(i, jj, k, GDW) + + fracx * fracy * q3(ii, jj, k, GDW); + } + } else { + fzm = amrex::Real(0.0); + } + + if (apz(i, j, k + 1) > 0.0) { + fzp = q3(i, j, k + 1, GDW); + if (apz(i, j, k + 1) != 0.0 && apz(i, j, k + 1) != 1.0) { + int ii = + i + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k + 1, 0))); + int jj = + j + static_cast(amrex::Math::copysign(1.0, fcz(i, j, k + 1, 1))); + amrex::Real fracx = flag(i, j, k + 1).isConnected(ii - i, 0, 0) + ? amrex::Math::abs(fcz(i, j, k + 1, 0)) + : 0.0; + amrex::Real fracy = flag(i, j, k + 1).isConnected(0, jj - j, 0) + ? amrex::Math::abs(fcz(i, j, k + 1, 1)) + : 0.0; + fzp = (1.0 - fracx) * (1.0 - fracy) * fzp + + fracx * (1.0 - fracy) * q3(ii, j, k + 1, GDW) + + fracy * (1.0 - fracx) * q3(i, jj, k + 1, GDW) + + fracx * fracy * q3(ii, jj, k + 1, GDW); + } + } else { + fzp = amrex::Real(0.0); + } + + divu(i, j, k, UEINT) -= + (1.0 / vfrc(i, j, k)) * q(i, j, k, QPRES) * + (dxinv[0] * (apx(i + 1, j, k) * fxp - apx(i, j, k) * fxm) + + dxinv[1] * (apy(i, j + 1, k) * fyp - apy(i, j, k) * fym) + + dxinv[2] * (apz(i, j, k + 1) * fzp - apz(i, j, k) * fzm)); + } +#endif +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -144,6 +859,24 @@ pc_ext_flx( } } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +pc_ext_flx_eb( + const int i, + const int j, + const int k, + amrex::Array4 const& flx, + amrex::Real area, + amrex::Array4 const& aper) +{ + // The fluxes are multiplied by area fractions inside amrex + // so here we just weight them by the regular areas + if (aper(i, j, k) > 0.) { + for (int n = 0; n < NVAR; ++n) { + flx(i, j, k, n) *= area; + } + } +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -228,6 +961,9 @@ void pc_umdrv_eb( amrex::BCRec const* bcs_d_ptr, const std::string& redistribution_type, const int eb_weights_type, + const int eb_srd_max_order, + const bool eb_clean_massfrac, + const amrex::Real eb_clean_massfrac_threshold, amrex::Real cflLoc); void pc_adjust_fluxes( @@ -238,12 +974,60 @@ void pc_adjust_fluxes( a, const amrex::Array4& divu, const amrex::GpuArray& del, + const int* domlo, + const int* domhi, + const int* bclo, + const int* bchi, const amrex::Real difmag); +void pc_adjust_fluxes_eb( + const amrex::Box& /*bx*/, + amrex::Array4 const& q_arr, + amrex::Array4 const& u_arr, + AMREX_D_DECL( + amrex::Array4 const& apx, + amrex::Array4 const& apy, + amrex::Array4 const& apz), + amrex::Array4 const& vfrac, + const amrex::GpuArray dx, + const amrex::GpuArray dxinv, + const amrex::GpuArray, AMREX_SPACEDIM> flux, + const amrex::Geometry& geom, + const int* bclo, + const int* bchi, + amrex::Real difmag); + void pc_consup( const amrex::Box& bx, const amrex::Array4& update, const amrex::GpuArray, AMREX_SPACEDIM>& flx, const amrex::Array4& vol, const amrex::Array4& pdivu); + +void pc_consup_eb( + const amrex::Box& bx, + amrex::Array4 const& q_arr, + amrex::Array4 const& qaux_arr, + amrex::Array4 const& divc_arr, + amrex::Array4 const& redistwgt_arr, + AMREX_D_DECL( + amrex::Array4 const& q1, + amrex::Array4 const& q2, + amrex::Array4 const& q3), + AMREX_D_DECL( + amrex::Array4 const& apx, + amrex::Array4 const& apy, + amrex::Array4 const& apz), + AMREX_D_DECL( + amrex::Array4 const& fcx, + amrex::Array4 const& fcy, + amrex::Array4 const& fcz), + amrex::Array4 const& vfrac, + amrex::Array4 const& flag, + const amrex::GpuArray dxinv, + const amrex::GpuArray, AMREX_SPACEDIM> + flux_tmp, + const amrex::GpuArray, AMREX_SPACEDIM> flux, + const int eb_weights_type); + #endif diff --git a/Source/Hydro.cpp b/Source/Hydro.cpp index 80cf331a7..fc61b591c 100644 --- a/Source/Hydro.cpp +++ b/Source/Hydro.cpp @@ -20,7 +20,7 @@ PeleC::construct_hydro_source( BL_PROFILE("PeleC::construct_hydro_source()"); if ((verbose != 0) && amrex::ParallelDescriptor::IOProcessor()) { - amrex::Print() << "... Computing hydro advance" << std::endl; + amrex::Print() << "... Computing Godunov hydro advance" << std::endl; } AMREX_ASSERT(S.nGrow() >= numGrow() + nGrowF); @@ -188,6 +188,8 @@ PeleC::construct_hydro_source( const amrex::Box& fbxg_i = grow(fbx, ngrow_bx); if (flag_fab.getType(fbxg_i) == amrex::FabType::singlevalued) { + BL_PROFILE("PeleC::umdrv_eb()") + auto const& vfrac_arr = vfrac.const_array(mfi); amrex::EBFluxRegister* fr_as_crse = @@ -238,7 +240,8 @@ PeleC::construct_hydro_source( p_rrflag_as_crse->array(), as_fine, dm_as_fine.array(), level_mask.const_array(mfi), dt, ppm_type, plm_iorder, use_flattening, difmag, bcs_d.data(), redistribution_type, - eb_weights_type, cflLoc); + eb_weights_type, eb_srd_max_order, eb_clean_massfrac, + eb_clean_massfrac_threshold, cflLoc); } else if (flag_fab.getType(fbxg_i) == amrex::FabType::regular) { BL_PROFILE("PeleC::umdrv()"); @@ -371,7 +374,8 @@ pc_umdrv( pc_divu(i, j, k, q, AMREX_D_DECL(dx[0], dx[1], dx[2]), divuarr); }); - pc_adjust_fluxes(bx, uin, flx, a, divuarr, dx, difmag); + pc_adjust_fluxes( + bx, uin, flx, a, divuarr, dx, domlo, domhi, bclo, bchi, difmag); // consup pc_consup(bx, uout, flx, vol, pdivuarr); @@ -386,6 +390,10 @@ pc_adjust_fluxes( a, const amrex::Array4& divu, const amrex::GpuArray& del, + const int* domlo, + const int* domhi, + const int* bclo, + const int* bchi, const amrex::Real difmag) { BL_PROFILE("PeleC::pc_adjust_fluxes()"); @@ -393,8 +401,14 @@ pc_adjust_fluxes( for (int dir = 0; dir < AMREX_SPACEDIM; dir++) { amrex::Box const& fbx = surroundingNodes(bx, dir); const amrex::Real dx = del[dir]; + const int domlo_dir = domlo[dir]; + const int domhi_dir = domhi[dir]; + const int bclo_dir = bclo[dir]; + const int bchi_dir = bchi[dir]; amrex::ParallelFor(fbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - pc_artif_visc(AMREX_D_DECL(i, j, k), flx[dir], divu, u, dx, difmag, dir); + pc_artif_visc( + AMREX_D_DECL(i, j, k), flx[dir], divu, u, dx, difmag, dir, domlo_dir, + domhi_dir, bclo_dir, bchi_dir); // Normalize Species Flux pc_norm_spec_flx(i, j, k, flx[dir]); // Make flux extensive @@ -403,6 +417,136 @@ pc_adjust_fluxes( } } +void +pc_adjust_fluxes_eb( + const amrex::Box& /*bx*/, + amrex::Array4 const& q_arr, + amrex::Array4 const& u_arr, + AMREX_D_DECL( + amrex::Array4 const& apx, + amrex::Array4 const& apy, + amrex::Array4 const& apz), + amrex::Array4 const& vfrac, + const amrex::GpuArray dx, + const amrex::GpuArray dxinv, + const amrex::GpuArray, AMREX_SPACEDIM> flux, + const amrex::Geometry& geom, + const int* bclo, + const int* bchi, + amrex::Real difmag) +{ + BL_PROFILE("PeleC::pc_adjust_fluxes_eb()"); + amrex::Real areafac; + + // These are the fluxes on face centroids -- they are defined in + // eb_compute_div + // and are the fluxes that go into the flux registers + AMREX_D_TERM(auto const& fx_arr = flux[0];, auto const& fy_arr = flux[1]; + , auto const& fz_arr = flux[2];); + + // Compute divu to be used in artificial viscosity + auto bx_divu = amrex::Box(q_arr); + bx_divu.grow(-1); + bx_divu.surroundingNodes(); + amrex::FArrayBox divu(bx_divu, NVAR, amrex::The_Async_Arena()); + divu.setVal(0.0); + + amrex::Box nddom = + amrex::convert(geom.growPeriodicDomain(16), amrex::IntVect(1)); + bx_divu &= nddom; + + auto const& divu_arr = divu.array(); + amrex::ParallelFor( + bx_divu, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + eb_divu(i, j, k, q_arr, divu_arr, vfrac, dxinv); + }); + + // Flux alterations + amrex::Box const& fbx = amrex::Box(fx_arr); +#if (AMREX_SPACEDIM == 2) + areafac = dx[1]; +#elif (AMREX_SPACEDIM == 3) + areafac = dx[1] * dx[2]; +#endif + + auto const* domlo = geom.Domain().loVect(); + auto const* domhi = geom.Domain().hiVect(); + + int domlo_dir = domlo[0]; + int domhi_dir = domhi[0]; + int bclo_dir = bclo[0]; + int bchi_dir = bchi[0]; + + amrex::ParallelFor(fbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + pc_artif_visc( + AMREX_D_DECL(i, j, k), fx_arr, divu_arr, u_arr, dx[0], difmag, 0, + domlo_dir, domhi_dir, bclo_dir, bchi_dir); + + // Normalize Species Flux + if (apx(i, j, k) > 0.) { + pc_norm_spec_flx(i, j, k, fx_arr); + } + + // Make flux extensive + pc_ext_flx_eb(i, j, k, fx_arr, areafac, apx); + }); + + // Flux alterations + amrex::Box const& fby = amrex::Box(fy_arr); +#if (AMREX_SPACEDIM == 2) + areafac = dx[0]; +#elif (AMREX_SPACEDIM == 3) + areafac = dx[0] * dx[2]; +#endif + + domlo_dir = domlo[1]; + domhi_dir = domhi[1]; + bclo_dir = bclo[1]; + bchi_dir = bchi[1]; + + amrex::ParallelFor(fby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + pc_artif_visc( + AMREX_D_DECL(i, j, k), fy_arr, divu_arr, u_arr, dx[1], difmag, 1, + domlo_dir, domhi_dir, bclo_dir, bchi_dir); + + // Normalize Species Flux + if (apy(i, j, k) > 0.) { + pc_norm_spec_flx(i, j, k, fy_arr); + } + + // Make flux extensive + pc_ext_flx_eb(i, j, k, fy_arr, areafac, apy); + }); + +#if (AMREX_SPACEDIM == 3) + // Flux alterations + amrex::Box const& fbz = amrex::Box(fz_arr); + areafac = dx[1] * dx[0]; + + domlo_dir = domlo[2]; + domhi_dir = domhi[2]; + bclo_dir = bclo[2]; + bchi_dir = bchi[2]; + + amrex::ParallelFor(fbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + pc_artif_visc( + i, j, k, fz_arr, divu_arr, u_arr, dx[2], difmag, 2, domlo_dir, domhi_dir, + bclo_dir, bchi_dir); + + // Normalize Species Flux + if (apz(i, j, k) > 0.) { + pc_norm_spec_flx(i, j, k, fz_arr); + } + + // Make flux extensive + pc_ext_flx_eb(i, j, k, fz_arr, areafac, apz); + }); + +#endif + + amrex::Gpu::streamSynchronize(); +} + void pc_consup( const amrex::Box& bx, @@ -418,41 +562,235 @@ pc_consup( }); } +void +pc_consup_eb( + const amrex::Box& bx, + amrex::Array4 const& q_arr, + amrex::Array4 const& qaux_arr, + amrex::Array4 const& divc_arr, + amrex::Array4 const& redistwgt_arr, + AMREX_D_DECL( + amrex::Array4 const& q1, + amrex::Array4 const& q2, + amrex::Array4 const& q3), + AMREX_D_DECL( + amrex::Array4 const& apx, + amrex::Array4 const& apy, + amrex::Array4 const& apz), + AMREX_D_DECL( + amrex::Array4 const& fcx, + amrex::Array4 const& fcy, + amrex::Array4 const& fcz), + amrex::Array4 const& vfrac, + amrex::Array4 const& flag, + const amrex::GpuArray dxinv, + const amrex::GpuArray, AMREX_SPACEDIM> + flux_tmp, + const amrex::GpuArray, AMREX_SPACEDIM> flux, + + const int eb_weights_type) +{ + const amrex::Box& bxg_i = amrex::Box(divc_arr); + + // These are the fluxes we computed in MOL_umeth_eb and modified in + // adjust_fluxes_eb + // -- they live at face centers + AMREX_D_TERM(auto const& fx_in_arr = flux_tmp[0]; + , auto const& fy_in_arr = flux_tmp[1]; + , auto const& fz_in_arr = flux_tmp[2];); + + // These are the fluxes on face centroids -- they are defined in + // eb_compute_div + // and are the fluxes that go into the flux registers + AMREX_D_TERM(auto const& fx_out_arr = flux[0]; + , auto const& fy_out_arr = flux[1]; + , auto const& fz_out_arr = flux[2];); + + auto const& blo = bx.smallEnd(); + auto const& bhi = bx.bigEnd(); + + amrex::ParallelFor( + bxg_i, NVAR, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { + // This does the divergence but not the redistribution -- we will do that + // later We do compute the weights here though + if (!flag(i, j, k).isCovered()) { + eb_compute_div( + i, j, k, n, blo, bhi, q_arr, qaux_arr, divc_arr, + AMREX_D_DECL(fx_in_arr, fy_in_arr, fz_in_arr), + AMREX_D_DECL(fx_out_arr, fy_out_arr, fz_out_arr), flag, vfrac, + redistwgt_arr, AMREX_D_DECL(apx, apy, apz), + AMREX_D_DECL(fcx, fcy, fcz), dxinv, eb_weights_type); + } + }); + + // Update UEINT update with pdivu term + amrex::ParallelFor(bxg_i, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (!flag(i, j, k).isCovered()) { + pc_eb_pdivu( + i, j, k, q_arr, AMREX_D_DECL(q1, q2, q3), divc_arr, flag, vfrac, + AMREX_D_DECL(apx, apy, apz), AMREX_D_DECL(fcx, fcy, fcz), dxinv); + } + }); +} + void pc_umdrv_eb( - const amrex::Box& /*bx*/, - const amrex::Box& /*bxg_i*/, - const amrex::MFIter& /*mfi*/, - const amrex::Geometry& /*geom*/, - const amrex::EBFArrayBoxFactory* /*fact*/, - const int* /*bclo*/, - const int* /*bchi*/, - const amrex::Array4& /*uin*/, - const amrex::Array4& /*uout*/, - const amrex::Array4& /*q*/, - const amrex::Array4& /*qaux*/, - const amrex::Array4& /*src_q*/, - const amrex::Array4& /*vf*/, - const amrex::Array4& /*flag*/, - const amrex::GpuArray& /*dx*/, - const amrex::GpuArray& /*dxInv*/, - const amrex:: - GpuArray, AMREX_SPACEDIM>& /*flx*/, - const int /*as_crse*/, - const amrex::Array4& /*drho_as_crse*/, - const amrex::Array4& /*rrflag_as_crse*/, - const int /*as_fine*/, - const amrex::Array4& /*dm_as_fine*/, - const amrex::Array4& /*lev_mask*/, - const amrex::Real /*dt*/, - const int /*ppm_type*/, - const int /*plm_iorder*/, - const bool /*use_flattening*/, - const amrex::Real /*difmag*/, - amrex::BCRec const* /*bcs_d_ptr*/, - const std::string& /*redistribution_type*/, - const int /*eb_weights_type*/, + const amrex::Box& bx, + const amrex::Box& bxg_i, + const amrex::MFIter& mfi, + const amrex::Geometry& geom, + const amrex::EBFArrayBoxFactory* fact, + const int* bclo, + const int* bchi, + const amrex::Array4& uin, + const amrex::Array4& uout, + const amrex::Array4& q, + const amrex::Array4& qaux, + const amrex::Array4& src_q, + const amrex::Array4& vf, + const amrex::Array4& flag, + const amrex::GpuArray& dx, + const amrex::GpuArray& dxInv, + const amrex::GpuArray, AMREX_SPACEDIM>& flx, + const int as_crse, + const amrex::Array4& drho_as_crse, + const amrex::Array4& rrflag_as_crse, + const int as_fine, + const amrex::Array4& dm_as_fine, + const amrex::Array4& lev_mask, + const amrex::Real dt, + const int ppm_type, + const int plm_iorder, + const bool use_flattening, + const amrex::Real difmag, + amrex::BCRec const* bcs_d_ptr, + const std::string& redistribution_type, + const int eb_weights_type, + const int eb_srd_max_order, + const bool eb_clean_massfrac, + const amrex::Real eb_clean_massfrac_threshold, amrex::Real /*cflLoc*/) { - amrex::Abort("EB Godunov not implemented yet"); + BL_PROFILE("PeleC::pc_umdrv_eb()"); + + const amrex::Box& bxg_ii = grow(bxg_i, 1); + + amrex::Array4 AMREX_D_DECL(fcx, fcy, fcz), + AMREX_D_DECL(apx, apy, apz), ccc; + AMREX_D_TERM(fcx = fact->getFaceCent()[0]->const_array(mfi); + , fcy = fact->getFaceCent()[1]->const_array(mfi); + , fcz = fact->getFaceCent()[2]->const_array(mfi);); + AMREX_D_TERM(apx = fact->getAreaFrac()[0]->const_array(mfi); + , apy = fact->getAreaFrac()[1]->const_array(mfi); + , apz = fact->getAreaFrac()[2]->const_array(mfi);); + ccc = fact->getCentroid().const_array(mfi); + + const auto* domlo = geom.Domain().loVect(); + const auto* domhi = geom.Domain().hiVect(); + + // Temporary FArrayBoxes + amrex::FArrayBox divu(bxg_ii, 1, amrex::The_Async_Arena()); + auto const& divuarr = divu.array(); + + amrex::FArrayBox qec[AMREX_SPACEDIM]; + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) { + const amrex::Box eboxes = amrex::surroundingNodes(grow(bxg_ii, 1), dir); + qec[dir].resize(eboxes, NGDNV, amrex::The_Async_Arena()); + } + const amrex::GpuArray, AMREX_SPACEDIM> + qec_arr{{AMREX_D_DECL(qec[0].array(), qec[1].array(), qec[2].array())}}; + + // Quantities for redistribution + amrex::FArrayBox divc, redistwgt; + if (redistribution_type == "StateRedist") { + divc.resize(bxg_i, NVAR); // This will hold "dUdt" before redistribution + redistwgt.resize( + bxg_i, NVAR); // This will be "scratch" which holds "Uold + dt*dUdt" + } else { + divc.resize(bxg_i, NVAR); // This will hold "dUdt" before redistribution + redistwgt.resize( + bxg_i, 1); // This will hold the weights used in flux redistribution + } + divc.setVal(0.0); + redistwgt.setVal(0.0); + + // Because we are going to redistribute, we put the divergence into divc + // rather than directly into dsdt_arr + auto const& divc_arr = divc.array(); + auto const& redistwgt_arr = redistwgt.array(); + + // We need one extra in flux_tmp so we can tangentially interpolate fluxes + amrex::FArrayBox flux_tmp[AMREX_SPACEDIM]; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + flux_tmp[idim].resize( + amrex::surroundingNodes(bxg_ii, idim), NVAR, amrex::The_Async_Arena()); + flux_tmp[idim].setVal(0.); + } + + const amrex::GpuArray, AMREX_SPACEDIM> + flux_tmp_arr{{AMREX_D_DECL( + flux_tmp[0].array(), flux_tmp[1].array(), flux_tmp[2].array())}}; + const amrex::GpuArray, AMREX_SPACEDIM> + ap{{AMREX_D_DECL(apx, apy, apz)}}; + + // Define divc, the update before redistribution + // Also construct the redistribution weights for flux redistribution if + // necessary +#if AMREX_SPACEDIM == 1 + amrex::Abort("PLM isn't implemented in 1D."); +#elif AMREX_SPACEDIM == 2 + pc_umeth_eb_2D( + amrex::Box(divc_arr), bclo, bchi, domlo, domhi, q, qaux, src_q, + flux_tmp_arr, qec_arr, ap, flag, dx, dt, ppm_type, use_flattening, + plm_iorder); +#elif AMREX_SPACEDIM == 3 + pc_umeth_eb_3D( + amrex::Box(divc_arr), bclo, bchi, domlo, domhi, q, qaux, src_q, + flux_tmp_arr, qec_arr, ap, flag, dx, dt, ppm_type, use_flattening, + plm_iorder); +#endif + + // Construct divu + AMREX_D_TERM(const amrex::Real dx0 = dx[0];, const amrex::Real dx1 = dx[1]; + , const amrex::Real dx2 = dx[2];); + + amrex::ParallelFor( + bxg_ii, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (flag(i, j, k).isRegular()) { + pc_divu(i, j, k, q, AMREX_D_DECL(dx0, dx1, dx2), divuarr); + } else { + divuarr(i, j, k) = amrex::Real(0.0); + } + }); + + pc_adjust_fluxes_eb( + bx, q, uin, AMREX_D_DECL(apx, apy, apz), vf, dx, dxInv, flux_tmp_arr, geom, + bclo, bchi, difmag); + + pc_consup_eb( + bx, q, qaux, divc_arr, redistwgt_arr, + AMREX_D_DECL(qec_arr[0], qec_arr[1], qec_arr[2]), + AMREX_D_DECL(apx, apy, apz), AMREX_D_DECL(fcx, fcy, fcz), vf, flag, dxInv, + flux_tmp_arr, flx, eb_weights_type); + + const int l_ncomp = uout.nComp(); + const int level_mask_not_covered = constants::level_mask_notcovered(); + const bool use_wts_in_divnc = false; + + const amrex::Real fac_for_redist = 1.0; + { + BL_PROFILE("ApplyMLRedistribution()"); + ApplyMLRedistribution( + bx, l_ncomp, uout, divc_arr, uin, redistwgt_arr, flag, + AMREX_D_DECL(apx, apy, apz), vf, AMREX_D_DECL(fcx, fcy, fcz), ccc, + bcs_d_ptr, geom, dt, redistribution_type, as_crse, drho_as_crse, + rrflag_as_crse, as_fine, dm_as_fine, lev_mask, level_mask_not_covered, + fac_for_redist, use_wts_in_divnc, 0, eb_srd_max_order); + } + auto const& flags = fact->getMultiEBCellFlagFab(); + const auto& flag_fab = flags[mfi]; + amrex::FabType typ = flag_fab.getType(bx); + pc_post_eb_redistribution( + bx, dt, eb_clean_massfrac, eb_clean_massfrac_threshold, uin, typ, flag, + redistwgt_arr, uout); } diff --git a/Source/PLM.H b/Source/PLM.H index 8ca79b428..eabac5115 100644 --- a/Source/PLM.H +++ b/Source/PLM.H @@ -23,6 +23,7 @@ plm_slope( const int n, const int dir, amrex::Array4 const& q, + const amrex::Real flat, const int order) { if (order == 1) { @@ -72,8 +73,45 @@ plm_slope( dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - // Flattening could be done here (see Nyx if we want to do it) - return dsgn * amrex::min(dlim, std::abs(dtemp)); + return flat * dsgn * amrex::min(dlim, std::abs(dtemp)); +} + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +plm_slope_eb( + AMREX_D_DECL(const int i, const int j, const int k), + const int n, + const int dir, + amrex::Array4 const& flags, + amrex::Array4 const& q, + const amrex::Real flat, + const int order) +{ + const amrex::IntVect iv{AMREX_D_DECL(i, j, k)}; + const amrex::IntVect dvec(amrex::IntVect::TheDimensionVector(dir)); + const amrex::IntVect ivm2(iv - 2 * dvec); + const amrex::IntVect ivm(iv - dvec); + const amrex::IntVect ivp(iv + dvec); + const amrex::IntVect ivp2(iv + 2 * dvec); + + // We have enough cells in the x-direction to do 4th order slopes + // centered on (i,j,k) + amrex::Real slope = 0.0; + if ( + (order == 4) && !flags(iv).isCovered() && !flags(ivm).isCovered() && + !flags(ivm2).isCovered() && !flags(ivp).isCovered() && + !flags(ivp2).isCovered()) { + slope = plm_slope(AMREX_D_DECL(i, j, k), n, dir, q, flat, order); + } + // We have enough cells in the x-direction to do 2nd order slopes + // centered on (i,j,k) + else if ( + (order > 1) && !flags(iv).isCovered() && !flags(ivm).isCovered() && + !flags(ivp).isCovered()) { + slope = plm_slope(AMREX_D_DECL(i, j, k), n, dir, q, flat, 2); + } + return slope; } AMREX_GPU_DEVICE @@ -97,6 +135,32 @@ pc_plm_d_passive( qm(ivp, n) = q(iv, n) + acmpleft; } +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +pc_plm_d_eb_passive( + const amrex::IntVect iv, + const amrex::IntVect ivp, + const int n, + const amrex::Real spzerom, + const amrex::Real spzerop, + amrex::Real const slope[QVAR], + amrex::Array4 const& q, + amrex::Array4 const& qp, + amrex::Array4 const& qm, + amrex::Array4 const& area) +{ + if (area(iv) > 0.0) { + const amrex::Real acmprght = 0.5 * (-1.0 - spzerom) * slope[n]; + qp(iv, n) = q(iv, n) + acmprght; + } + + if (area(ivp) > 0.0) { + const amrex::Real acmpleft = 0.5 * (1.0 - spzerop) * slope[n]; + qm(ivp, n) = q(iv, n) + acmpleft; + } +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -248,4 +312,186 @@ pc_plm_d( #endif } +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +pc_plm_d_eb( + const int i, + const int j, + const int k, + const int dir, + amrex::Array4 const& qm, + amrex::Array4 const& qp, + amrex::Real const slope[QVAR], + amrex::Array4 const& q, + amrex::Real const cc, + amrex::Real const dx, + amrex::Real const dt, + amrex::Array4 const& area) +{ + amrex::ignore_unused(k); + const amrex::IntVect iv{AMREX_D_DECL(i, j, k)}; + const amrex::IntVect ivp(iv + amrex::IntVect::TheDimensionVector(dir)); + const amrex::GpuArray bdim{ + {static_cast(dir == 0), static_cast(dir == 1), + static_cast(dir == 2)}}; + const amrex::GpuArray l_idx{ + {bdim[0] * 0 + bdim[1] * 1 + bdim[2] * 2, + bdim[0] * 1 + bdim[1] * 0 + bdim[2] * 0, + bdim[0] * 2 + bdim[1] * 2 + bdim[2] * 1}}; + + const amrex::Real dtdx = dt / dx; + const amrex::Real cs = cc * cc; + const amrex::Real rho = q(iv, QRHO); + const amrex::GpuArray vel{ + AMREX_D_DECL(q(iv, QU), q(iv, QV), q(iv, QW))}; + const amrex::Real p = q(iv, QPRES); + const amrex::Real rhoe = q(iv, QREINT); + const amrex::Real enth = ((rhoe + p) / rho) / cs; + const amrex::Real drho = slope[QRHO]; + const amrex::GpuArray dvel{ + AMREX_D_DECL(slope[QU], slope[QV], slope[QW])}; + const amrex::Real dp = slope[QPRES]; + const amrex::Real drhoe = slope[QREINT]; + const amrex::Real alpham = 0.5 * (dp / (rho * cc) - dvel[dir]) * rho / cc; + const amrex::Real alphap = 0.5 * (dp / (rho * cc) + dvel[dir]) * rho / cc; + const amrex::Real alpha0r = drho - dp / cs; + const amrex::Real alpha0e = drhoe - dp * enth; + AMREX_D_TERM(, const amrex::Real alpha0v = dvel[l_idx[1]]; + , const amrex::Real alpha0w = dvel[l_idx[2]];) + const amrex::GpuArray wv = { + vel[dir] - cc, vel[dir], vel[dir] + cc}; + + // Construct the right state on the i-1/2 interface + if (area(iv) > 0.0) { + const amrex::Real rho_ref = + rho - 0.5 * (1.0 + dtdx * amrex::min(wv[0], 0.0)) * drho; + amrex::GpuArray vel_ref{ + AMREX_D_DECL(0.0, 0.0, 0.0)}; + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + vel_ref[dim] = + vel[dim] - + 0.5 * (1.0 + dtdx * amrex::min(wv[0], 0.0)) * dvel[dim]; + } + const amrex::Real p_ref = + p - 0.5 * (1.0 + dtdx * amrex::min(wv[0], 0.0)) * dp; + const amrex::Real rhoe_ref = + rhoe - 0.5 * (1.0 + dtdx * amrex::min(wv[0], 0.0)) * drhoe; + + const amrex::Real apright = 0.25 * dtdx * (wv[0] - wv[2]) * + (1.0 - amrex::Math::copysign(1.0, wv[2])) * + alphap; + const amrex::Real amright = 0.0; + + const amrex::Real azrright = 0.25 * dtdx * (wv[0] - wv[1]) * + (1.0 - amrex::Math::copysign(1.0, wv[1])) * + alpha0r; + const amrex::Real azeright = 0.25 * dtdx * (wv[0] - wv[1]) * + (1.0 - amrex::Math::copysign(1.0, wv[1])) * + alpha0e; + AMREX_D_TERM(, const amrex::Real azv1rght = + 0.25 * dtdx * (wv[0] - wv[1]) * + (1.0 - amrex::Math::copysign(1.0, wv[1])) * alpha0v; + , const amrex::Real azw1rght = + 0.25 * dtdx * (wv[0] - wv[1]) * + (1.0 - amrex::Math::copysign(1.0, wv[1])) * alpha0w;) + + qp(iv, QRHO) = rho_ref + apright + amright + azrright; + qp(iv, QRHO) = + amrex::max(qp(iv, QRHO), constants::very_small_num()); + AMREX_D_TERM(qp(iv, QU + l_idx[0]) = + vel_ref[l_idx[0]] + (apright - amright) * cc / rho; + qp(iv, QU + l_idx[1]) = 0.; qp(iv, QU + l_idx[2]) = 0.; + , qp(iv, QU + l_idx[1]) = vel_ref[l_idx[1]] + azv1rght; + , qp(iv, QU + l_idx[2]) = vel_ref[l_idx[2]] + azw1rght;); + qp(iv, QPRES) = p_ref + (apright + amright) * cs; + qp(iv, QPRES) = + amrex::max(qp(iv, QPRES), constants::very_small_num()); + qp(iv, QREINT) = rhoe_ref + (apright + amright) * enth * cs + azeright; + } else { + for (int n = 0; n < QVAR; n++) { + qp(iv, n) = 0.0; + } + } + + // Construct the left state on the i+1/2 interface + if (area(ivp) > 0.0) { + const amrex::Real rho_ref = + rho + 0.5 * (1.0 - dtdx * amrex::max(wv[2], 0.0)) * drho; + amrex::GpuArray vel_ref{ + AMREX_D_DECL(0.0, 0.0, 0.0)}; + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + vel_ref[dim] = + vel[dim] + + 0.5 * (1.0 - dtdx * amrex::max(wv[2], 0.0)) * dvel[dim]; + } + const amrex::Real p_ref = + p + 0.5 * (1.0 - dtdx * amrex::max(wv[2], 0.0)) * dp; + + const amrex::Real rhoe_ref = + rhoe + 0.5 * (1.0 - dtdx * amrex::max(wv[2], 0.0)) * drhoe; + + const amrex::Real apleft = 0.0; + const amrex::Real amleft = 0.25 * dtdx * (wv[2] - wv[0]) * + (1.0 + amrex::Math::copysign(1.0, wv[0])) * + alpham; + + const amrex::Real azrleft = 0.25 * dtdx * (wv[2] - wv[1]) * + (1.0 + amrex::Math::copysign(1.0, wv[1])) * + alpha0r; + const amrex::Real azeleft = 0.25 * dtdx * (wv[2] - wv[1]) * + (1.0 + amrex::Math::copysign(1.0, wv[1])) * + alpha0e; + AMREX_D_TERM(, const amrex::Real azv1left = + 0.25 * dtdx * (wv[2] - wv[1]) * + (1.0 + amrex::Math::copysign(1.0, wv[1])) * alpha0v; + , const amrex::Real azw1left = + 0.25 * dtdx * (wv[2] - wv[1]) * + (1.0 + amrex::Math::copysign(1.0, wv[1])) * alpha0w;) + qm(ivp, QRHO) = rho_ref + apleft + amleft + azrleft; + qm(ivp, QRHO) = + amrex::max(qm(ivp, QRHO), constants::very_small_num()); + AMREX_D_TERM(qm(ivp, QU + l_idx[0]) = + vel_ref[l_idx[0]] + (apleft - amleft) * cc / rho; + qm(ivp, QU + l_idx[1]) = 0.; qm(ivp, QU + l_idx[2]) = 0.; + , qm(ivp, QU + l_idx[1]) = vel_ref[l_idx[1]] + azv1left; + , qm(ivp, QU + l_idx[2]) = vel_ref[l_idx[2]] + azw1left;); + qm(ivp, QPRES) = p_ref + (apleft + amleft) * cs; + qm(ivp, QPRES) = + amrex::max(qm(ivp, QPRES), constants::very_small_num()); + qm(ivp, QREINT) = rhoe_ref + (apleft + amleft) * enth * cs + azeleft; + } else { + for (int n = 0; n < QVAR; n++) { + qm(ivp, n) = 0.0; + } + } + + // Upwind the passive variables + const amrex::Real vel_ad = q(iv, QU + l_idx[0]); + const amrex::Real spzerom = (vel_ad > 0) ? -1.0 : vel_ad * dtdx; + const amrex::Real spzerop = (vel_ad >= 0) ? vel_ad * dtdx : 1.0; +#if NUM_ADV > 0 + for (int n = 0; n < NUM_ADV; n++) { + pc_plm_d_eb_passive( + iv, ivp, QFA + n, spzerom, spzerop, slope, q, qp, qm, area); + } +#endif + for (int n = 0; n < NUM_SPECIES; n++) { + pc_plm_d_eb_passive( + iv, ivp, QFS + n, spzerom, spzerop, slope, q, qp, qm, area); + } +#if NUM_AUX > 0 + for (int n = 0; n < NUM_AUX; n++) { + pc_plm_d_eb_passive( + iv, ivp, QFX + n, spzerom, spzerop, slope, q, qp, qm, area); + } +#endif +#if NUM_LIN > 0 + for (int n = 0; n < NUM_LIN; n++) { + pc_plm_d_eb_passive( + iv, ivp, QLIN + n, spzerom, spzerop, slope, q, qp, qm, area); + } +#endif +} + #endif diff --git a/Source/PeleC.cpp b/Source/PeleC.cpp index 8f7cc3a47..1b7c0dde6 100644 --- a/Source/PeleC.cpp +++ b/Source/PeleC.cpp @@ -403,10 +403,6 @@ PeleC::read_params() soot_model.readSootParams(); #endif - if ((!do_mol) && eb_in_domain) { - amrex::Abort("Must do_mol = 1 when using EB\n"); - } - // TODO: What is this? amrex::StateDescriptor::setBndryFuncThreadSafety( static_cast(bndry_func_thread_safe)); diff --git a/Source/main.cpp b/Source/main.cpp index cd0b5b659..6f84835fd 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -125,7 +125,8 @@ main(int argc, char* argv[]) amrex::AmrLevel::SetEBSupportLevel( amrex::EBSupport::full); // need both area and volume fractions - amrex::AmrLevel::SetEBMaxGrowCells(7, 7, 7); + amrex::AmrLevel::SetEBMaxGrowCells( + PeleC::numGrow() + 1, PeleC::numGrow() + 1, PeleC::numGrow() + 1); initialize_EB2( amrptr->Geom(PeleC::getEBMaxLevel()), PeleC::getEBMaxLevel(), diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 34f151457..b14c6a580 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -122,10 +122,10 @@ function(add_test_rn TEST_NAME TEST_EXE_DIR) endfunction(add_test_rn) # Verification test with 1 resolution -function(add_test_v1 TEST_NAME TEST_EXE_DIR) +function(add_test_v1 TEST_NAME TEST_SCRIPT_NAME TEST_EXE_DIR) setup_test() set(RUN_COMMAND "rm -f mmslog datlog && ${MPI_COMMANDS} ${CURRENT_TEST_EXE} ${MPIEXEC_POSTFLAGS} ${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.inp") - add_test(${TEST_NAME} sh -c "${RUN_COMMAND} ${RUNTIME_OPTIONS} > ${TEST_NAME}.log && nosetests ${TEST_NAME}.py") + add_test(${TEST_NAME} sh -c "${RUN_COMMAND} ${RUNTIME_OPTIONS} > ${TEST_NAME}.log && nosetests ${TEST_SCRIPT_NAME}.py") set_tests_properties(${TEST_NAME} PROPERTIES TIMEOUT 18000 PROCESSORS ${PELE_NP} WORKING_DIRECTORY "${CURRENT_TEST_BINARY_DIR}/" LABELS "verification" ATTACHED_FILES_ON_FAIL "${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.log") endfunction(add_test_v1) @@ -223,7 +223,9 @@ add_test_r(pmf-lidryer-arkode PMF) add_test_r(pmf-srk-1 PMF-SRK) add_test_rv(masscons-mol-1 MassCons) add_test_rv(masscons-mol-2 MassCons) +add_test_rv(masscons-mol-eb MassCons) add_test_rv(masscons-plm MassCons) +add_test_rv(masscons-plm-eb MassCons) add_test_rv(masscons-ppm MassCons) add_test_rv(masscons-isothermal MassCons) add_test_r(masscons-isothermal-whydro MassCons) @@ -246,6 +248,9 @@ add_test_rv(eb-c11 EB-C11) add_test_rv(eb-c12 EB-C12) # add_test_r(eb-c14 EB-C14) # disable due to FPE in ghost cells add_test_r(eb-converging-nozzle EB-ConvergingNozzle) +if(PELE_DIM GREATER 2) + add_test_r(shock-cylinder Sod) # can run in 2D but needs input file change +endif() if(PELE_ENABLE_AMREX_PARTICLES) if(PELE_DIM EQUAL 2) add_test_spray(Spray-Conv) @@ -257,6 +262,7 @@ if(PELE_ENABLE_MASA) add_test_r(mms-2 MMS) add_test_r(mms-4 MMS) add_test_r(ebmms-1 EB-C1) + add_test_r(ebmms-2 EB-C1) if(PELE_DIM GREATER 2) add_test_r(mms-5 MMS) endif() @@ -284,7 +290,8 @@ add_test_re(spray-eb Spray-EB) add_test_re(spray-a-wbreakup Spray-A-Wbreakup) add_test_re(soot-flame Soot-Flame) add_test_re(eb-c8 EB-C8) -add_test_re(eb-c8-rere EB-C8) +add_test_re(eb-c8-rere-mol EB-C8) +add_test_re(eb-c8-rere-plm EB-C8) if(PELE_DIM GREATER 2) add_test_re(tg-3 TG) add_test_re(tg-4 TG) @@ -296,8 +303,9 @@ endif() # Verification tests #============================================================================= if(PELE_ENABLE_MASA AND (PELE_DIM GREATER 2)) - add_test_v1(symmetry MMS) - add_test_v1(eb-symmetry EB-C1) + add_test_v1(symmetry symmetry MMS) + add_test_v1(eb-symmetry-mol eb-symmetry EB-C1) + add_test_v1(eb-symmetry-plm eb-symmetry EB-C1) # Created an option to exclude these tests when using ASAN and UBSAN because they take too long if(NOT PELE_ENABLE_SANITIZE_FOR_TESTS)