diff --git a/Source/ERF.H b/Source/ERF.H index d5ed5053d..c69751c69 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -614,9 +614,9 @@ private: InputSpongeData input_sponge_data; // Vector (6 planes) of DeviceVectors (ncell in plane) for Dirichlet BC data - amrex::Gpu::DeviceVector xvel_bc_data; - amrex::Gpu::DeviceVector yvel_bc_data; - amrex::Gpu::DeviceVector zvel_bc_data; + amrex::Vector> xvel_bc_data; + amrex::Vector> yvel_bc_data; + amrex::Vector> zvel_bc_data; // Function to read and populate above vectors (if input file exists) void init_Dirichlet_bc_data (const std::string input_file); diff --git a/Source/ERF.cpp b/Source/ERF.cpp index b4bc194e3..c0415401f 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -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 diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index f77624a2d..434d42c10 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -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); } @@ -569,13 +559,13 @@ ERF::make_physbcs (int lev) z_phys_nd[lev], use_real_bcs); physbcs_u[lev] = std::make_unique (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 (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 (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 (lev, geom[lev], domain_bcs_type, domain_bcs_type_d); } diff --git a/Source/Initialization/ERF_init_bcs.cpp b/Source/Initialization/ERF_init_bcs.cpp index 11b78be73..449349b06 100644 --- a/Source/Initialization/ERF_init_bcs.cpp +++ b/Source/Initialization/ERF_init_bcs.cpp @@ -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 zcc_inp(Nz ); - Vector znd_inp(Nz+1); - Vector u_inp(Nz ); xvel_bc_data.resize(Nz ,0.0); - Vector v_inp(Nz ); yvel_bc_data.resize(Nz ,0.0); - Vector 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 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 zcc_inp(Nz ); + Vector znd_inp(Nz+1); + Vector u_inp(Nz ); xvel_bc_data[lev].resize(Nz ,0.0); + Vector v_inp(Nz ); yvel_bc_data[lev].resize(Nz ,0.0); + Vector 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 @@ -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 }