ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_SuperDropletPC.H
Go to the documentation of this file.
1 #ifndef SUPERDROPLET_PC_H_
2 #define SUPERDROPLET_PC_H_
3 
4 #ifdef ERF_USE_PARTICLES
5 
6 #include <random>
7 #include "ERF_Constants.H"
8 #include "ERF_MaterialProperties.H"
9 #include "ERF_SuperDropletPCDefinitions.H"
10 #include "ERF_SuperDropletPCProcess.H"
11 #include <AMReX_StructOfArrays.H>
12 
13 /*! \brief Ice category for density calculations */
14 enum class IceCategory { Ice, Snow, Graupel, Total };
15 
16 /**
17  * Particle container for the super-droplet method.
18  */
19 class SuperDropletPC : public ERFPC
20 {
21  using MFPtr = std::unique_ptr<amrex::MultiFab>;
22  using BCTypeArr = amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2>;
23 
24  public:
25 
26  /*! \brief Constructor */
27  SuperDropletPC ( amrex::ParGDBBase* a_gdb,
28  const std::vector<Species::Name>& a_species_mat,
29  const std::vector<Species::Name>& a_aerosol_mat,
30  const double a_dt,
31  const std::string& a_name = "super_droplets" )
32  : ERFPC (a_gdb, a_name, ERFPCOptions{ /*read_particle_init=*/false, /*advect_with_gravity=*/true })
33  {
34  define( a_species_mat,
35  a_aerosol_mat,
36  a_gdb->ParticleBoxArray(0),
37  a_gdb->ParticleDistributionMap(0),
38  a_dt );
39  }
40 
41  /*! \brief Constructor */
42  SuperDropletPC ( const amrex::Geometry& a_geom,
43  const amrex::DistributionMapping& a_dmap,
44  const amrex::BoxArray& a_ba,
45  const std::vector<Species::Name>& a_species_mat,
46  const std::vector<Species::Name>& a_aerosol_mat,
47  const double a_dt,
48  const std::string& a_name = "super_droplets" )
49  : ERFPC (a_geom, a_dmap, a_ba, a_name, ERFPCOptions{ /*read_particle_init=*/false, /*advect_with_gravity=*/true })
50  {
51  define(a_species_mat, a_aerosol_mat, a_ba, a_dmap, a_dt);
52  }
53 
54  /*! \brief Destructor */
55  ~SuperDropletPC()
56  {
57  if (m_mass_change_logging) {
58  fclose(m_mass_change_log);
59  }
60  }
61 
62  /*! \brief Set vapour species material */
63  AMREX_FORCE_INLINE
64  virtual void setSpeciesMaterial ( const Species::Name& a_name )
65  {
66  if (a_name == Species::Name::H2O) { m_idx_w = m_species_mat.size(); }
67  if (a_name == Species::Name::ice) { m_idx_i = m_species_mat.size(); }
68  m_species_mat.push_back(std::make_unique<MaterialProperties>(a_name));
69  m_device_props_initialized = false; // Properties have changed, need to reinitialize
70  }
71 
72  /*! \brief Set species material */
73  AMREX_FORCE_INLINE
74  virtual void setSpeciesMaterial ( const std::vector<Species::Name>& a_names )
75  {
76  for (auto& name : a_names) { setSpeciesMaterial(name); }
77  }
78 
79  /*! \brief Get vapour material */
80  AMREX_FORCE_INLINE
81  virtual const MaterialProperties& getSpeciesMaterial(const Species::Name& a_name) const
82  {
83  for (auto& species : m_species_mat) {
84  if (species->m_name == a_name) { return *species; }
85  }
86  amrex::Abort("SuperDropletPC::getSpeciesMaterial() - species not found");
87  return *m_species_mat[0];
88  }
89 
90 
91  /*! \brief Set aerosol material */
92  AMREX_FORCE_INLINE
93  virtual void setAerosolMaterial ( const Species::Name& a_name )
94  {
95  m_aerosol_mat.push_back(std::make_unique<MaterialProperties>(a_name));
96  m_device_props_initialized = false; // Properties have changed, need to reinitialize
97  }
98 
99  /*! \brief Set aerosol material */
100  AMREX_FORCE_INLINE
101  virtual void setAerosolMaterial ( const std::vector<Species::Name>& a_names )
102  {
103  for (auto& name : a_names) { setAerosolMaterial(name); }
104  }
105 
106  /*! \brief Update device properties if material properties change */
107  void updateDeviceProperties();
108 
109  /*! \brief Setup species and aerosol mass pointers from SOA */
110  template<typename SOAType>
111  void setupMassPointers(SOAType& soa, SDPCDefn::SDSpeciesMassArr& sp_mass_ptrs,
112  SDPCDefn::SDAerosolMassArr& ae_mass_ptrs) const
113  {
114  for (int i = 0; i < m_num_species; i++) {
115  sp_mass_ptrs[i] = soa.GetRealData(idx_s(i, m_num_aerosols, m_num_species)).data();
116  }
117  for (int i = 0; i < m_num_aerosols; i++) {
118  ae_mass_ptrs[i] = soa.GetRealData(idx_a(i, m_num_aerosols, m_num_species)).data();
119  }
120  }
121 
122  /*! \brief Update particle radius and total mass from species/aerosol masses */
123  AMREX_GPU_DEVICE AMREX_FORCE_INLINE
124  static void updateParticleAttributes(
125  int particle_idx,
126  amrex::ParticleReal* radius_ptr,
127  amrex::ParticleReal* mass_ptr,
128  int idx_w,
129  amrex::ParticleReal rho_w,
130  int num_sp, int num_ae,
131  const int* sp_sol_arr,
132  const int* ae_sol_arr,
133  const SDPCDefn::SDSpeciesMassArr sp_mass_ptrs,
134  const SDPCDefn::SDAerosolMassArr ae_mass_ptrs,
135  const amrex::ParticleReal* sp_rho_arr,
136  const amrex::ParticleReal* ae_rho_arr)
137  {
138  // Update effective radius
139  radius_ptr[particle_idx] = SD_effective_radius(
140  particle_idx, idx_w, rho_w,
141  num_sp, num_ae,
142  sp_sol_arr, ae_sol_arr,
143  sp_mass_ptrs, ae_mass_ptrs,
144  sp_rho_arr, ae_rho_arr);
145 
146  // Update total mass
147  mass_ptr[particle_idx] = SD_total_mass(
148  particle_idx, num_sp, num_ae,
149  sp_mass_ptrs, ae_mass_ptrs);
150  }
151 
152  /*! \brief Get real-type particle attribute names */
153  [[nodiscard]] virtual amrex::Vector<std::string> varNames () const override;
154 
155  /*! \brief Get Eulerian plot variable names */
156  [[nodiscard]] virtual amrex::Vector<std::string> meshPlotVarNames () const override;
157 
158  /*! \brief Compute mesh variable from particles */
159  virtual void computeMeshVar( const std::string&,
160  amrex::MultiFab&,
161  const amrex::MultiFab&,
162  const int ) const override;
163 
164  /*! \brief Initialize super-droplets in domain at AMR level a_lev */
165  void InitializeParticles (const int a_lev, const double a_t, const MFPtr& a_ptr);
166 
167  using ERFPC::InitializeParticles;
168 
169  /*! \brief Inject super-droplets in domain */
170  virtual void InjectParticles (const double, const MFPtr& a_ptr, const double);
171 
172  /*! \brief add particles for a given set of initialization parameters */
173  void addParticles ( int a_lev, const MFPtr&, const SDInitProperties& );
174 
175  /*! \brief Set initial number of superdroplets per cell as a box with uniform density */
176  void setNumSDBoxDistribution( int a_lev,
177  amrex::iMultiFab&,
178  const int,
179  const MFPtr&,
180  const amrex::RealBox&,
181  const bool );
182 
183  /*! \brief Scatter a box-level count of super-droplets over the injection box */
184  void setNumSDPerBox( int a_lev,
185  amrex::iMultiFab&,
186  const int,
187  const amrex::RealBox&,
188  const unsigned int );
189 
190  /*! \brief Set initial number of superdroplets per cell as a bubble with uniform density */
191  void setNumSDBubbleDistribution( int a_lev,
192  amrex::iMultiFab&,
193  const int,
194  const MFPtr&,
195  const amrex::RealBox&,
196  const bool );
197 
198  /*! \brief Set super-droplets attributes from a given condensate mass density
199  * \param[in] a_mf Multifab containing the initial condensate mass density
200  */
201  virtual void SetAttributes ( amrex::MultiFab& a_mf );
202 
203  /*! \brief Scale the SDM number density with air density
204  * \param[in] a_mf Multifab containing the air density field
205  */
206  virtual void DensityScaling ( const amrex::MultiFab& a_mf );
207 
208  /*! \brief Evolve particles for one time step */
209  virtual void EvolveParticles ( int,
210  double,
211  amrex::Vector<amrex::Vector<amrex::MultiFab>>&,
212  const amrex::Vector<MFPtr>& ) override
213  {
214  amrex::Abort("SuperDropletPC::EvolveParticles() is intentionally disabled.");
215  }
216 
217  /*! \brief Advect particles for one time step
218  * \param[in] a_lev AMR level
219  * \param[in] a_time Current simulation time
220  * \param[in] a_dt Timestep for advection
221  * \param[in] a_flow_vel Array of face-based velocities
222  * \param[in] a_density Density field
223  * \param[in] a_pressure Pressure field
224  * \param[in] a_temperature Temperature field
225  * \param[in] a_z_phys_nd Array of terrain heights
226  * \param[in] a_bctypes Array of boundary condition types
227  * \param[in] a_recycle Flag to enable particle recycling
228  */
229  virtual void AdvectParticles ( int a_lev,
230  double a_time,
231  double a_dt,
232  const amrex::MultiFab* const a_flow_vel,
233  const amrex::MultiFab& a_density,
234  const amrex::MultiFab& a_pressure,
235  const amrex::MultiFab& a_temperature,
236  const amrex::Vector<MFPtr>& a_z_phys_nd,
237  const BCTypeArr& a_bctypes,
238  const bool a_recycle);
239 
240  /*! \brief Condensation and evaporation growth/shrinking for one time step */
241  virtual void MassChange_LV (int,
242  double,
243  const Species::Name&,
244  const amrex::MultiFab&,
245  const amrex::MultiFab&,
246  const amrex::MultiFab&,
247  const amrex::MultiFab&,
248  const amrex::Vector<MFPtr>&,
249  const bool );
250 
251  /*! \brief Freezing/melting for one time step */
252  virtual void MassChange_SL (int,
253  double,
254  const amrex::MultiFab&,
255  const amrex::MultiFab&,
256  const amrex::MultiFab&,
257  const amrex::Vector<MFPtr>& );
258 
259  /*! \brief Deposition and sublimation growth/shrinking for one time step */
260  virtual void MassChange_SV (int,
261  double,
262  const amrex::MultiFab&,
263  const amrex::MultiFab&,
264  const amrex::MultiFab&,
265  const amrex::MultiFab&,
266  const amrex::MultiFab&,
267  const amrex::MultiFab&,
268  const amrex::MultiFab&,
269  const amrex::Vector<MFPtr>& );
270 
271  /*! \brief Coalescence of super-droplets */
272  virtual void Coalescence ( int,
273  double,
274  const amrex::MultiFab&,
275  const amrex::MultiFab&,
276  const amrex::MultiFab&,
277  const amrex::MultiFab&,
278  const amrex::Vector<MFPtr>& );
279 
280  /*! \brief Recycle super-droplets
281  * \param[in] a_lev AMR level
282  * \param[in] a_z_phys_nd Array of terrain heights
283  * \param[in] a_iter Current iteration number
284  * \param[in] a_dt Timestep size
285  * \param[in] a_recycle Flag to enable recycling
286  */
287  virtual void Recycle ( const int a_lev,
288  const amrex::Vector<MFPtr>& a_z_phys_nd,
289  const int a_iter,
290  const double a_dt,
291  const bool a_recycle);
292 
293  /*! \brief Return the number of super-droplets
294  This returns the number of AMReX particles being used in the simulation. */
295  [[nodiscard]] inline virtual amrex::Long NumSuperDroplets ()
296  {
297  return ERFPC::TotalNumberOfParticles();
298  }
299 
300  /*! \brief Return the total number of physical particles */
301  [[nodiscard]] virtual amrex::Real TotalNumberOfParticles ();
302 
303  /*! \brief Return the number of deactivated super-droplets */
304  [[nodiscard]] virtual amrex::Long NumSDDeactivated ();
305 
306  /*! \brief Compute and print diagnostics */
307  virtual void Diagnostics (int, int, double, bool);
308  /*! \brief Compute mass density distribution */
309  virtual void ComputeDistributions ( int, int,
310  amrex::ParticleReal,
311  amrex::ParticleReal );
312 #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS
313  /*! \brief Compute binned size distributions */
314  virtual void ComputeBinnedDistributions ( int, int );
315  /*! \brief Compute cell-wise binned size distributions */
316  virtual void ComputeBinnedDistributionsCell ( int, int, amrex::Real );
317 #endif
318 
319  /*! \brief Compute the SD number density on a mesh */
320  virtual void SDNumberDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int a_comp=0 ) const;
321  /*! \brief Compute the particle number density on a mesh */
322  virtual void numberDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int a_comp=0 ) const;
323  /*! \brief Compute the particle mass density on a mesh */
324  virtual void massDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, const int& a_lev, const int& a_comp=0 ) const override;
325  /*! \brief Compute the mass flux component on a mesh */
326  virtual void massFlux ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int, const int a_comp=0 ) const;
327 
328  /*! \brief Compute the species mass density on a mesh */
329  virtual void speciesMassDensity ( amrex::MultiFab&,
330  const amrex::MultiFab& a_z_phys_nd,
331  int a_lev,
332  int,
333  const int a_comp = 0) const;
334 
335  /*! \brief Compute the cloud/rain mass density on a mesh */
336  virtual void cloudRainDensity ( amrex::MultiFab&,
337  const amrex::MultiFab& a_z_phys_nd,
338  int a_lev,
339  amrex::Real,
340  amrex::Real,
341  const int a_comp = 0) const;
342 
343  /*! \brief Compute ice category mass density on a mesh (ice, snow, graupel, or total) */
344  virtual void iceCategoryDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
345  int a_lev, IceCategory, amrex::Real,
346  const int a_comp = 0) const;
347 
348  /*! \brief Compute the ice mass density on a mesh */
349  virtual void iceDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
350  int a_lev, amrex::Real, const int a_comp = 0) const;
351 
352  /*! \brief Compute the snow mass density on a mesh */
353  virtual void snowDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
354  int a_lev, amrex::Real, const int a_comp = 0) const;
355 
356  /*! \brief Compute the graupel mass density on a mesh */
357  virtual void graupelDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
358  int a_lev, amrex::Real, const int a_comp = 0) const;
359 
360  /*! \brief Compute the frozen water mass density on a mesh */
361  virtual void totalIceDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
362  int a_lev, const int a_comp = 0) const;
363 
364  /*! \brief Compute the species mass flux component on a mesh */
365  virtual void speciesMassFlux ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int, const int, const int a_comp=0 ) const;
366 
367  /*! \brief Compute the aerosol mass density on a mesh */
368  virtual void aerosolMassDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int, const int a_comp=0 ) const;
369  /*! \brief Compute the aerosol mass flux component on a mesh */
370  virtual void aerosolMassFlux ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int, const int, const int a_comp=0 ) const;
371 
372  /*! \brief Compute the particle effective radius on a mesh */
373  virtual void effectiveRadius ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int a_comp=0 ) const;
374 
375  /*! \brief Applies boundary treatment to the particle container */
376  virtual void applyBoundaryTreatment ( int, const amrex::Vector<MFPtr>&, const BCTypeArr&, const bool );
377 
378  /*! \brief Split coarse-level particles in newly-refined cells after AMR
379  * regrid to restore per-cell SD count. Cascades correctly across
380  * multiple new levels: a coarse particle under L2 is split by the
381  * cumulative refRatio(0)*refRatio(1) and tagged for L2. Skips
382  * particles whose multiplicity is below the relevant split factor. */
383  void SplitParticlesForRefinement ( int a_finest_level ) override;
384 
385  /*! \brief Continuously split new entrants on fine levels and merge
386  * excess-tag super-droplets back into native hosts (tag = level + 1).
387  * Per-step departees and regrid-time leftovers from disappeared finer
388  * levels are folded into the same per-level tag-normalizing sweep. */
389  void SplitMergeAtLevelBoundary () override;
390 
391  /*! \brief Regrid-time merge in cells that lost fine-level coverage.
392  * Reduces super-droplet count to `np_bin / merge_factor` so the
393  * coarse-level density is restored after de-refinement. */
394  void MergeParticlesAtDerefinement ( int a_lev,
395  const amrex::BoxArray& a_old_fine_ba,
396  const amrex::IntVect& a_ref_ratio ) override;
397 
398  /*! \brief Whether split/merge across AMR levels is enabled. */
399  bool splitMergeAMR () const { return m_split_merge_amr; }
400 
401  protected:
402 
403  SDCoalescenceKernelType m_coalescence_kernel; /*!< Coalescence kernel */
404  amrex::Real m_ice_agg_eff; /*!< Ice-ice aggregation (collection) efficiency */
405  bool m_include_brownian_coalescence; /*!< Include Brownian coalescence? */
406  SDKernelRelativeVelocityType m_kernel_relative_velocity; /*!< Kernel relative velocity type */
407 
408  SDTerminalVelocityType m_term_vel_type_w; /*!< Terminal velocity model for water droplets */
409  SDTerminalVelocityType m_term_vel_type_i; /*!< Terminal velocity model for ice particles */
410 
411  int m_num_sd_per_cell; /*!< Number of super-droplets per cell */
412  bool m_density_scaling; /*!< Scale initial number density with air density */
413  bool m_nucleate_particles; /*!< nucleate new superdroplets from vapour */
414  bool m_prescribed_advection; /*!< move particles with prescribed vertical velocity */
415  amrex::Real m_prescribed_w = 0.0; /*!< constant prescribed updraft [m/s]; 0 -> sinusoidal pulse */
416  bool m_split_merge_amr = false; /*!< split/merge particles across AMR level boundaries */
417 
418  int m_num_species; /*!< Number of vapour/condensate species */
419  std::vector<std::unique_ptr<MaterialProperties>> m_species_mat; /*!< Vapour/condensate material */
420  int m_num_aerosols; /*!< Number of aerosols */
421  std::vector<std::unique_ptr<MaterialProperties>> m_aerosol_mat; /*!< Aerosol materials */
422 
423  /* Mass change equation Newton solver parameters */
424  amrex::Real m_newton_rtol; /*!< Newton solver - relative tolerance */
425  amrex::Real m_newton_atol; /*!< Newton solver - absolute tolerance */
426  amrex::Real m_newton_stol; /*!< Newton solver - step tolerance */
427  int m_newton_maxits; /*!< Newton solver - maximum iterations */
428 
429  /* Mass change equation (un)convergence logs */
430  bool m_mass_change_logging; /*!< Whether to log unconverged particles */
431  FILE* m_mass_change_log; /*!< File handle for log with unconverged info */
432  std::string m_mass_change_log_fname; /*!< Name of unconverged info log file */
433 
434  amrex::Real m_mass_change_cfl; /*!< CFL for phase change equation */
435  SDMassChangeTIMethod m_mass_change_ti; /*!< time integrator for mass change ODE */
436  long m_num_unconverged_particles; /*!< total number of unconverged particles */
437 
438  /* Ventilation factor for phase change (water species) */
439  bool m_mass_change_ventilation; /*!< include ventilation factor in phase change */
440  amrex::Real m_vent_alpha1; /*!< ventilation factor fit coefficient alpha_1 */
441  amrex::Real m_vent_beta1; /*!< ventilation factor fit exponent beta_1 */
442  amrex::Real m_vent_alpha2; /*!< ventilation factor fit coefficient alpha_2 */
443  amrex::Real m_vent_beta2; /*!< ventilation factor fit exponent beta_2 */
444  amrex::Real m_vent_fcap; /*!< maximum ventilation factor */
445 
446  amrex::IntVect m_coalescence_bin_size; /*!< bin size for coalescence */
447 
448  int m_distribution_grid_size; /*!< Size of 1D grid to compute distributions on */
449 
450  amrex::MultiFab m_mf_buf;
451 
452 #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS
453  amrex::Real m_bindist_rmin; /*!< min radius for binned distribution */
454  amrex::Real m_bindist_rmax; /*!< max radius for binned distribution */
455 
456  /*! Cell-wise mass distribution */
457  amrex::MultiFab m_mass_ln_R_mf;
458  /*! Cell-wise number distribution */
459  amrex::MultiFab m_num_ln_R_mf;
460 #endif
461 
462  /*! maximum coalescence time scale */
463  amrex::Real m_t_coalescence;
464 
465  /*! Initialization parameters objects */
466  int m_num_initializations = 1;
467  /*! vector of initializations */
468  std::vector< std::unique_ptr<SDInitialization> > m_initializations;
469 
470  /*! Injection parameters objects */
471  int m_num_injections = 0;
472  /*! vector of injections */
473  std::vector< std::unique_ptr<SDInjection> > m_injections;
474 
475  /*! sigma_0 parameter for mass distribution calculation */
476  amrex::Real m_sigma0;
477 
478  /*! place the initial SDs randomly within cells */
479  bool m_place_randomly_in_cells;
480 
481  /*! random engine */
482  std::mt19937 m_rndeng;
483 
484  /*! species index of water */
485  int m_idx_w = -1;
486 
487  /*! species index of ice */
488  int m_idx_i = -1;
489 
490  /*! inactive particles threshold */
491  amrex::Real m_deac_threshold;
492 
493  /*! save inactive particles to file */
494  bool m_save_inactive;
495 
496  /*! Persistent device vectors for material properties */
497  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sp_density; // Species densities
498  amrex::Gpu::DeviceVector<int> m_sp_solubility; // Species solubilities
499  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sp_ionization; // Species ionization
500  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sp_mol_weight; // Species molecular weight
501  amrex::Gpu::DeviceVector<int> m_sp_is_INP; // Species ice nucleating particle flags
502 
503  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_ae_density; // Aerosol densities
504  amrex::Gpu::DeviceVector<int> m_ae_solubility; // Aerosol solubilities
505  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_ae_ionization; // Aerosol ionization
506  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_ae_mol_weight; // Aerosol molecular weight
507  amrex::Gpu::DeviceVector<int> m_ae_is_INP; // Aerosol ice nucleating particle flags
508 
509  /*! Flag to track if device properties are initialized */
510  bool m_device_props_initialized = false;
511 
512  /* recycled particle position bounds */
513  amrex::Real m_recyc_xmin;
514  amrex::Real m_recyc_xmax;
515  amrex::Real m_recyc_ymin;
516  amrex::Real m_recyc_ymax;
517  amrex::Real m_recyc_zmin;
518  amrex::Real m_recyc_zmax;
519 
520  /*! Method to initialize device properties */
521  void initializeDeviceProperties();
522 
523  /*! \brief Build ProcessContext with geometry and species info */
524  SDProcess::ProcessContext buildProcessContext(int a_lev) const
525  {
526  SDProcess::ProcessContext ctx;
527  const amrex::Geometry& geom = m_gdb->Geom(a_lev);
528  ctx.plo = geom.ProbLoArray();
529  ctx.phi = geom.ProbHiArray();
530  ctx.dxi = geom.InvCellSizeArray();
531  ctx.dx = geom.CellSizeArray();
532  ctx.domain = geom.Domain();
533  for (int d = 0; d < AMREX_SPACEDIM; d++) {
534  ctx.is_periodic[d] = geom.isPeriodic(d) ? 1 : 0;
535  }
536  const auto cell_size = geom.CellSize();
537  ctx.cell_volume = AMREX_D_TERM(cell_size[0], *cell_size[1], *cell_size[2]);
538  ctx.num_species = m_num_species;
539  ctx.num_aerosols = m_num_aerosols;
540  ctx.idx_water = m_idx_w;
541  ctx.idx_ice = m_idx_i;
542  ctx.rho_water = m_species_mat[m_idx_w]->m_density;
543  if (m_idx_i >= 0) {
544  ctx.rho_ice = m_species_mat[m_idx_i]->m_density;
545  }
546  return ctx;
547  }
548 
549  /*! \brief Setup particle attribute pointers for a tile */
550  template<typename SOAType, typename AOSType>
551  void setupParticlePointers(
552  SOAType& soa,
553  AOSType& aos,
554  SDProcess::ParticlePointers& ptrs) const
555  {
556  using namespace SDPCDefn;
557 
558  constexpr int rtoff_i = SuperDropletsIntIdx::ncomps;
559  constexpr int rtoff_r = SuperDropletsRealIdx::ncomps;
560 
561  ptrs.num_particles = aos.numParticles();
562  ptrs.mass_ptr = soa.GetRealData(SuperDropletsRealIdx::mass).data();
563  ptrs.radius_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::radius).data();
564  ptrs.active_ptr = soa.GetIntData(rtoff_i + SuperDropletsIntIdxSoA_RT::active).data();
565  ptrs.v_ptr[0] = soa.GetRealData(SuperDropletsRealIdx::vx).data();
566  ptrs.v_ptr[1] = soa.GetRealData(SuperDropletsRealIdx::vy).data();
567  ptrs.v_ptr[2] = soa.GetRealData(SuperDropletsRealIdx::vz).data();
568  ptrs.vterm_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::term_vel).data();
569  ptrs.mult_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::multiplicity).data();
570  // Always set ice-related pointers: SoA always has these comps; water-only coalescence uses Tfz_ptr
571  ptrs.Tfz_ptr = soa.GetRealData(idx_ice_Tfz(m_num_aerosols, m_num_species)).data();
572  ptrs.a_ptr = soa.GetRealData(idx_ice_a(m_num_aerosols, m_num_species)).data();
573  ptrs.c_ptr = soa.GetRealData(idx_ice_c(m_num_aerosols, m_num_species)).data();
574  ptrs.mrime_ptr = soa.GetRealData(idx_ice_mrime(m_num_aerosols, m_num_species)).data();
575  ptrs.nmono_ptr = soa.GetRealData(idx_ice_nmono(m_num_aerosols, m_num_species)).data();
576  setupMassPointers(soa, ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs);
577  if (!m_device_props_initialized) {
578  const_cast<SuperDropletPC*>(this)->initializeDeviceProperties();
579  }
580  ptrs.sp_rho_arr = m_sp_density.data();
581  ptrs.sp_sol_arr = m_sp_solubility.data();
582  ptrs.sp_ion_arr = m_sp_ionization.data();
583  ptrs.sp_mw_arr = m_sp_mol_weight.data();
584  ptrs.sp_INP_arr = m_sp_is_INP.data();
585  ptrs.ae_rho_arr = m_ae_density.data();
586  ptrs.ae_sol_arr = m_ae_solubility.data();
587  ptrs.ae_ion_arr = m_ae_ionization.data();
588  ptrs.ae_mw_arr = m_ae_mol_weight.data();
589  ptrs.ae_INP_arr = m_ae_is_INP.data();
590  }
591 
592  /*! \brief Common body for forEachParticleTile iterations */
593  template<typename TileFunc>
594  void forEachParticleTileBody(
595  ParIterType& pti, int a_lev,
596  const SDProcess::ProcessContext& ctx,
597  TileFunc&& func)
598  {
599  int grid = pti.index();
600  auto& ptile = ParticlesAt(a_lev, pti);
601  auto& aos = ptile.GetArrayOfStructs();
602  auto& soa = ptile.GetStructOfArrays();
603  if (aos.numParticles() == 0) { return; }
604  auto* p_pbox = aos().data();
605  SDProcess::ParticlePointers ptrs;
606  setupParticlePointers(soa, aos, ptrs);
607  func(pti, grid, p_pbox, ptrs, ctx);
608  }
609 
610  /*! \brief Iterate over particle tiles with common setup (parallel) */
611  template<typename TileFunc>
612  AMREX_FORCE_INLINE
613  void forEachParticleTile(
614  int a_lev,
615  const SDProcess::ProcessContext& ctx,
616  TileFunc&& func)
617  {
618 #ifdef AMREX_USE_OMP
619 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
620 #endif
621  for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) {
622  forEachParticleTileBody(pti, a_lev, ctx, std::forward<TileFunc>(func));
623  }
624  }
625 
626  /*! \brief Iterate over particle tiles WITHOUT OpenMP (serial)
627  * Use when per-tile operations are not thread-safe (e.g., DenseBins).
628  */
629  template<typename TileFunc>
630  AMREX_FORCE_INLINE
631  void forEachParticleTileSerial(int a_lev, const SDProcess::ProcessContext& ctx, TileFunc&& func)
632  {
633  for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) {
634  forEachParticleTileBody(pti, a_lev, ctx, std::forward<TileFunc>(func));
635  }
636  }
637 
638  /*! \brief Lightweight iteration (no ProcessContext) */
639  template<typename TileFunc>
640  void forEachParticleTile(int a_lev, TileFunc&& func)
641  {
642 #ifdef AMREX_USE_OMP
643 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
644 #endif
645  for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) {
646  int grid = pti.index();
647  auto& ptile = ParticlesAt(a_lev, pti);
648  auto& aos = ptile.GetArrayOfStructs();
649  auto& soa = ptile.GetStructOfArrays();
650  const int num_particles = aos.numParticles();
651  if (num_particles == 0) { continue; }
652  auto* p_pbox = aos().data();
653  SDProcess::ParticlePointers ptrs;
654  setupParticlePointers(soa, aos, ptrs);
655  func(pti, grid, p_pbox, ptrs, num_particles);
656  }
657  }
658 
659  /*! \brief read inputs from file */
660  virtual void readInputs () override
661  {
662  amrex::Abort("SuperDropletPC::readInputs(): Do not use this interface.");
663  }
664 
665  /*! \brief read inputs from file */
666  virtual void readInputs (const double);
667 
668  /*! \brief Particle initialization - null (no particles are initialized) */
669  void initializeParticlesNull (const MFPtr&) { }
670 
671  private:
672 
673  /*! \brief define super-droplets */
674  void define ( const std::vector<Species::Name>&,
675  const std::vector<Species::Name>&,
676  const amrex::BoxArray&,
677  const amrex::DistributionMapping&,
678  const double );
679 
680  /*! \brief add super-droplet method-specific particle attributes */
681  void add_superdroplet_attributes();
682 
683 };
684 
685 #endif
686 #endif
ERF_Constants.H
m_num_species
int m_num_species
Definition: ERF_InitCustomPert_MultiSpeciesBubble.H:28
ERF_MaterialProperties.H
name
std::string name
Definition: ERF_Plotfile2DCatalog.cpp:101
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
ERF_SuperDropletPCDefinitions.H
ERF_SuperDropletPCProcess.H
Common data structures for SuperDroplet physical processes.
SDInitProperties
Super-droplets initial properties.
Definition: ERF_SDInitialization.H:220
MaterialProperties
Definition: ERF_MaterialProperties.H:187