Skip to content

Commit

Permalink
make xvel_bcs_dir etc multi-level (erf-model#1851)
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren authored Oct 7, 2024
1 parent 11c1678 commit 0f89138
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 76 deletions.
6 changes: 3 additions & 3 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -614,9 +614,9 @@ private:
InputSpongeData input_sponge_data;

// Vector (6 planes) of DeviceVectors (ncell in plane) for Dirichlet BC data
amrex::Gpu::DeviceVector<amrex::Real> xvel_bc_data;
amrex::Gpu::DeviceVector<amrex::Real> yvel_bc_data;
amrex::Gpu::DeviceVector<amrex::Real> zvel_bc_data;
amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> xvel_bc_data;
amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> yvel_bc_data;
amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> zvel_bc_data;

// Function to read and populate above vectors (if input file exists)
void init_Dirichlet_bc_data (const std::string input_file);
Expand Down
5 changes: 5 additions & 0 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,6 +219,11 @@ ERF::ERF_shared ()
physbcs_w.resize(nlevs_max);
physbcs_base.resize(nlevs_max);

// Planes to hold Dirichlet values at boundaries
xvel_bc_data.resize(nlevs_max);
yvel_bc_data.resize(nlevs_max);
zvel_bc_data.resize(nlevs_max);

advflux_reg.resize(nlevs_max);

// Stresses
Expand Down
16 changes: 3 additions & 13 deletions Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -550,16 +550,6 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf)
void
ERF::make_physbcs (int lev)
{
// Dirichlet BC data only lives on level 0
Real* u_bc_tmp(nullptr);
Real* v_bc_tmp(nullptr);
Real* w_bc_tmp(nullptr);
if (lev==0) {
u_bc_tmp = xvel_bc_data.data();
v_bc_tmp = yvel_bc_data.data();
w_bc_tmp = zvel_bc_data.data();
}

if (solverChoice.use_terrain) {
AMREX_ALWAYS_ASSERT(z_phys_nd[lev] != nullptr);
}
Expand All @@ -569,13 +559,13 @@ ERF::make_physbcs (int lev)
z_phys_nd[lev], use_real_bcs);
physbcs_u[lev] = std::make_unique<ERFPhysBCFunct_u> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
m_bc_extdir_vals, m_bc_neumann_vals,
z_phys_nd[lev], use_real_bcs, u_bc_tmp);
z_phys_nd[lev], use_real_bcs, xvel_bc_data[lev].data());
physbcs_v[lev] = std::make_unique<ERFPhysBCFunct_v> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
m_bc_extdir_vals, m_bc_neumann_vals,
z_phys_nd[lev], use_real_bcs, v_bc_tmp);
z_phys_nd[lev], use_real_bcs, yvel_bc_data[lev].data());
physbcs_w[lev] = std::make_unique<ERFPhysBCFunct_w> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
m_bc_extdir_vals, m_bc_neumann_vals,
solverChoice.terrain_type, z_phys_nd[lev],
use_real_bcs, w_bc_tmp);
use_real_bcs, zvel_bc_data[lev].data());
physbcs_base[lev] = std::make_unique<ERFPhysBCFunct_base> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d);
}
120 changes: 60 additions & 60 deletions Source/Initialization/ERF_init_bcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -602,63 +602,66 @@ void ERF::init_bcs ()

