diff --git a/Docs/sphinx_doc/MeshRefinement.rst b/Docs/sphinx_doc/MeshRefinement.rst index 2a0ae0a5d..aa385172f 100644 --- a/Docs/sphinx_doc/MeshRefinement.rst +++ b/Docs/sphinx_doc/MeshRefinement.rst @@ -124,7 +124,7 @@ computed by dividing the variable named rhotheta by the variable named density. Coupling Types -------------- -ERF supports one-way, two-way, and "mixed" coupling; this is a run-time input +ERF supports one-way, two-way, and "mixed" coupling between levels; this is a run-time input :: @@ -136,8 +136,9 @@ for the time advance of the fine solution . For cell-centered quantities, and face-baced normal momenta on the coarse-fine interface, the coarse data is conservatively interpolated to the fine mesh. The interpolated data is utilized to specify ghost cell data (outside of the valid fine region) as well as specified and relaxation data inside the lateral boundaries -of the fine region. More specifically, a user may specify the total width of the interior -Dirichlet and relaxation region with ``erf.cf_width = `` (yellow + blue) +of the fine region. More specifically, similarly to how the lateral boundaries are treated, +a user may specify the total width of the interior Dirichlet and relaxation region with +``erf.cf_width = `` (yellow + blue) and analogously the width of the interior Dirichlet region may be specified with ``erf.cf_set_width = `` (yellow). @@ -173,7 +174,9 @@ where :math:`G` is the RHS of the evolution equations, :math:`\psi^{\prime}` is relaxation, :math:`\psi^{FP}` is the fine data obtained from spatial and temporal interpolation of the coarse data, and :math:`n` is the minimum number of grid points from a lateral boundary. The specified and relaxation regions are applied to all dycore variables :math:`\left[\rho \; \rho\Theta \; U\; V\; W \right]` -on the fine mesh. Finally, we note that time dependent Dirichlet data, provided via an external boundary file, +on the fine mesh. + +Finally, we note that time dependent Dirichlet data, provided via an external boundary file, may be enforced on the lateral boundary conditions of the domain (coarsest mesh). For such cases, the relaxation region width at the domain edges may be specified with ``erf.wrfbdy_width = `` (yellow + blue) while the interior Dirichlet region may be specified with ``erf.wrfbdy_set_width = `` @@ -187,17 +190,11 @@ the fine mesh communicates data back to the coarse mesh in two ways: - A "reflux" operation is performed for all cell-centered data. -Because the normal momentum at the fine level is itself interpolated from the coarse level, the -difference between fluxes -- along the coarse-fine interfaces -- used to update the coarse data and fluxes -used to update the fine data is due to the difference in the averaging of the advected quantity to the face -where the flux is defined. - -We note that both coupling schemes are conservative for mass because the fluxes for the continuity -equation are the momenta themselves, which are interpolated on faces at the coarse-fine interface. Other advected -quantities which are advanced in conservation form will lose conservation with one-way coupling. -Two-way coupling is conservative for these scalars as long as the refluxing operation is included with the -averaging down. - We define "mixed" coupling as using the two-way coupling algorithm for all cell-centered quantities except for :math:`\rho` and :math:`\rho \theta.` +We note that all three coupling schemes are conservative for mass because the fluxes for the continuity +equation are the momenta themselves, which are interpolated on faces at the coarse-fine interface. Other advected +quantities which are advanced in conservation form will lose conservation with one-way coupling. +Two-way coupling ensures conservation of the advective contribution to all scalar updates but +does not account for loss of conservation due to diffusive or source terms. diff --git a/Source/Advection/AdvectionSrcForScalars.H b/Source/Advection/AdvectionSrcForScalars.H index 5b99c5cbf..c09f15c81 100644 --- a/Source/Advection/AdvectionSrcForScalars.H +++ b/Source/Advection/AdvectionSrcForScalars.H @@ -72,28 +72,28 @@ AdvectionSrcForScalarsVert(const amrex::Box& bx, switch(vert_adv_type) { case AdvType::Centered_2nd: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, - flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Upwind_3rd: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, - flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Centered_4th: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, - flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Upwind_5th: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, - flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Centered_6th: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, - flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); break; default: AMREX_ASSERT_WITH_MESSAGE(false, "Unknown vertical advection scheme!"); diff --git a/Source/Advection/AdvectionSrcForState.cpp b/Source/Advection/AdvectionSrcForState.cpp index 286a4fecd..015107e59 100644 --- a/Source/Advection/AdvectionSrcForState.cpp +++ b/Source/Advection/AdvectionSrcForState.cpp @@ -176,43 +176,43 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, switch(horiz_adv_type) { case AdvType::Centered_2nd: AdvectionSrcForScalarsVert(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); break; case AdvType::Upwind_3rd: AdvectionSrcForScalarsVert(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); break; case AdvType::Centered_4th: AdvectionSrcForScalarsVert(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); break; case AdvType::Upwind_5th: AdvectionSrcForScalarsVert(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); break; case AdvType::Centered_6th: AdvectionSrcForScalarsVert(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); break; case AdvType::Weno_3: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Weno_5: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Weno_3Z: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Weno_3MZQ: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + avg_xmom, avg_ymom, avg_zmom); break; case AdvType::Weno_5Z: AdvectionSrcForScalarsWrapper(bx, ncomp, icomp, flx_arr, cell_prim, - avg_xmom, avg_ymom, avg_zmom); + avg_xmom, avg_ymom, avg_zmom); break; default: AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); diff --git a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp index afb303590..4fe58d7f0 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp @@ -553,7 +553,7 @@ void erf_fast_rhs_MT (int step, int /*level*/, // Note that in the solve we effectively impose new_drho_w(i,j,vbx_hi.z+1)=0 // so we don't update avg_zmom at k=vbx_hi.z+1 - avg_zmom(i,j,k) += facinv*zflux_lo; + avg_zmom(i,j,k) += facinv*zflux_lo / (mf_m(i,j,0) * mf_m(i,j,0)); // Note that the factor of (1/J) in the fast source term is canceled // when we multiply old and new by detJ_old and detJ_new , respectively diff --git a/Source/TimeIntegration/ERF_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index 98cde76a9..ab3e13279 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -484,7 +484,7 @@ void erf_fast_rhs_N (int step, int /*level*/, Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * old_drho_w(i,j,k ); Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * old_drho_w(i,j,k+1); - avg_zmom(i,j,k) += facinv*zflux_lo; + avg_zmom(i,j,k) += facinv*zflux_lo / (mf_m(i,j,0) * mf_m(i,j,0)); // Note that in the solve we effectively impose soln_a(i,j,vbx_hi.z+1)=0 // so we don't update avg_zmom at k=vbx_hi.z+1 diff --git a/Source/TimeIntegration/ERF_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index 8cdd9be40..d42560d1c 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -584,7 +584,7 @@ void erf_fast_rhs_T (int step, int /*level*/, // Note that in the solve we effectively impose new_drho_w(i,j,vbx_hi.z+1)=0 // so we don't update avg_zmom at k=vbx_hi.z+1 - avg_zmom(i,j,k) += facinv*zflux_lo; + avg_zmom(i,j,k) += facinv*zflux_lo / (mf_m(i,j,0) * mf_m(i,j,0)); Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k); cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho);