diff --git a/cases/dust_settling/README.md b/cases/drycblles_particle_bin/README.md similarity index 100% rename from cases/dust_settling/README.md rename to cases/drycblles_particle_bin/README.md diff --git a/cases/dust_settling/dust.ini.base b/cases/drycblles_particle_bin/particle_bin.ini.base similarity index 94% rename from cases/dust_settling/dust.ini.base rename to cases/drycblles_particle_bin/particle_bin.ini.base index 770c26dbe..13d45c273 100644 --- a/cases/dust_settling/dust.ini.base +++ b/cases/drycblles_particle_bin/particle_bin.ini.base @@ -57,10 +57,10 @@ flow_direction[north]=outflow flow_direction[east]=outflow flow_direction[south]=outflow -[dust] -swdust=1 -dustlist=None -w_terminal=None +[particle_bin] +sw_particle=1 +particle_list=None +w_particle=None [fields] visc=1.e-5 diff --git a/cases/dust_settling/dust_input.py b/cases/drycblles_particle_bin/particle_bin_input.py similarity index 70% rename from cases/dust_settling/dust_input.py rename to cases/drycblles_particle_bin/particle_bin_input.py index 99bf52475..f556a4aea 100644 --- a/cases/dust_settling/dust_input.py +++ b/cases/drycblles_particle_bin/particle_bin_input.py @@ -9,7 +9,7 @@ """ float_type = np.float32 # np.float32 for -USESP=true, else np.float64. -dust_bins = np.array([0, 2, 10, 20, 58, 83, 440]) # (μm) +particle_bins = np.array([0, 2, 10, 20, 58, 83, 440]) # (μm) xsize = 12800 ysize = 6400 @@ -42,10 +42,10 @@ # All scalars start at concentration zero, with # inflow of air with concentration zero at lateral boundaries. -dust_list = [f'{dust_bins[i]}-{dust_bins[i+1]}um' for i in range(dust_bins.size - 1)] +particle_list = [f'{particle_bins[i]}-{particle_bins[i+1]}um' for i in range(particle_bins.size - 1)] scalars = {} -for scalar in dust_list: +for scalar in particle_list: scalars[scalar] = np.zeros(ktot) @@ -58,9 +58,9 @@ nu = 1e-5 # Kinematic viscosity air [m2 s-1] g = 9.81 # Gravitational acceleration [m s-2] -dust_diameter = 0.5*(dust_bins[1:] + dust_bins[:-1]) * 1e-6 +particle_diameter = 0.5*(particle_bins[1:] + particle_bins[:-1]) * 1e-6 -tau_p = dust_diameter**2 * rho_p / (18 * nu * rho_a) +tau_p = particle_diameter**2 * rho_p / (18 * nu * rho_a) w_terminal = -tau_p * g # Create circular field with dust emissions. @@ -74,14 +74,14 @@ field_flux = np.zeros((jtot, itot), dtype=float_type) field_flux[field_mask] = 1. -for scalar in dust_list: +for scalar in particle_list: field_flux.tofile('{}_bot_in.0000000'.format(scalar)) """ Set/write new namelist. """ -ini = mht.Read_namelist('dust.ini.base') +ini = mht.Read_namelist('particle_bin.ini.base') ini['grid']['itot'] = itot ini['grid']['jtot'] = jtot @@ -95,23 +95,23 @@ ini['time']['endtime'] = endtime -ini['fields']['slist'] = dust_list -ini['advec']['fluxlimit_list'] = dust_list -ini['limiter']['limitlist'] = dust_list -ini['boundary']['scalar_outflow'] = dust_list -ini['boundary']['sbot_2d_list'] = dust_list +ini['fields']['slist'] = particle_list +ini['advec']['fluxlimit_list'] = particle_list +ini['limiter']['limitlist'] = particle_list +ini['boundary']['scalar_outflow'] = particle_list +ini['boundary']['sbot_2d_list'] = particle_list -ini['dust']['dustlist'] = dust_list -for i in range(len(dust_list)): - ini['dust'][f'w_terminal[{dust_list[i]}]'] = w_terminal[i] +ini['particle_bin']['particle_list'] = particle_list +for i in range(len(particle_list)): + ini['particle_bin'][f'w_particle[{particle_list[i]}]'] = w_terminal[i] # Statistics/crosses/... -scalar_crosses = dust_list + [s+'_path' for s in dust_list] +scalar_crosses = particle_list + [s+'_path' for s in particle_list] ini['cross']['crosslist'] = scalar_crosses + ['th', 'u', 'v', 'w'] ini['cross']['xz'] = y0 ini['cross']['yz'] = x0 -ini.save('dust.ini', allow_overwrite=True) +ini.save('particle_bin.ini', allow_overwrite=True) """ @@ -121,7 +121,7 @@ def add_var(name, dims, values, nc_group): nc_var = nc_group.createVariable(name, float_type, dims) nc_var[:] = values -nc_file = nc.Dataset('dust_input.nc', mode='w', datamodel='NETCDF4') +nc_file = nc.Dataset('particle_bin_input.nc', mode='w', datamodel='NETCDF4') nc_file.createDimension('z', ktot) add_var('z', ('z'), z, nc_file) diff --git a/cases/dust_settling/plot_sedimentation_velocity.py b/cases/drycblles_particle_bin/plot_sedimentation_velocity.py similarity index 100% rename from cases/dust_settling/plot_sedimentation_velocity.py rename to cases/drycblles_particle_bin/plot_sedimentation_velocity.py diff --git a/cases/dust_settling/sedimentation_velocity.png b/cases/drycblles_particle_bin/sedimentation_velocity.png similarity index 100% rename from cases/dust_settling/sedimentation_velocity.png rename to cases/drycblles_particle_bin/sedimentation_velocity.png diff --git a/include/model.h b/include/model.h index 236354848..f83d1da38 100644 --- a/include/model.h +++ b/include/model.h @@ -46,7 +46,7 @@ template class Pres; template class Force; template class Aerosol; template class Background; -template class Dust; +template class Particle_bin; template class Thermo; template class Microphys; template class Radiation; @@ -107,7 +107,8 @@ class Model std::shared_ptr> decay; std::shared_ptr> limiter; std::shared_ptr> source; - std::shared_ptr> dust; + + std::shared_ptr> particle_bin; std::shared_ptr> stats; std::shared_ptr> budget; diff --git a/include/dust.h b/include/particle_bin.h similarity index 85% rename from include/dust.h rename to include/particle_bin.h index ad81671cd..fa06503d3 100644 --- a/include/dust.h +++ b/include/particle_bin.h @@ -20,8 +20,8 @@ * along with MicroHH. If not, see . */ -#ifndef DUST_H -#define DUST_H +#ifndef PARTICLE_BIN_H +#define PARTICLE_BIN_H class Master; class Input; @@ -31,11 +31,11 @@ template class Stats; template class Timeloop; template -class Dust +class Particle_bin { public: - Dust(Master&, Grid&, Fields&, Input&); - ~Dust(); + Particle_bin(Master&, Grid&, Fields&, Input&); + ~Particle_bin(); void exec(Stats&); void create(Timeloop&); @@ -46,11 +46,11 @@ class Dust Grid& grid; Fields& fields; - bool sw_dust; + bool sw_particle; TF cfl_max; unsigned long idt_max; // Gravitational settling velocities, negative downward. - std::map w_terminal; + std::map w_particle; }; #endif diff --git a/src/model.cxx b/src/model.cxx index 4e7d5e7c5..de859774f 100644 --- a/src/model.cxx +++ b/src/model.cxx @@ -42,7 +42,7 @@ #include "diff.h" #include "pres.h" #include "force.h" -#include "dust.h" +#include "particle_bin.h" #include "thermo.h" #include "radiation.h" #include "microphys.h" @@ -141,10 +141,11 @@ Model::Model(Master& masterin, int argc, char *argv[]) : decay = std::make_shared>(master, *grid, *fields, *input); limiter = std::make_shared>(master, *grid, *fields, *diff, *input); source = std::make_shared>(master, *grid, *fields, *input); - dust = std::make_shared>(master, *grid, *fields, *input); aerosol = std::make_shared>(master, *grid, *fields, *input); background= std::make_shared>(master, *grid, *fields, *input); + particle_bin = std::make_shared>(master, *grid, *fields, *input); + ib = std::make_shared>(master, *grid, *fields, *input); stats = std::make_shared>(master, *grid, *soil_grid, *background, *fields, *advec, *diff, *input); @@ -267,7 +268,7 @@ void Model::load() buffer->create(*input, *input_nc, *stats); force->create(*input, *input_nc, *stats); source->create(*input, *input_nc); - dust->create(*timeloop); + particle_bin->create(*timeloop); aerosol->create(*input, *input_nc, *stats); background->create(*input, *input_nc, *stats); @@ -417,8 +418,8 @@ void Model::exec() // Add point and line sources of scalars. source->exec(*timeloop); - // Gravitational settling of dust. - dust->exec(*stats); + // Gravitational settling of binned dust types. + particle_bin->exec(*stats); // Apply the large scale forcings. Keep this one always right before the pressure. force->exec(timeloop->get_sub_time_step(), *thermo, *stats); @@ -765,16 +766,16 @@ void Model::set_time_step() // Retrieve the maximum allowed time step per class. timeloop->set_time_step_limit(); - timeloop->set_time_step_limit(advec ->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); - timeloop->set_time_step_limit(diff ->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); - timeloop->set_time_step_limit(thermo ->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); - timeloop->set_time_step_limit(microphys->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); - timeloop->set_time_step_limit(radiation->get_time_limit(timeloop->get_itime())); - timeloop->set_time_step_limit(stats ->get_time_limit(timeloop->get_itime())); - timeloop->set_time_step_limit(cross ->get_time_limit(timeloop->get_itime())); - timeloop->set_time_step_limit(dump ->get_time_limit(timeloop->get_itime())); - timeloop->set_time_step_limit(column ->get_time_limit(timeloop->get_itime())); - timeloop->set_time_step_limit(dust ->get_time_limit()); + timeloop->set_time_step_limit(advec ->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); + timeloop->set_time_step_limit(diff ->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); + timeloop->set_time_step_limit(thermo ->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); + timeloop->set_time_step_limit(microphys ->get_time_limit(timeloop->get_idt(), timeloop->get_dt())); + timeloop->set_time_step_limit(radiation ->get_time_limit(timeloop->get_itime())); + timeloop->set_time_step_limit(stats ->get_time_limit(timeloop->get_itime())); + timeloop->set_time_step_limit(cross ->get_time_limit(timeloop->get_itime())); + timeloop->set_time_step_limit(dump ->get_time_limit(timeloop->get_itime())); + timeloop->set_time_step_limit(column ->get_time_limit(timeloop->get_itime())); + timeloop->set_time_step_limit(particle_bin->get_time_limit()); // Set the time step. timeloop->set_time_step(); diff --git a/src/dust.cu b/src/particle_bin.cu similarity index 85% rename from src/dust.cu rename to src/particle_bin.cu index 1e17f94bd..58c0ec603 100644 --- a/src/dust.cu +++ b/src/particle_bin.cu @@ -22,17 +22,17 @@ #include "grid.h" #include "fields.h" -#include "dust.h" +#include "particle_bin.h" #include "tools.h" namespace { template __global__ - void settle_dust_g( + void settle_particles_g( TF* const __restrict__ st, const TF* const __restrict__ s, const TF* const __restrict__ dzhi, - const TF w_terminal, + const TF w_particles, const int istart, const int iend, const int jstart, const int jend, const int kstart, const int kend, @@ -46,16 +46,16 @@ namespace if (i < iend && j < jend && k < kend) { const int ijk = i + j*jstride + k*kstride; - st[ijk] -= w_terminal * (s[ijk+kstride]-s[ijk])*dzhi[k+1]; + st[ijk] -= w_particles * (s[ijk+kstride]-s[ijk])*dzhi[k+1]; } } } #ifdef USECUDA template -void Dust::exec(Stats& stats) +void Particle_bin::exec(Stats& stats) { - if (!sw_dust) + if (!sw_particle) return; auto& gd = grid.get_grid_data(); @@ -68,8 +68,8 @@ void Dust::exec(Stats& stats) dim3 gridGPU(gridi, gridj, gd.ktot); dim3 blockGPU(blocki, blockj, 1); - for (auto& w : w_terminal) - settle_dust_g<<>>( + for (auto& w : w_particle) + settle_particles_g<<>>( fields.st.at(w.first)->fld_g, fields.sp.at(w.first)->fld_g, gd.dzhi_g, @@ -84,7 +84,7 @@ void Dust::exec(Stats& stats) #endif #ifdef FLOAT_SINGLE -template class Dust; +template class Particle_bin; #else -template class Dust; +template class Particle_bin; #endif diff --git a/src/dust.cxx b/src/particle_bin.cxx similarity index 74% rename from src/dust.cxx rename to src/particle_bin.cxx index 9899df2d0..be1dd2454 100644 --- a/src/dust.cxx +++ b/src/particle_bin.cxx @@ -27,19 +27,20 @@ #include "input.h" #include "grid.h" #include "fields.h" -#include "dust.h" #include "constants.h" #include "timeloop.h" #include "constants.h" +#include "particle_bin.h" + namespace { template - void settle_dust( + void settle_particles( TF* const restrict st, const TF* const restrict s, const TF* const dzhi, - const TF w_terminal, + const TF w_particle, const int istart, const int iend, const int jstart, const int jend, const int kstart, const int kend, @@ -51,46 +52,46 @@ namespace for (int i=istart; i -Dust::Dust(Master& masterin, Grid& gridin, Fields& fieldsin, Input& inputin) : +Particle_bin::Particle_bin(Master& masterin, Grid& gridin, Fields& fieldsin, Input& inputin) : master(masterin), grid(gridin), fields(fieldsin) { - sw_dust = inputin.get_item("dust", "swdust", "", false); + sw_particle = inputin.get_item("particle_bin", "sw_particle", "", false); - if (sw_dust) + if (sw_particle) { - std::vector scalars = inputin.get_list("dust", "dustlist", "", std::vector()); + std::vector scalars = inputin.get_list("particle_bin", "particle_list", "", std::vector()); // Read gravitational settling velocities. for (auto& scalar : scalars) { - w_terminal.emplace(scalar, inputin.get_item("dust", "w_terminal", scalar)); + w_particle.emplace(scalar, inputin.get_item("particle_bin", "w_particle", scalar)); // Raise error if any of the velocities is positive. - if (w_terminal.at(scalar) > 0) + if (w_particle.at(scalar) > 0) throw std::runtime_error("Gravitational settling velocities need to be negative!"); } // Constraint on time stepping. - cfl_max = inputin.get_item("dust", "cfl_max", "", 1.2); + cfl_max = inputin.get_item("particle_bin", "cfl_max", "", 1.2); } } template -Dust::~Dust() +Particle_bin::~Particle_bin() { } template -void Dust::create(Timeloop& timeloop) +void Particle_bin::create(Timeloop& timeloop) { - if (!sw_dust) + if (!sw_particle) return; auto& gd = grid.get_grid_data(); @@ -102,7 +103,7 @@ void Dust::create(Timeloop& timeloop) // Find maximum gravitational settling velocity. TF w_max = -TF(Constants::dbig); - for (auto& w : w_terminal) + for (auto& w : w_particle) w_max = std::max(w_max, std::abs(w.second)); // Calculate maximum time step. @@ -112,9 +113,9 @@ void Dust::create(Timeloop& timeloop) } template -unsigned long Dust::get_time_limit() +unsigned long Particle_bin::get_time_limit() { - if (!sw_dust) + if (!sw_particle) return Constants::ulhuge; return idt_max; @@ -122,15 +123,15 @@ unsigned long Dust::get_time_limit() #ifndef USECUDA template -void Dust::exec(Stats& stats) +void Particle_bin::exec(Stats& stats) { - if (!sw_dust) + if (!sw_particle) return; auto& gd = grid.get_grid_data(); - for (auto& w : w_terminal) - settle_dust( + for (auto& w : w_particle) + settle_particles( fields.st.at(w.first)->fld.data(), fields.sp.at(w.first)->fld.data(), gd.dzhi.data(), @@ -143,7 +144,7 @@ void Dust::exec(Stats& stats) #endif #ifdef FLOAT_SINGLE -template class Dust; +template class Particle_bin; #else -template class Dust; +template class Particle_bin; #endif