Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove flux limiter on large velocities #2132

Merged
merged 2 commits into from
Apr 24, 2022
Merged

Remove flux limiter on large velocities #2132

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
37a6908
Remove flux limiter on large velocities
maxpkatz Apr 24, 2022
cd2cde5
Update CHANGES
maxpkatz Apr 24, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
# 22.05

* The option castro.limit_fluxes_on_large_vel has been removed. (#2132)

# 22.04

* We now abort on GPUs if species do not sum to 1 (#2099)
Expand Down
3 changes: 0 additions & 3 deletions Source/driver/_cpp_parameters
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,9 +167,6 @@ pslope_cutoff_density Real -1.e20
# Should we limit the density fluxes so that we do not create small densities?
limit_fluxes_on_small_dens int 0

# Should we limit the momentum fluxes so that we do not create large velocities?
limit_fluxes_on_large_vel int 0

# Enforce the magnitude of the velocity to be no larger than this number (and
# optionally limit the fluxes as well). Only applies if it is greater than 0.
speed_limit Real 0.0
Expand Down
11 changes: 0 additions & 11 deletions Source/hydro/Castro_ctu_hydro.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1215,17 +1215,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
dt);
}

if (limit_fluxes_on_large_vel == 1) {
limit_hydro_fluxes_on_large_vel
(nbx, idir,
Sborder.array(mfi),
q.array(),
volume.array(mfi),
flux[idir].array(),
area[idir].array(mfi),
dt);
}

normalize_species_fluxes(nbx, flux_arr);

}
Expand Down
10 changes: 0 additions & 10 deletions Source/hydro/Castro_hydro.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -781,16 +781,6 @@
amrex::Array4<amrex::Real const> const& area,
amrex::Real dt);

void
limit_hydro_fluxes_on_large_vel(const amrex::Box& bx,
int idir,
amrex::Array4<amrex::Real const> const& u,
amrex::Array4<amrex::Real const> const& q,
amrex::Array4<amrex::Real const> const& vol,
amrex::Array4<amrex::Real> const& flux,
amrex::Array4<amrex::Real const> const& area,
amrex::Real dt);

