Merge pull request microhh#292 from bartvstratum/develop_dust

Renamed `dust` to `particle_bin`.

cases/dust_settling/dust.ini.base → 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

cases/dust_settling/dust_input.py → 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)
 

include/model.h

@@ -46,7 +46,7 @@ template<typename> class Pres;
 template<typename> class Force;
 template<typename> class Aerosol;
 template<typename> class Background;
-template<typename> class Dust;
+template<typename> class Particle_bin;
 template<typename> class Thermo;
 template<typename> class Microphys;
 template<typename> class Radiation;
@@ -107,7 +107,8 @@ class Model
         std::shared_ptr<Decay<TF>> decay;
         std::shared_ptr<Limiter<TF>> limiter;
         std::shared_ptr<Source<TF>> source;
-        std::shared_ptr<Dust<TF>> dust;
+
+        std::shared_ptr<Particle_bin<TF>> particle_bin;
 
         std::shared_ptr<Stats<TF>> stats;
         std::shared_ptr<Budget<TF>> budget;

include/dust.h → include/particle_bin.h

@@ -20,8 +20,8 @@
  * along with MicroHH.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#ifndef DUST_H
-#define DUST_H
+#ifndef PARTICLE_BIN_H
+#define PARTICLE_BIN_H
 
 class Master;
 class Input;
@@ -31,11 +31,11 @@ template<typename> class Stats;
 template<typename> class Timeloop;
 
 template<typename TF>
-class Dust
+class Particle_bin
 {
     public:
-        Dust(Master&, Grid<TF>&, Fields<TF>&, Input&);
-        ~Dust();
+        Particle_bin(Master&, Grid<TF>&, Fields<TF>&, Input&);
+        ~Particle_bin();
 
         void exec(Stats<TF>&);
         void create(Timeloop<TF>&);
@@ -46,11 +46,11 @@ class Dust
         Grid<TF>& grid;
         Fields<TF>& fields;
 
-        bool sw_dust;
+        bool sw_particle;
         TF cfl_max;
         unsigned long idt_max;
 
         // Gravitational settling velocities, negative downward.
-        std::map<std::string, TF> w_terminal;
+        std::map<std::string, TF> w_particle;
 };
 #endif

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<TF>::Model(Master& masterin, int argc, char *argv[]) :
         decay     = std::make_shared<Decay  <TF>>(master, *grid, *fields, *input);
         limiter   = std::make_shared<Limiter<TF>>(master, *grid, *fields, *diff, *input);
         source    = std::make_shared<Source <TF>>(master, *grid, *fields, *input);
-        dust      = std::make_shared<Dust   <TF>>(master, *grid, *fields, *input);
         aerosol   = std::make_shared<Aerosol<TF>>(master, *grid, *fields, *input);
         background= std::make_shared<Background<TF>>(master, *grid, *fields, *input);
 
+        particle_bin = std::make_shared<Particle_bin<TF>>(master, *grid, *fields, *input);
+
         ib        = std::make_shared<Immersed_boundary<TF>>(master, *grid, *fields, *input);
 
         stats     = std::make_shared<Stats <TF>>(master, *grid, *soil_grid, *background, *fields, *advec, *diff, *input);
@@ -267,7 +268,7 @@ void Model<TF>::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<TF>::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<TF>::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();

src/dust.cu → 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<typename TF> __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<typename TF>
-void Dust<TF>::exec(Stats<TF>& stats)
+void Particle_bin<TF>::exec(Stats<TF>& stats)
 {
-    if (!sw_dust)
+    if (!sw_particle)
         return;
 
     auto& gd = grid.get_grid_data();
@@ -68,8 +68,8 @@ void Dust<TF>::exec(Stats<TF>& stats)
     dim3 gridGPU(gridi, gridj, gd.ktot);
     dim3 blockGPU(blocki, blockj, 1);
 
-    for (auto& w : w_terminal)
-        settle_dust_g<TF><<<gridGPU, blockGPU>>>(
+    for (auto& w : w_particle)
+        settle_particles_g<TF><<<gridGPU, blockGPU>>>(
                 fields.st.at(w.first)->fld_g,
                 fields.sp.at(w.first)->fld_g,
                 gd.dzhi_g,
@@ -84,7 +84,7 @@ void Dust<TF>::exec(Stats<TF>& stats)
 #endif
 
 #ifdef FLOAT_SINGLE
-template class Dust<float>;
+template class Particle_bin<float>;
 #else
-template class Dust<double>;
+template class Particle_bin<double>;
 #endif