void ERF::init_Dirichlet_bc_data (const std::string input_file)
{
// Only on coarsest level!
int lev = 0;
const bool use_terrain = solverChoice.use_terrain;

const int klo = 0;
const int khi = geom[lev].Domain().bigEnd()[2];
const int Nz = geom[lev].Domain().size()[2];
const Real dz = geom[lev].CellSize()[2];

const bool use_terrain = (zlevels_stag[0].size() > 0);
const Real zbot = (use_terrain) ? zlevels_stag[0][klo] : geom[lev].ProbLo(2);
const Real ztop = (use_terrain) ? zlevels_stag[0][khi+1] : geom[lev].ProbHi(2);
// Read the dirichlet_input file
Print() << "dirichlet_input file location : " << input_file << std::endl;
std::ifstream input_reader(input_file);
if (!input_reader.is_open()) {
amrex::Abort("Error opening the dirichlet_input file.\n");
}

// Size of Nz (domain grid)
Vector<Real> zcc_inp(Nz );
Vector<Real> znd_inp(Nz+1);
Vector<Real> u_inp(Nz ); xvel_bc_data.resize(Nz ,0.0);
Vector<Real> v_inp(Nz ); yvel_bc_data.resize(Nz ,0.0);
Vector<Real> w_inp(Nz+1); zvel_bc_data.resize(Nz+1,0.0);
Print() << "Successfully opened the dirichlet_input file. Now reading... " << std::endl;
std::string line;

// Size of Ninp (number of z points in input file)
Vector<Real> z_inp_tmp, u_inp_tmp, v_inp_tmp, w_inp_tmp;

// Read the dirichlet_input file
Print() << "dirichlet_input file location : " << input_file << std::endl;
std::ifstream input_reader(input_file);
if(!input_reader.is_open()) {
Error("Error opening the dirichlet_input file.\n");
} else {
Print() << "Successfully opened the dirichlet_input file. Now reading... " << std::endl;
std::string line;

// First, read the input data into temp vectors; (Ninp from size of input file)

// Add surface
z_inp_tmp.push_back(zbot); // height above sea level [m]
u_inp_tmp.push_back(0.);
v_inp_tmp.push_back(0.);
w_inp_tmp.push_back(0.);

// Read the vertical profile at each given height
Real z, u, v, w;
while(std::getline(input_reader, line)) {
std::istringstream iss_z(line);
iss_z >> z >> u >> v >> w;
if (z == zbot) {
u_inp_tmp[0] = u;
v_inp_tmp[0] = v;
w_inp_tmp[0] = w;
} else {
AMREX_ALWAYS_ASSERT(z > z_inp_tmp[z_inp_tmp.size()-1]); // sounding is increasing in height
z_inp_tmp.push_back(z);
u_inp_tmp.push_back(u);
v_inp_tmp.push_back(v);
w_inp_tmp.push_back(w);
if (z >= ztop) break;
}
const int klo = geom[0].Domain().smallEnd()[2];
const int khi = geom[0].Domain().bigEnd()[2];

const Real zbot = (use_terrain) ? zlevels_stag[0][klo] : geom[0].ProbLo(2);
const Real ztop = (use_terrain) ? zlevels_stag[0][khi+1] : geom[0].ProbHi(2);

// Add surface
z_inp_tmp.push_back(zbot); // height above sea level [m]
u_inp_tmp.push_back(0.);
v_inp_tmp.push_back(0.);
w_inp_tmp.push_back(0.);

// Read the vertical profile at each given height
Real z, u, v, w;
while(std::getline(input_reader, line)) {
std::istringstream iss_z(line);
iss_z >> z >> u >> v >> w;
if (z == zbot) {
u_inp_tmp[0] = u;
v_inp_tmp[0] = v;
w_inp_tmp[0] = w;
} else {
AMREX_ALWAYS_ASSERT(z > z_inp_tmp[z_inp_tmp.size()-1]); // sounding is increasing in height
z_inp_tmp.push_back(z);
u_inp_tmp.push_back(u);
v_inp_tmp.push_back(v);
w_inp_tmp.push_back(w);
if (z >= ztop) break;
}
}

amrex::Print() << "Successfully read and interpolated the dirichlet_input file..." << std::endl;
input_reader.close();

for (int lev = 0; lev <= max_level; lev++) {

const int Nz = geom[lev].Domain().size()[2];
const Real dz = geom[lev].CellSize()[2];

// Size of Nz (domain grid)
Vector<Real> zcc_inp(Nz );
Vector<Real> znd_inp(Nz+1);
Vector<Real> u_inp(Nz ); xvel_bc_data[lev].resize(Nz ,0.0);
Vector<Real> v_inp(Nz ); yvel_bc_data[lev].resize(Nz ,0.0);
Vector<Real> w_inp(Nz+1); zvel_bc_data[lev].resize(Nz+1,0.0);

// At this point, we have an input from zbot up to
// z_inp_tmp[N-1] >= ztop. Now, interpolate to grid level 0 heights
Expand All @@ -673,16 +676,13 @@ void ERF::init_Dirichlet_bc_data (const std::string input_file)
}
znd_inp[Nz] = ztop;
w_inp[Nz] = interpolate_1d(z_inp_tmp.dataPtr(), w_inp_tmp.dataPtr(), ztop, Ninp);
}

amrex::Print() << "Successfully read and interpolated the dirichlet_input file..." << std::endl;
input_reader.close();

// Copy host data to the device
Gpu::copy(Gpu::hostToDevice, u_inp.begin(), u_inp.end(), xvel_bc_data.begin());
Gpu::copy(Gpu::hostToDevice, v_inp.begin(), v_inp.end(), yvel_bc_data.begin());
Gpu::copy(Gpu::hostToDevice, w_inp.begin(), w_inp.end(), zvel_bc_data.begin());
// Copy host data to the device
Gpu::copy(Gpu::hostToDevice, u_inp.begin(), u_inp.end(), xvel_bc_data[lev].begin());
Gpu::copy(Gpu::hostToDevice, v_inp.begin(), v_inp.end(), yvel_bc_data[lev].begin());
Gpu::copy(Gpu::hostToDevice, w_inp.begin(), w_inp.end(), zvel_bc_data[lev].begin());

// NOTE: These device vectors are passed to the PhysBC constructors when that
// class is instantiated in ERF_make_new_arrays.cpp.
// NOTE: These device vectors are passed to the PhysBC constructors when that
// class is instantiated in ERF_make_new_arrays.cpp.
} // lev
}

0 comments on commit 0f89138

Please sign in to comment.