Skip to content

Commit

Permalink
remove files with duplicate code (#1504)
Browse files Browse the repository at this point in the history 
* remove files with duplicate code

* fix oops
  • Loading branch information
@asalmgren
asalmgren authored Mar 18, 2024
1 parent 8bd9c3f commit 03f3255
Show file tree
Hide file tree
Showing 11 changed files with 87 additions and 177 deletions.
2 changes: 1 addition & 1 deletion Exec/DevTests/MovingTerrain/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ class Problem : public ProblemBase
public:
Problem();

#include "Prob/init_density_hse_dry_terrain.H"
#include "Prob/init_density_hse_dry.H"

void init_custom_pert (
const amrex::Box& bx,
Expand Down
2 changes: 1 addition & 1 deletion Exec/DevTests/MultiBlock/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ class Problem : public ProblemBase
public:
Problem();

#include "Prob/init_density_hse_dry_noterrain.H"
#include "Prob/init_density_hse_dry.H"

void init_custom_pert (
const amrex::Box& bx,
Expand Down
2 changes: 1 addition & 1 deletion Exec/DevTests/ParticlesOverWoA/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ class Problem : public ProblemBase
public:
Problem();

#include "Prob/init_density_hse_dry_terrain.H"
#include "Prob/init_density_hse_dry.H"

void init_custom_pert (
const amrex::Box& bx,
Expand Down
2 changes: 1 addition & 1 deletion Exec/RegTests/Terrain2d_Cylinder/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ class Problem : public ProblemBase
public:
Problem();

#include "Prob/init_density_hse_dry_terrain.H"
#include "Prob/init_density_hse_dry.H"

void init_custom_pert (
const amrex::Box& bx,
Expand Down
2 changes: 1 addition & 1 deletion Exec/RegTests/Terrain3d_Hemisphere/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ class Problem : public ProblemBase
public:
Problem (const amrex::Real* problo, const amrex::Real* probhi);

#include "Prob/init_density_hse_dry_terrain.H"
#include "Prob/init_density_hse_dry.H"

void init_custom_pert (
const amrex::Box& bx,
Expand Down
2 changes: 1 addition & 1 deletion Exec/RegTests/WitchOfAgnesi/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ class Problem : public ProblemBase
public:
Problem();

#include "Prob/init_density_hse_dry_terrain.H"
#include "Prob/init_density_hse_dry.H"
#include "Prob/init_rayleigh_damping.H"

void init_custom_pert (
Expand Down
142 changes: 72 additions & 70 deletions Source/HSE_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,88 +138,90 @@ namespace HSEutils
amrex::Real* r,
amrex::Real* p,
const amrex::Array4<amrex::Real const> z_cc,
const int& khi)
{
// r_sfc / p_0 are the density / pressure at the surface
int k0 = 0;
Real half_dz = z_cc(i,j,k0);

// Initial guess
r[k0] = r_sfc;
p[k0] = p_0 - half_dz * r[k0] * CONST_GRAV;

{
// We do a Newton iteration to satisfy the EOS & HSE (with constant theta)
bool converged_hse = false;
Real p_hse;
Real p_eos;

for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
{
p_hse = p_0 - half_dz * r[k0] * CONST_GRAV;
p_eos = getPgivenRTh(r[k0]*theta);

Real A = p_hse - p_eos;

Real dpdr = getdPdRgivenConstantTheta(r[k0],theta);

Real drho = A / (dpdr + half_dz * CONST_GRAV);

r[k0] = r[k0] + drho;
p[k0] = getPgivenRTh(r[k0]*theta);

if (std::abs(drho) < TOL)
{
converged_hse = true;
break;
}
}
const int& klo, const int& khi)
{
if (klo == 0) {
//
// r_sfc / p_0 are the density / pressure at the surface
//
// Initial guess
int k0 = 0;
Real half_dz = z_cc(i,j,k0);
r[k0] = r_sfc;
p[k0] = p_0 - half_dz * r[k0] * CONST_GRAV;
{
// We do a Newton iteration to satisfy the EOS & HSE (with constant theta)
bool converged_hse = false;
Real p_hse;
Real p_eos;

for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
{
p_hse = p_0 - half_dz * r[k0] * CONST_GRAV;
p_eos = getPgivenRTh(r[k0]*theta);

Real A = p_hse - p_eos;

Real dpdr = getdPdRgivenConstantTheta(r[k0],theta);

Real drho = A / (dpdr + half_dz * CONST_GRAV);

r[k0] = r[k0] + drho;
p[k0] = getPgivenRTh(r[k0]*theta);

if (std::abs(drho) < TOL)
{
converged_hse = true;
break;
}
}

//if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << k0 << std::endl;
//if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
}
}

//if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << k0 << std::endl;
//if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
}
// To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
for (int k = klo+1; k <= khi; k++)
{
// To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
// to discretely satisfy HSE -- here we assume spatial_order = 2 -- we can generalize this later if needed
bool converged_hse = false;

// To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
for (int k = 1; k <= khi; k++)
{
// To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
// to discretely satisfy HSE -- here we assume spatial_order = 2 -- we can generalize this later if needed
bool converged_hse = false;
Real dz_loc = (z_cc(i,j,k) - z_cc(i,j,k-1));

Real dz_loc = (z_cc(i,j,k) - z_cc(i,j,k-1));
r[k] = r[k-1];

r[k] = r[k-1];
Real p_eos = getPgivenRTh(r[k]*theta);
Real p_hse;

Real p_eos = getPgivenRTh(r[k]*theta);
Real p_hse;
for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
{
p_hse = p[k-1] - dz_loc * 0.5 * (r[k-1]+r[k]) * CONST_GRAV;
p_eos = getPgivenRTh(r[k]*theta);

for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
{
Real A = p_hse - p_eos;

p_hse = p[k-1] - dz_loc * 0.5 * (r[k-1]+r[k]) * CONST_GRAV;
p_eos = getPgivenRTh(r[k]*theta);
Real dpdr = getdPdRgivenConstantTheta(r[k],theta);

Real A = p_hse - p_eos;

Real dpdr = getdPdRgivenConstantTheta(r[k],theta);
Real drho = A / (dpdr + dz_loc * CONST_GRAV);

Real drho = A / (dpdr + dz_loc * CONST_GRAV);
r[k] = r[k] + drho;
p[k] = getPgivenRTh(r[k]*theta);

r[k] = r[k] + drho;
p[k] = getPgivenRTh(r[k]*theta);
if (std::abs(drho) < TOL * r[k-1])
{
converged_hse = true;
//amrex::Print() << " converged " << std::endl;
break;
}
}

if (std::abs(drho) < TOL * r[k-1])
{
converged_hse = true;
//amrex::Print() << " converged " << std::endl;
break;
}
}

//if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << k << std::endl;
//if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
}
// if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << k << std::endl;
// if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
}
}

} // namespace

#endif
2 changes: 1 addition & 1 deletion Source/Microphysics/LagrangianMicrophysics.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ public:
/*! \brief Advance the moisture model for one time step */
void Advance ( const int& lev, /*!< AMR level */
const amrex::Real& dt_advance, /*!< Time step */
const SolverChoice &solverChoice, /*!< Solver choice object */
const SolverChoice& /*solverChoice*/, /*!< Solver choice object */
amrex::Vector<amrex::Vector<amrex::MultiFab>>& a_vars, /*!< Dycore state variables */
const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& a_z/*!< terrain */) override
{
Expand Down
11 changes: 8 additions & 3 deletions Source/Prob/init_density_hse_dry.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ erf_init_dens_hse (amrex::MultiFab& rho_hse,
{
const amrex::Real prob_lo_z = geom.ProbLo()[2];
const amrex::Real dz = geom.CellSize()[2];
const int khi = geom.Domain().bigEnd()[2];

const amrex::Real T_sfc = 300.;
const amrex::Real rho_sfc = p_0 / (R_d*T_sfc);
Expand All @@ -24,7 +23,7 @@ erf_init_dens_hse (amrex::MultiFab& rho_hse,
// use_terrain = 1
if (z_phys_nd) {

if (khi > 255) amrex::Abort("1D Arrays are hard-wired to only 256 high");
if (geom.Domain().bigEnd(2) > 255) amrex::Abort("1D Arrays are hard-wired to only 256 high");

for ( amrex::MFIter mfi(rho_hse, TileNoZ()); mfi.isValid(); ++mfi )
{
Expand All @@ -37,11 +36,14 @@ erf_init_dens_hse (amrex::MultiFab& rho_hse,
b2d.grow(0,1); b2d.grow(1,1); // Grow by one in the lateral directions
b2d.setRange(2,0);

const int klo = tbz.smallEnd(2);
const int khi = tbz.bigEnd(2)-1;

amrex::ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
amrex::Array1D<amrex::Real,0,255> r;;
amrex::Array1D<amrex::Real,0,255> p;;

HSEutils::init_isentropic_hse_terrain(i,j,rho_sfc,Thetabar,&(r(0)),&(p(0)),z_cc_arr,khi);
HSEutils::init_isentropic_hse_terrain(i,j,rho_sfc,Thetabar,&(r(0)),&(p(0)),z_cc_arr,klo,khi);

for (int k = 0; k <= khi; k++) {
rho_arr(i,j,k) = r(k);
Expand All @@ -53,6 +55,9 @@ erf_init_dens_hse (amrex::MultiFab& rho_hse,

} else { // use_terrain = 0

const int klo = 0;
const int khi = geom.Domain().bigEnd(2);

// These are at cell centers (unstaggered)
amrex::Vector<amrex::Real> h_r(khi+2);
amrex::Vector<amrex::Real> h_p(khi+2);
Expand Down
50 changes: 0 additions & 50 deletions Source/Prob/init_density_hse_dry_noterrain.H
View file

This file was deleted.

47 changes: 0 additions & 47 deletions Source/Prob/init_density_hse_dry_terrain.H
View file

This file was deleted.

0 comments on commit 03f3255

Please sign in to comment.