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_NumericalConstants.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  bool composite_init = false );
175 
176  /*! \brief Set initial number of superdroplets per cell as a box with uniform density */
177  void setNumSDBoxDistribution ( int a_lev,
178  amrex::iMultiFab&,
179  const int,
180  const MFPtr&,
181  const amrex::RealBox&,
182  const bool );
183 
184  /*! \brief Scatter a box-level count of super-droplets over the injection box */
185  void setNumSDPerBox ( int a_lev,
186  amrex::iMultiFab&,
187  const int,
188  const amrex::RealBox&,
189  const unsigned int );
190 
191  /*! \brief Set initial number of superdroplets per cell as a bubble with uniform density */
192  void setNumSDBubbleDistribution ( int a_lev,
193  amrex::iMultiFab&,
194  const int,
195  const MFPtr&,
196  const amrex::RealBox&,
197  const bool );
198 
199  /*! \brief Set super-droplets attributes from a given condensate mass density
200  * \param[in] a_mf Multifab containing the initial condensate mass density
201  */
202  virtual void SetAttributes ( amrex::MultiFab& a_mf, int a_lev = 0 );
203 
204  /*! \brief Scale the SDM number density with air density
205  * \param[in] a_mf Multifab containing the air density field
206  */
207  virtual void DensityScaling ( const amrex::MultiFab& a_mf, int a_lev = 0 );
208 
209  /*! \brief Evolve particles for one time step */
210  virtual void EvolveParticles ( int,
211  double,
212  amrex::Vector<amrex::Vector<amrex::MultiFab>>&,
213  const amrex::Vector<MFPtr>& ) override
214  {
215  amrex::Abort("SuperDropletPC::EvolveParticles() is intentionally disabled.");
216  }
217 
218  /*! \brief Advect particles for one time step
219  * \param[in] a_lev AMR level
220  * \param[in] a_time Current simulation time
221  * \param[in] a_dt Timestep for advection
222  * \param[in] a_flow_vel Array of face-based velocities
223  * \param[in] a_density Density field
224  * \param[in] a_pressure Pressure field
225  * \param[in] a_temperature Temperature field
226  * \param[in] a_z_phys_nd Array of terrain heights
227  * \param[in] a_bctypes Array of boundary condition types
228  * \param[in] a_recycle Flag to enable particle recycling
229  */
230  virtual void AdvectParticles ( int a_lev,
231  double a_time,
232  double a_dt,
233  const amrex::MultiFab* const a_flow_vel,
234  const amrex::MultiFab& a_density,
235  const amrex::MultiFab& a_pressure,
236  const amrex::MultiFab& a_temperature,
237  const amrex::Vector<MFPtr>& a_z_phys_nd,
238  const BCTypeArr& a_bctypes,
239  const bool a_recycle);
240 
241  /*! \brief Condensation and evaporation growth/shrinking for one time step */
242  virtual void MassChange_LV (int,
243  double,
244  const Species::Name&,
245  const amrex::MultiFab&,
246  const amrex::MultiFab&,
247  const amrex::MultiFab&,
248  const amrex::MultiFab&,
249  const amrex::Vector<MFPtr>&,
250  const bool );
251 
252  /*! \brief Freezing/melting for one time step */
253  virtual void MassChange_SL (int,
254  double,
255  const amrex::MultiFab&,
256  const amrex::MultiFab&,
257  const amrex::MultiFab&,
258  const amrex::Vector<MFPtr>& );
259 
260  /*! \brief Deposition and sublimation growth/shrinking for one time step */
261  virtual void MassChange_SV (int,
262  double,
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::MultiFab&,
270  const amrex::Vector<MFPtr>& );
271 
272  /*! \brief Coalescence of super-droplets */
273  virtual void Coalescence ( int,
274  double,
275  const amrex::MultiFab&,
276  const amrex::MultiFab&,
277  const amrex::MultiFab&,
278  const amrex::MultiFab&,
279  const amrex::Vector<MFPtr>& );
280 
281  /*! \brief Recycle super-droplets
282  * \param[in] a_lev AMR level
283  * \param[in] a_z_phys_nd Array of terrain heights
284  * \param[in] a_iter Current iteration number
285  * \param[in] a_dt Timestep size
286  * \param[in] a_recycle Flag to enable recycling
287  */
288  virtual void Recycle ( const int a_lev,
289  const amrex::Vector<MFPtr>& a_z_phys_nd,
290  const int a_iter,
291  const double a_dt,
292  const bool a_recycle);
293 
294  /*! \brief Set the vertical coordinate of recycled super-droplets on their owning grids
295  * \param[in] a_z_phys_nd Array of terrain heights
296  */
297  void placeRecycledParticles ( const amrex::Vector<MFPtr>& a_z_phys_nd );
298 
299  /*! \brief Return the number of super-droplets
300  This returns the number of AMReX particles being used in the simulation. */
301  [[nodiscard]] inline virtual amrex::Long NumSuperDroplets ()
302  {
303  return ERFPC::TotalNumberOfParticles();
304  }
305 
306  /*! \brief Return the total number of physical particles */
307  [[nodiscard]] virtual amrex::Real TotalNumberOfParticles ();
308 
309  /*! \brief Return the number of deactivated super-droplets */
310  [[nodiscard]] virtual amrex::Long NumSDDeactivated ();
311 
312  /*! \brief Compute and print diagnostics */
313  virtual void Diagnostics (int, int, double, bool);
314  /*! \brief Compute mass density distribution */
315  virtual void ComputeDistributions ( int, int,
316  amrex::ParticleReal,
317  amrex::ParticleReal );
318 #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS
319  /*! \brief Compute binned size distributions */
320  virtual void ComputeBinnedDistributions ( int, int );
321  /*! \brief Compute cell-wise binned size distributions */
322  virtual void ComputeBinnedDistributionsCell ( int, int, amrex::Real );
323 #endif
324 
325  /*! \brief Compute the SD number density on a mesh */
326  virtual void SDNumberDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int a_comp=0 ) const;
327  /*! \brief Compute the particle number density on a mesh */
328  virtual void numberDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int a_comp=0 ) const;
329  /*! \brief Compute the particle mass density on a mesh */
330  virtual void massDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, const int& a_lev, const int& a_comp=0 ) const override;
331  /*! \brief Compute the mass flux component on a mesh */
332  virtual void massFlux ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int, const int a_comp=0 ) const;
333 
334  /*! \brief Compute the species mass density on a mesh */
335  virtual void speciesMassDensity ( amrex::MultiFab&,
336  const amrex::MultiFab& a_z_phys_nd,
337  int a_lev,
338  int,
339  const int a_comp = 0) const;
340 
341  /*! \brief Compute the cloud/rain mass density on a mesh */
342  virtual void cloudRainDensity ( amrex::MultiFab&,
343  const amrex::MultiFab& a_z_phys_nd,
344  int a_lev,
345  amrex::Real,
346  amrex::Real,
347  const int a_comp = 0) const;
348 
349  /*! \brief Compute ice category mass density on a mesh (ice, snow, graupel, or total) */
350  virtual void iceCategoryDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
351  int a_lev, IceCategory, amrex::Real,
352  const int a_comp = 0) const;
353 
354  /*! \brief Compute the ice mass density on a mesh */
355  virtual void iceDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
356  int a_lev, amrex::Real, const int a_comp = 0) const;
357 
358  /*! \brief Compute the snow mass density on a mesh */
359  virtual void snowDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
360  int a_lev, amrex::Real, const int a_comp = 0) const;
361 
362  /*! \brief Compute the graupel mass density on a mesh */
363  virtual void graupelDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
364  int a_lev, amrex::Real, const int a_comp = 0) const;
365 
366  /*! \brief Compute the frozen water mass density on a mesh */
367  virtual void totalIceDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd,
368  int a_lev, const int a_comp = 0) const;
369 
370  /*! \brief Compute the species mass flux component on a mesh */
371  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;
372 
373  /*! \brief Compute the aerosol mass density on a mesh */
374  virtual void aerosolMassDensity ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int, const int a_comp=0 ) const;
375  /*! \brief Compute the aerosol mass flux component on a mesh */
376  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;
377 
378  /*! \brief Compute the particle effective radius on a mesh */
379  virtual void effectiveRadius ( amrex::MultiFab&, const amrex::MultiFab& a_z_phys_nd, int a_lev, const int a_comp=0 ) const;
380 
381  /*! \brief Applies boundary treatment to the particle container */
382  virtual void applyBoundaryTreatment ( int, const amrex::Vector<MFPtr>&, const BCTypeArr&, const bool );
383 
384  /*! \brief Split coarse-level particles in newly-refined cells after AMR
385  * regrid to restore per-cell SD count. Cascades correctly across
386  * multiple new levels: a coarse particle under L2 is split by the
387  * cumulative refRatio(0)*refRatio(1) and tagged for L2. Skips
388  * particles whose multiplicity is below the relevant split factor. */
389  void SplitParticlesForRefinement ( int a_finest_level ) override;
390 
391  /*! \brief Continuously split new entrants on fine levels and merge
392  * excess-tag super-droplets back into native hosts (tag = level + 1).
393  * Per-step departees and regrid-time leftovers from disappeared finer
394  * levels are folded into the same per-level tag-normalizing sweep. */
395  void SplitMergeAtLevelBoundary () override;
396 
397  /*! \brief Regrid-time merge in cells that lost fine-level coverage.
398  * Reduces super-droplet count to `np_bin / merge_factor` so the
399  * coarse-level density is restored after de-refinement. */
400  void MergeParticlesAtDerefinement ( int a_lev,
401  const amrex::BoxArray& a_old_fine_ba,
402  const amrex::IntVect& a_ref_ratio ) override;
403 
404  /*! \brief Whether split/merge across AMR levels is enabled. */
405  bool splitMergeAMR () const { return m_split_merge_amr; }
406 
407  protected:
408 
409  SDCoalescenceKernelType m_coalescence_kernel; /*!< Coalescence kernel */
410  amrex::Real m_ice_agg_eff; /*!< Ice-ice aggregation (collection) efficiency */
411  bool m_include_brownian_coalescence; /*!< Include Brownian coalescence? */
412  SDKernelRelativeVelocityType m_kernel_relative_velocity; /*!< Kernel relative velocity type */
413 
414  SDTerminalVelocityType m_term_vel_type_w; /*!< Terminal velocity model for water droplets */
415  SDTerminalVelocityType m_term_vel_type_i; /*!< Terminal velocity model for ice particles */
416 
417  int m_num_sd_per_cell; /*!< Number of super-droplets per cell */
418  bool m_density_scaling; /*!< Scale initial number density with air density */
419  bool m_nucleate_particles; /*!< nucleate new superdroplets from vapour */
420  bool m_prescribed_advection; /*!< move particles with prescribed vertical velocity */
421  amrex::Real m_prescribed_w = 0.0; /*!< constant prescribed updraft [m/s]; 0 -> sinusoidal pulse */
422  bool m_split_merge_amr = false; /*!< split/merge particles across AMR level boundaries */
423 
424  int m_num_species; /*!< Number of vapour/condensate species */
425  std::vector<std::unique_ptr<MaterialProperties>> m_species_mat; /*!< Vapour/condensate material */
426  int m_num_aerosols; /*!< Number of aerosols */
427  std::vector<std::unique_ptr<MaterialProperties>> m_aerosol_mat; /*!< Aerosol materials */
428 
429  /* Mass change equation Newton solver parameters */
430  amrex::Real m_newton_rtol; /*!< Newton solver - relative tolerance */
431  amrex::Real m_newton_atol; /*!< Newton solver - absolute tolerance */
432  amrex::Real m_newton_stol; /*!< Newton solver - step tolerance */
433  int m_newton_maxits; /*!< Newton solver - maximum iterations */
434 
435  /* Mass change equation (un)convergence logs */
436  bool m_mass_change_logging; /*!< Whether to log unconverged particles */
437  FILE* m_mass_change_log; /*!< File handle for log with unconverged info */
438  std::string m_mass_change_log_fname; /*!< Name of unconverged info log file */
439 
440  amrex::Real m_mass_change_cfl; /*!< CFL for phase change equation */
441  SDMassChangeTIMethod m_mass_change_ti; /*!< time integrator for mass change ODE */
442  long m_num_unconverged_particles; /*!< total number of unconverged particles */
443 
444  /* Ventilation factor for phase change (water species) */
445  bool m_mass_change_ventilation; /*!< include ventilation factor in phase change */
446  amrex::Real m_vent_alpha1; /*!< ventilation factor fit coefficient alpha_1 */
447  amrex::Real m_vent_beta1; /*!< ventilation factor fit exponent beta_1 */
448  amrex::Real m_vent_alpha2; /*!< ventilation factor fit coefficient alpha_2 */
449  amrex::Real m_vent_beta2; /*!< ventilation factor fit exponent beta_2 */
450  amrex::Real m_vent_fcap; /*!< maximum ventilation factor */
451 
452  amrex::IntVect m_coalescence_bin_size; /*!< bin size for coalescence */
453 
454  int m_distribution_grid_size; /*!< Size of 1D grid to compute distributions on */
455 
456  amrex::MultiFab m_mf_buf;
457 
458 #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS
459  amrex::Real m_bindist_rmin; /*!< min radius for binned distribution */
460  amrex::Real m_bindist_rmax; /*!< max radius for binned distribution */
461 
462  /*! Cell-wise mass distribution */
463  amrex::MultiFab m_mass_ln_R_mf;
464  /*! Cell-wise number distribution */
465  amrex::MultiFab m_num_ln_R_mf;
466 #endif
467 
468  /*! maximum coalescence time scale */
469  amrex::Real m_t_coalescence;
470 
471  /*! Initialization parameters objects */
472  int m_num_initializations = 1;
473  /*! vector of initializations */
474  std::vector< std::unique_ptr<SDInitialization> > m_initializations;
475 
476  /*! Injection parameters objects */
477  int m_num_injections = 0;
478  /*! vector of injections */
479  std::vector< std::unique_ptr<SDInjection> > m_injections;
480 
481  /*! sigma_0 parameter for mass distribution calculation */
482  amrex::Real m_sigma0;
483 
484  /*! place the initial SDs randomly within cells */
485  bool m_place_randomly_in_cells;
486 
487  /*! random engine */
488  std::mt19937 m_rndeng;
489 
490  /*! species index of water */
491  int m_idx_w = -1;
492 
493  /*! species index of ice */
494  int m_idx_i = -1;
495 
496  /*! inactive particles threshold */
497  amrex::Real m_deac_threshold;
498 
499  /*! save inactive particles to file */
500  bool m_save_inactive;
501 
502  /*! Persistent device vectors for material properties */
503  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sp_density; // Species densities
504  amrex::Gpu::DeviceVector<int> m_sp_solubility; // Species solubilities
505  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sp_ionization; // Species ionization
506  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sp_mol_weight; // Species molecular weight
507  amrex::Gpu::DeviceVector<int> m_sp_is_INP; // Species ice nucleating particle flags
508 
509  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_ae_density; // Aerosol densities
510  amrex::Gpu::DeviceVector<int> m_ae_solubility; // Aerosol solubilities
511  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_ae_ionization; // Aerosol ionization
512  amrex::Gpu::DeviceVector<amrex::ParticleReal> m_ae_mol_weight; // Aerosol molecular weight
513  amrex::Gpu::DeviceVector<int> m_ae_is_INP; // Aerosol ice nucleating particle flags
514 
515  /*! Flag to track if device properties are initialized */
516  bool m_device_props_initialized = false;
517 
518  /* recycled particle position bounds */
519  amrex::Real m_recyc_xmin;
520  amrex::Real m_recyc_xmax;
521  amrex::Real m_recyc_ymin;
522  amrex::Real m_recyc_ymax;
523  amrex::Real m_recyc_zmin;
524  amrex::Real m_recyc_zmax;
525 
526  /*! Method to initialize device properties */
527  void initializeDeviceProperties ();
528 
529  /*! \brief Build ProcessContext with geometry and species info */
530  SDProcess::ProcessContext buildProcessContext (int a_lev) const
531  {
532  SDProcess::ProcessContext ctx;
533  const amrex::Geometry& geom = m_gdb->Geom(a_lev);
534  ctx.plo = geom.ProbLoArray();
535  ctx.phi = geom.ProbHiArray();
536  ctx.dxi = geom.InvCellSizeArray();
537  ctx.dx = geom.CellSizeArray();
538  ctx.domain = geom.Domain();
539  for (int d = 0; d < AMREX_SPACEDIM; d++) {
540  ctx.is_periodic[d] = geom.isPeriodic(d) ? 1 : 0;
541  }
542  const auto cell_size = geom.CellSize();
543  ctx.cell_volume = AMREX_D_TERM(cell_size[0], *cell_size[1], *cell_size[2]);
544  ctx.num_species = m_num_species;
545  ctx.num_aerosols = m_num_aerosols;
546  ctx.idx_water = m_idx_w;
547  ctx.idx_ice = m_idx_i;
548  ctx.rho_water = m_species_mat[m_idx_w]->m_density;
549  if (m_idx_i >= 0) {
550  ctx.rho_ice = m_species_mat[m_idx_i]->m_density;
551  }
552  return ctx;
553  }
554 
555  /*! \brief Setup particle attribute pointers for a tile */
556  template<typename SOAType, typename AOSType>
557  void setupParticlePointers (
558  SOAType& soa,
559  AOSType& aos,
560  SDProcess::ParticlePointers& ptrs) const
561  {
562  using namespace SDPCDefn;
563 
564  constexpr int rtoff_i = SuperDropletsIntIdx::ncomps;
565  constexpr int rtoff_r = SuperDropletsRealIdx::ncomps;
566 
567  ptrs.num_particles = aos.numParticles();
568  ptrs.mass_ptr = soa.GetRealData(SuperDropletsRealIdx::mass).data();
569  ptrs.radius_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::radius).data();
570  ptrs.active_ptr = soa.GetIntData(rtoff_i + SuperDropletsIntIdxSoA_RT::active).data();
571  ptrs.v_ptr[0] = soa.GetRealData(SuperDropletsRealIdx::vx).data();
572  ptrs.v_ptr[1] = soa.GetRealData(SuperDropletsRealIdx::vy).data();
573  ptrs.v_ptr[2] = soa.GetRealData(SuperDropletsRealIdx::vz).data();
574  ptrs.vterm_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::term_vel).data();
575  ptrs.mult_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::multiplicity).data();
576  // Always set ice-related pointers: SoA always has these comps; water-only coalescence uses Tfz_ptr
577  ptrs.Tfz_ptr = soa.GetRealData(idx_ice_Tfz(m_num_aerosols, m_num_species)).data();
578  ptrs.a_ptr = soa.GetRealData(idx_ice_a(m_num_aerosols, m_num_species)).data();
579  ptrs.c_ptr = soa.GetRealData(idx_ice_c(m_num_aerosols, m_num_species)).data();
580  ptrs.mrime_ptr = soa.GetRealData(idx_ice_mrime(m_num_aerosols, m_num_species)).data();
581  ptrs.nmono_ptr = soa.GetRealData(idx_ice_nmono(m_num_aerosols, m_num_species)).data();
582  setupMassPointers(soa, ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs);
583  if (!m_device_props_initialized) {
584  const_cast<SuperDropletPC*>(this)->initializeDeviceProperties();
585  }
586  ptrs.sp_rho_arr = m_sp_density.data();
587  ptrs.sp_sol_arr = m_sp_solubility.data();
588  ptrs.sp_ion_arr = m_sp_ionization.data();
589  ptrs.sp_mw_arr = m_sp_mol_weight.data();
590  ptrs.sp_INP_arr = m_sp_is_INP.data();
591  ptrs.ae_rho_arr = m_ae_density.data();
592  ptrs.ae_sol_arr = m_ae_solubility.data();
593  ptrs.ae_ion_arr = m_ae_ionization.data();
594  ptrs.ae_mw_arr = m_ae_mol_weight.data();
595  ptrs.ae_INP_arr = m_ae_is_INP.data();
596  }
597 
598  /*! \brief Common body for forEachParticleTile iterations */
599  template<typename TileFunc>
600  void forEachParticleTileBody (
601  ParIterType& pti, int a_lev,
602  const SDProcess::ProcessContext& ctx,
603  TileFunc&& func)
604  {
605  int grid = pti.index();
606  auto& ptile = ParticlesAt(a_lev, pti);
607  auto& aos = ptile.GetArrayOfStructs();
608  auto& soa = ptile.GetStructOfArrays();
609  if (aos.numParticles() == 0) { return; }
610  auto* p_pbox = aos().data();
611  SDProcess::ParticlePointers ptrs;
612  setupParticlePointers(soa, aos, ptrs);
613  func(pti, grid, p_pbox, ptrs, ctx);
614  }
615 
616  /*! \brief Iterate over particle tiles with common setup (parallel) */
617  template<typename TileFunc>
618  AMREX_FORCE_INLINE
619  void forEachParticleTile (
620  int a_lev,
621  const SDProcess::ProcessContext& ctx,
622  TileFunc&& func)
623  {
624 #ifdef AMREX_USE_OMP
625 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
626 #endif
627  for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) {
628  forEachParticleTileBody(pti, a_lev, ctx, std::forward<TileFunc>(func));
629  }
630  }
631 
632  /*! \brief Iterate over particle tiles WITHOUT OpenMP (serial)
633  * Use when per-tile operations are not thread-safe (e.g., DenseBins).
634  */
635  template<typename TileFunc>
636  AMREX_FORCE_INLINE
637  void forEachParticleTileSerial (int a_lev, const SDProcess::ProcessContext& ctx, TileFunc&& func)
638  {
639  for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) {
640  forEachParticleTileBody(pti, a_lev, ctx, std::forward<TileFunc>(func));
641  }
642  }
643 
644  /*! \brief Lightweight iteration (no ProcessContext) */
645  template<typename TileFunc>
646  void forEachParticleTile (int a_lev, TileFunc&& func)
647  {
648 #ifdef AMREX_USE_OMP
649 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
650 #endif
651  for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) {
652  int grid = pti.index();
653  auto& ptile = ParticlesAt(a_lev, pti);
654  auto& aos = ptile.GetArrayOfStructs();
655  auto& soa = ptile.GetStructOfArrays();
656  const int num_particles = aos.numParticles();
657  if (num_particles == 0) { continue; }
658  auto* p_pbox = aos().data();
659  SDProcess::ParticlePointers ptrs;
660  setupParticlePointers(soa, aos, ptrs);
661  func(pti, grid, p_pbox, ptrs, num_particles);
662  }
663  }
664 
665  /*! \brief read inputs from file */
666  virtual void readInputs () override
667  {
668  amrex::Abort("SuperDropletPC::readInputs(): Do not use this interface.");
669  }
670 
671  /*! \brief read inputs from file */
672  virtual void readInputs (const double);
673 
674  /*! \brief Particle initialization - null (no particles are initialized) */
675  void initializeParticlesNull (const MFPtr&) { }
676 
677  private:
678 
679  /*! \brief define super-droplets */
680  void define ( const std::vector<Species::Name>&,
681  const std::vector<Species::Name>&,
682  const amrex::BoxArray&,
683  const amrex::DistributionMapping&,
684  const double );
685 
686  /*! \brief add super-droplet method-specific particle attributes */
687  void add_superdroplet_attributes ();
688 
689 };
690 
691 #endif
692 #endif
m_num_species
int m_num_species
Definition: ERF_InitCustomPert_MultiSpeciesBubble.H:32
ERF_MaterialProperties.H
ERF_NumericalConstants.H
Dimensionless numeric literals and pure mathematical constants.
name
std::string name
Definition: ERF_Plotfile2DCatalog.cpp:116
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:223
erf_sbm::PopulationPhase::Ice
@ Ice
MaterialProperties
Definition: ERF_MaterialProperties.H:178