void scale_flux(const amrex::Box& bx,
#if AMREX_SPACEDIM == 1
amrex::Array4<amrex::Real const> const& qint,
Expand Down
184 changes: 0 additions & 184 deletions Source/hydro/advection_util.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -893,190 +893,6 @@ Castro::limit_hydro_fluxes_on_small_dens(const Box& bx,



void
Castro::limit_hydro_fluxes_on_large_vel(const Box& bx,
int idir,
Array4<Real const> const& u,
Array4<Real const> const& q,
Array4<Real const> const& vol,
Array4<Real> const& flux,
Array4<Real const> const& area_arr,
Real dt)
{

// This limiter is similar to the density-based limiter above, but limits
// on velocities that are too large instead. The comments are minimal since
// the algorithm is effectively the same.

if (castro::speed_limit <= 0.0_rt) return;

const Real* dx = geom.CellSize();

Real dtdx = dt / dx[idir];
Real lcfl = cfl;
Real alpha = 1.0_rt / AMREX_SPACEDIM;

auto coord = geom.Coord();
GeometryData geomdata = geom.data();

Real lspeed_limit = speed_limit / (2 * AMREX_SPACEDIM);

amrex::ParallelFor(bx,
[=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
{
GpuArray<Real, NUM_STATE> uR;
for (int n = 0; n < NUM_STATE; ++n) {
uR[n] = u(i,j,k,n);
}

GpuArray<Real, NQ> qR;
for (int n = 0; n < NQ; ++n) {
qR[n] = q(i,j,k,n);
}

Real volR = vol(i,j,k);

GpuArray<int, 3> idxR = {i,j,k};

GpuArray<Real, NUM_STATE> uL;
GpuArray<Real, NQ> qL;
Real volL;
GpuArray<int, 3> idxL;

if (idir == 0) {

for (int n = 0; n < NUM_STATE; ++n) {
uL[n] = u(i-1,j,k,n);
}

for (int n = 0; n < NQ; ++n) {
qL[n] = q(i-1,j,k,n);
}

volL = vol(i-1,j,k);

idxL = {i-1,j,k};

}
else if (idir == 1) {

for (int n = 0; n < NUM_STATE; ++n) {
uL[n] = u(i,j-1,k,n);
}

for (int n = 0; n < NQ; ++n) {
qL[n] = q(i,j-1,k,n);
}

volL = vol(i,j-1,k);

idxL = {i,j-1,k};

}
else {

for (int n = 0; n < NUM_STATE; ++n) {
uL[n] = u(i,j,k-1,n);
}

for (int n = 0; n < NQ; ++n) {
qL[n] = q(i,j,k-1,n);
}

volL = vol(i,j,k-1);

idxL = {i,j,k-1};

}

// Construct cell-centered fluxes.

GpuArray<Real, NUM_STATE> fluxL;
dflux(uL, qL, idir, coord, geomdata, idxL, fluxL);

GpuArray<Real, NUM_STATE> fluxR;
dflux(uR, qR, idir, coord, geomdata, idxR, fluxR);

// Construct the Lax-Friedrichs flux on the interface.

GpuArray<Real, NUM_STATE> fluxLF;
for (int n = 0; n < NUM_STATE; ++n) {
fluxLF[n] = 0.5_rt * (fluxL[n] + fluxR[n] + (lcfl / dtdx / alpha) * (uL[n] - uR[n]));
}

// Coefficients of fluxes on either side of the interface.

Real flux_coefR = 2.0_rt * (dt / alpha) * area_arr(i,j,k) / volR;
Real flux_coefL = 2.0_rt * (dt / alpha) * area_arr(i,j,k) / volL;

Real theta = 1.0_rt;

// Loop over all three momenta, and choose the strictest
// limiter among them.

for (int n = 0; n < 3; ++n) {

int UMOM = UMX + n;

// Obtain the one-sided update to the momentum.

Real drhouL = flux_coefL * flux(i,j,k,UMOM);
Real rhouL = std::abs(uL[UMOM] - drhouL);

Real drhoL = flux_coefL * flux(i,j,k,URHO);
Real rhoL = uL[URHO] - drhoL;

Real drhouR = flux_coefR * flux(i,j,k,UMOM);
Real rhouR = std::abs(uR[UMOM] + drhouR);

Real drhoR = flux_coefR * flux(i,j,k,URHO);
Real rhoR = uR[URHO] + drhoR;

if (std::abs(rhouL) > rhoL * lspeed_limit) {

// Obtain the final density corresponding to the LF flux.

Real drhouLF = flux_coefL * fluxLF[UMOM];
Real rhouLF = std::abs(uL[UMOM] - drhouLF);

// Solve for theta from (1 - theta) * rhouLF + theta * rhou = rhoL * speed_limit.

theta = amrex::min(theta, std::abs(rhoL * lspeed_limit - rhouLF) / std::abs(rhouL - rhouLF));

}
else if (std::abs(rhouR) > rhoR * lspeed_limit) {

Real drhouLF = flux_coefR * fluxLF[UMOM];
Real rhouLF = std::abs(uR[UMOM] + drhouLF);

theta = amrex::min(theta, std::abs(rhoR * lspeed_limit - rhouLF) / std::abs(rhouR - rhouLF));

}

}

// Limit theta to the valid range (this will deal with roundoff issues).

theta = amrex::min(1.0_rt, amrex::max(theta, 0.0_rt));

// Assemble the limited flux (Equation 16).

for (int n = 0; n < NUM_STATE; ++n) {
flux(i,j,k,n) = (1.0_rt - theta) * fluxLF[n] + theta * flux(i,j,k,n);
}

// Zero out fluxes for quantities that don't advect.

flux(i,j,k,UTEMP) = 0.0_rt;
#ifdef SHOCK_VAR
flux(i,j,k,USHK) = 0.0_rt;
#endif

});

}


void
Castro::do_enforce_minimum_density(const Box& bx,
Array4<Real> const& state_arr,
Expand Down