ERF_InitCustomTerrain.cpp File Reference

#include "AMReX_ParmParse.H"
#include "ERF_Constants.H"
#include "ERF_TerrainMetrics.H"

Include dependency graph for ERF_InitCustomTerrain.cpp:

Functions
void	init_my_custom_terrain (const Geometry &geom, FArrayBox &terrain_fab, const double &time_d)

Function Documentation

init_my_custom_terrain()

void init_my_custom_terrain	(	const Geometry &	geom,
		FArrayBox &	terrain_fab,
		const double &	time_d
	)

Initialize a problem-specific terrain profile on the nodal terrain FAB.

Parameters

geom	Geometry defining the domain and physical extents
terrain_fab	Terrain FAB to fill with surface height data
time_d	Current simulation time used by time-dependent terrain options

{
    //
    // We put this here as a convenience for testing the map factor implementation
    // Note that these factors must match those in Source/ERF_MakeNewArrays.cpp
    //
    ParmParse pp("erf");
    bool test_mapfactor = false;
    pp.queryAdd("test_mapfactor",test_mapfactor);
 
    Real mf_m;
    if (test_mapfactor) {
        mf_m = myhalf;
    } else {
        mf_m = one;
    }
 
    // Domain cell size and real bounds
    auto dx = geom.CellSizeArray();
    auto ProbLoArr = geom.ProbLoArray();
    auto ProbHiArr = geom.ProbHiArray();
 
    const amrex::Box& domain = geom.Domain();
    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
    int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
    int domlo_z = domain.smallEnd(2);
 
    // User function parameters
    Real a    = myhalf;
    Real num  = Real(8.) * a * a * a;
    Real xcen = myhalf * (ProbLoArr[0] + ProbHiArr[0]) / mf_m;
    Real ycen = myhalf * (ProbLoArr[1] + ProbHiArr[1]) / mf_m;
 
    // Populate bottom plane
    int k0 = domlo_z;
 
    std::string custom_terrain_type = "None";
    ParmParse pp_prob("prob"); pp_prob.queryAdd("custom_terrain_type", custom_terrain_type);
 
    amrex::Box zbx = terrain_fab.box();
    if (zbx.smallEnd(2) <= k0)
    {
        amrex::Array4<Real> const& z_arr = terrain_fab.array();
 
        if (custom_terrain_type == "WoA") {
 
            // Default to x-direction
            int dir = 0;  pp_prob.queryAdd("dir", dir);
 
            Real  L        = Real(100.0); pp_prob.queryAdd("L"        , L);
            Real  z_offset = zero;        pp_prob.queryAdd("z_offset" , z_offset);
 
            // If hm is nonzero, then use alternate hill definition
            Real  hm       = zero;        pp_prob.queryAdd("hmax" , hm);
 
            // This is a 2D hill with variation in only the x-direction
            if (dir == 0) {
                ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
                {
                    // Clip indices for ghost-cells
                    int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
 
                    // Location of nodes
                    Real x = (ProbLoArr[0] + ii * dx[0] - xcen) * mf_m;
 
                    // WoA Hill in x-direction
                    if (hm==0) {
                        z_arr(i,j,k0) = num / (x*x + 4 * a * a);
                    } else {
                        Real x_L = x / L;
                        z_arr(i,j,k0) = hm / (1 + x_L*x_L) + z_offset;
                    }
                });
            } else if (dir == 1) {
                ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
                {
                    // Clip indices for ghost-cells
                    int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
 
                    // Location of nodes
                    Real y = (ProbLoArr[1] + jj * dx[1] - ycen) * mf_m;
 
                    // WoA Hill in y-direction
                    if (hm==0) {
                        z_arr(i,j,k0) = num / (y*y + Real(4.0) * a * a);
                    } else {
                        Real y_L = y / L;
                        z_arr(i,j,k0) = hm / (one + y_L*y_L) + z_offset;
                    }
                });
            } else if (dir == 2) {
                ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
                {
                    // Clip indices for ghost-cells
                    int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
                    int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
 
                    // Location of nodes
                    Real x = (ProbLoArr[0] + ii * dx[0] - xcen) * mf_m;
                    Real y = (ProbLoArr[1] + jj * dx[1] - ycen) * mf_m;
                    Real r = std::sqrt(x*x + y*y);
 
                    // WoA Hill in radial direction
                    if (hm==0) {
                        z_arr(i,j,k0) = num / (r*r + Real(4.0) * a * a);
                    } else {
                        Real r_L = r / L;
                        z_arr(i,j,k0) = hm / (one + r_L*r_L) + z_offset;
                    }
                });
            } else {
                amrex::Abort("Unknown dir in ERF_Prob.cpp");
            }
 
        } else if (custom_terrain_type == "ScharMountain") {
 
            Real a_mtn  = Real(5000.0); pp_prob.queryAdd("a"     , a_mtn);
            Real Hm     = Real(250.0);  pp_prob.queryAdd("hm"    , Hm);
            Real lambda = Real(4000.0); pp_prob.queryAdd("lambda", lambda);
            Real asq    = a_mtn * a_mtn;
 
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                // Clip indices for ghost-cells
                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
 
                // Location of nodes
                Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
 
                Real cosx = std::cos(PI * x / lambda);
 
                z_arr(i,j,k0) = Hm * std::exp(-x*x/asq) * cosx * cosx;
            });
 
        } else if (custom_terrain_type == "HalfCylinder") {
 
            Real asq = myhalf * myhalf;
 
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                // Clip indices for ghost-cells
                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
 
                // Location of nodes
                Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
 
                Real rsq = x*x;
 
                if (rsq < asq) {
                    z_arr(i,j,k0) = std::sqrt(asq - rsq);
                } else {
                    z_arr(i,j,k0) = zero;
                }
            });
 
        } else if (custom_terrain_type == "Hemisphere") {
 
            Real asq = myhalf * myhalf;
 
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                // Clip indices for ghost-cells
                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
                int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
 
                // Location of nodes
                Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
                Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
 
                Real rsq = x*x + y*y;
 
                if (rsq < asq) {
                    z_arr(i,j,k0) = std::pow(asq-rsq, myhalf);
                } else {
                    z_arr(i,j,k0) = zero;
                }
            });
 
        } else if (custom_terrain_type == "MovingSineWave") {
 
            Real Ampl        = zero;  pp_prob.queryAdd("Ampl", Ampl);
            Real wavelength  = Real(100.); pp_prob.queryAdd("wavelength", wavelength);
            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(wavelength > zero, "prob.wavelength must be > 0");
 
            Real kp          = two * PI / wavelength;
            Real g           = CONST_GRAV;
            Real omega       = std::sqrt(g * kp);
 
            Real time = static_cast<Real>(time_d);
 
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                // Clip indices for ghost-cells
                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
 
                // Location of nodes, measured from geometry.prob_lo as in the
                // MovingTerrain initial velocities
                Real x = ProbLoArr[0] + ii * dx[0];
 
                // Wave height
                Real height = Ampl * std::sin(kp * x - omega * time);
 
                // Populate terrain height at the terrain plane, as every other option does
                z_arr(i,j,k0) = height;
            });
 
        } else if (custom_terrain_type == "WindFarmTest") {
 
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                // Clip indices for ghost-cells
                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
                int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
 
                // Location of nodes
                Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
                Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
 
                Real x_L = x/Real(100.0);
                Real y_L = y/Real(100.0);
 
                z_arr(i,j,k0) = Real(100.0) / (one + x_L*x_L + y_L*y_L);
            });
 
        } else if (custom_terrain_type == "RaisedFlat") {
            Real  z_offset = zero; pp_prob.queryAdd("z_offset" , z_offset);
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                z_arr(i,j,k0) = z_offset;
            });
        } else if (custom_terrain_type == "Cos4Hill") {
 
            // Get prob parameters (must be outside GPU kernel)
            Real hm = zero; pp_prob.queryAdd("hmax", hm);
            Real L  = Real(100.0); pp_prob.queryAdd("L", L);
            Real z_offset = zero; pp_prob.queryAdd("z_offset", z_offset);
            Real fourL = four * L;
 
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
 
                // Clip indices for ghost-cells
                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
                int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
 
                // Location of nodes
                Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
                Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
                Real r = std::sqrt(x*x + y*y);
 
                if (r < fourL) {
                    z_arr(i,j,k0) = z_offset + hm * Real(0.0625) * std::pow(one + std::cos(PI*r/fourL), four);
                } else {
                    z_arr(i,j,k0) = z_offset;
                }
            });
 
        } else if (custom_terrain_type == "None") {
            // A flat surface at the bottom of the domain
            Real z_bot = ProbLoArr[2];
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                z_arr(i,j,k0) = z_bot;
            });
        } else {
            Abort("Don't know this custom_terrain_type");
        }
    }
}

Here is the call graph for this function: