Peano 4
Loading...
Searching...
No Matches
SPHLeapfrogFixedSearchRadius.py
Go to the documentation of this file.
1# This file is part of the SWIFT2 project. For conditions of distribution and
2# use, please see the copyright notice at www.peano-framework.org
3from swift2.particle.SPHParticle import SPHParticle
4from swift2.particle.AlgorithmStep import AlgorithmStep
5from swift2.particle.AlgorithmStepLibrary import get_algorithm_step_dict
6
7import peano4
8import dastgen2
10
11from abc import ABC
12import copy
13
14
15from enum import Enum
16
17
19 """!
20
21 SPH particle with fixed search radius
22
23 Same as LeapfrogFixedSearchRadius but augmented for full SPH.
24 Implements the Minimal SPH scheme.
25
26 Simple particle with a fixed search radius h which moves according
27 to leapfrog KDK scheme.
28 By default, it uses global time stepping,
29 i.e. the combination of maximum velocity and minimal mesh size determines
30 the time step size of the subsequent time step. Besides the default
31 variables x and h, the particle has the following properties:
32
33 - a vector a which is the acceleration;
34 - a vector v which is the velocity.
35
36 You can add further properties via
37
38 myparticle.data.add_attribute( peano4.dastgen2.Peano4DoubleArray("myFancyArray","Dimensions") )
39
40 in your code. Or you can create a subclass which adds additional fields
41 after it has called the baseline constructor.
42
43 You will need to add further properties for any SPH project.
44
45 ## Force calculation
46
47 particle_particle_interaction_over_particle_sets_kernel is a C++ string which defines a force between two
48 particles. It has access to three important objects:
49
50 - localParticles is a container (list) over pointers to local particles
51 - activeParticles is a container (list) over pointers to active particles
52 - marker is a cell marker which identifies the current cell.
53
54 Please consult the guidebook for a definition of local and active particles
55 but take into account that the local particles always are a subset of the
56 active particles.
57
58 Besides the actual particle-to-particle calculation, i.e. a force
59 calculation, users can also provide kernels that kick in when you touch
60 particles for the first time before you actually compute any
61 particle-particle interaction, and there is a plug-in point what you do
62 just once the particle-particle interaction has terminated. The latter
63 point is reached before we do the actual time stepping. In both plug-in
64 points, you have a vertex marker which gives you the position of the
65 vertex to which a particular is associated, and you have the localParticles.
66 This is a vector of pointers in this particular case.
67
68 The force calculation does exist in two variants: There's the default one
69 and then we have one that is aggressively optimised. The latter one requires
70 us to have properly sorted particle-mesh associations.
71
72 ## Precision
73
74 By default, all particle attributes are modelled as doubles or arrays of
75 doubles. This might be overly precise. You can alter the precision through
76 our LLVM compiler modifications if you invoke the following commands on a
77 particle species:
78
79 ~~~~~~~~~~~~~~~~~~~~~~~~~~
80 particle.mass.valid_mantissa_bits = args.mantissa_size
81 particle.velocity.valid_mantissa_bits = args.mantissa_size
82 particle.acceleration.valid_mantissa_bits = args.mantissa_size
83 particle.density.valid_mantissa_bits = args.mantissa_size
84 particle.pressure.valid_mantissa_bits = args.mantissa_size
85 particle.smoothL.valid_mantissa_bits = args.mantissa_size
86 particle.u.valid_mantissa_bits = args.mantissa_size
87 particle.uDot.valid_mantissa_bits = args.mantissa_size
88 particle.f.valid_mantissa_bits = args.mantissa_size
89 particle.wcount_dh.valid_mantissa_bits = args.mantissa_size
90 particle.rho_dh.valid_mantissa_bits = args.mantissa_size
91 particle.wcount.valid_mantissa_bits = args.mantissa_size
92 particle.hDot.valid_mantissa_bits = args.mantissa_size
93 particle.balsara.valid_mantissa_bits = args.mantissa_size
94 particle.rot_v.valid_mantissa_bits = args.mantissa_size
95 particle.div_v.valid_mantissa_bits = args.mantissa_size
96 particle.v_sig_AV.valid_mantissa_bits = args.mantissa_size
97 particle.soundSpeed.valid_mantissa_bits = args.mantissa_size
98 particle.v_full.valid_mantissa_bits = args.mantissa_size
99 particle.u_full.valid_mantissa_bits = args.mantissa_size
100 ~~~~~~~~~~~~~~~~~~~~~~~~~~
101
102 This is an example. You can also add ony a few of these fields, and picking
103 a value for args.mantissa_size is obviously up to you. We refrain from
104 illustrating that you can also alter the precision of the position and
105 maximum search radius of a particle this way. This is on purpose: We have
106 made bad experience with such modifications.
107
108 ## Boundary conditions
109
110 Boundary conditions for Swift are @ref swift_boundary_conditions "discussed on a separate page".
111 If you decide that you want to plug into the drift mechanism of
112 this SPH flavour, you typically add something similar to the
113 code snippet below to your code:
114
115 ~~~~~~~~~~~~~~~~~~~~~~
116 particle.algorithm_steps_dict["Drift"].touch_vertex_last_time_kernel += " " "
117 if (::swift2::boundaryconditions::isVertexOnGlobalBoundary(marker,DomainOffset,DomainSize)) {
118 ::swift2::kernels::forAllParticles(
119 marker,
120 assignedParticles,
121 [&](const peano4::datamanagement::VertexMarker& marker, globaldata::" " " + name + " " "& assignedParticle)->void {
122 ::swift2::boundaryconditions::applyFixedBoundaryCondition(
123 assignedParticle,
124 marker,
125 DomainOffset,
126 DomainSize,
127 0.1,
128 _spacetreeId
129 );
130 }
131 );
132 }
133 " " "
134 particle.algorithm_steps_dict["Drift"].includes += " " "
135 #include "swift2/boundaryconditions/FixedBoundary.h"
136 #include "swift2/boundaryconditions/Utils.h"
137 " " "
138 particle._setup_algorithm_steps()
139 particle._setup_initialisation_steps()
140 ~~~~~~~~~~~~~~~~~~~~~~
141
142 The spaced-out syntax is just there to ensure that the Python documentation
143 syntax is not broken.
144
145
146 ## Attributes
147
148 name: String
149 To be in line with other conventions, I recommend you start with an
150 uppercase letter. This has to be a valid C++ identifier, i.e. don't
151 use any special symbols besides an underscore.
152
153 ## Vectorisation
154
155 This type can be used with aggressively optimised kernels if you follow
156 @ref page_swift_performance_optimisation Performance optimisation and use
157 the map_particle_steps_onto_separate_mesh_traversals_insert_dummy_sweeps_for_coalesced_memory_access
158 graph compiler from the Sequential suite; or literally any other
159 flavour which yields coalesced memory.
160
161 @param particle_interaction_kernel_realisation: ParticleKernelRealisation
162 Switch through various interaction variants.
163
164
165 @see peano4::datamanagement::CellMarker
166 @see peano4::datamanagement::VertexMarker
167
168 """
169
171 self,
172 name,
173 dimensions_hydro=2,
174 cfl_factor=0.1,
175 initial_time_step_size=1e-4,
176 constant_time_step_size=True,
177 swift_project_namespace="SPH",
178 particles_per_cell=0, # makes no sense, and should likely be forbidden
179 min_h=0.3,
180 max_h=0.3,
181 particle_interaction_kernel_realisation: SPHParticle.ParticleKernelRealisation = SPHParticle.ParticleKernelRealisation.USE_OUTER_GUARDS,
182 ):
183 super(SPHLeapfrogFixedSearchRadius, self).__init__(
184 name=name,
185 dimensions_hydro=dimensions_hydro,
186 cfl_factor=cfl_factor,
187 initial_time_step_size=initial_time_step_size,
188 constant_time_step_size=constant_time_step_size,
189 swift_project_namespace=swift_project_namespace,
190 particles_per_cell=particles_per_cell,
191 min_h=min_h,
192 max_h=max_h,
193 )
194
195 # Initialize default parameters.
196 # ---------------------------------------
197
198 # This is an SPH scheme.
199 # Hardcoded this for now until we can switch flavours.
201
202 # Minimal SPH parameters
203 self._alpha_av = 0.8
204 self._beta_av = 3.0
205
206 # now set these parameters as dastgen attributes.
207 # This needs to be done in a careful way so that
208 # changes to parameters outside of the initializer
209 # propagate into the generated c++ files.
211
212 # Construct particle properties
213 # --------------------------------
214
215 # Artificial viscosity scheme --------------------------------------------------
218 self.rot_v = peano4.dastgen2.Peano4DoubleArray("rot_v", "Dimensions")
219 else: # Curl is a scalar in lower dimensions
221 "rot_v"
222 ) # In 2D sims it only has 1 component (z)
226 self.data.add_attribute(self.balsara)
227 self.data.add_attribute(self.rot_v)
228 self.data.add_attribute(self.div_v)
229 self.data.add_attribute(self.v_sig_AV)
230 self.data.add_attribute(self.soundSpeed)
231
232 # Algorithm Steps ---------------------------------------------------------------
233 # actual list containing the algorithm steps
235 # actual list containing the algorithm steps for initialization
237 # the dict containing all possible algorithm steps
238 self.algorithm_steps_dict = get_algorithm_step_dict(self)
239
240 # Now store an independent instance of the algorithm steps in the object
241 # We need an independent instance for situations where we modify the
242 # algorithm steps later. For example when we are adding debugging or
243 # dependence checks for each algorithm step.
245 # same for initialization steps
247
249
250 return
251
252 def set_parameters(self):
253 """
254 This function translates "global" particle parameters which are
255 constant throughout the simulation (like CFL factor, minimal time step
256 size, viscosity parameters...) into dastgen attributes of the C++
257 particle class.
258 If you modify any of the attributes manually outside of the particle
259 initialisation, e.g. by invoking
260
261 ```
262 particle = SPHLeapfrogFixedSearchRadius(initial_time_step_size=ABC, ...)
263 particle.initial_time_step_size = XYZ
264 ```
265
266 you need to call this function manually so your changes propagate
267 into the generated C++ files.
268 """
269
270 super(SPHLeapfrogFixedSearchRadius, self).set_parameters()
271
272 const_static = dastgen2.attributes.Attribute.Qualifier.CONST_STATIC
273
274 alpha_av_attr = dastgen2.attributes.Double(
275 "alphaAV", qualifier=const_static, initval=self._alpha_av
276 )
277 self.data.add_or_replace_attribute(alpha_av_attr)
278
279 beta_av_attr = dastgen2.attributes.Double(
280 "betaAV", qualifier=const_static, initval=self._beta_av
281 )
282 self.data.add_or_replace_attribute(beta_av_attr)
283
284 return
285
287 """!
288 Set up the internal list of algorithm steps for this particle.
289 We need to maintain an individual instance of this list for cases
290 where we modify the algorithm steps later on. This can happen for
291 example when adding debugging/dependency checks later on.
292
293
294 The algorithm steps shall be a list of AlgorithmStep objects to be
295 executed in that order.
296
297 Definitions of the algorithm steps stored in the algorithm_steps_dict are
298 placed in get_algorithm_steps_dict(self). The dict is set up during
299 instantiation of this class, and reads in all available algorithm steps
300 defined in AlgorithmStepLibrary.py.
301
302 Leapfrog consists basically of four steps per particle. We first
303 determine the force. Then we update the velocity by half a timestep and move
304 the particle by a full timestep. Then the force is re-computed and the
305 second half of the velocity update is done.
306 Some variations of this KDK form re-arrange the steps executed per timestep
307 to avoid a second force loop. We can't do this here because we're not planning
308 on running with fixed time step sizes throughout the run.
309 """
310
311 steps = [
312 # self.algorithm_steps_dict["FlagBoundaryParticles"],
313 self.algorithm_steps_dict["Drift"],
314 self.algorithm_steps_dict["PredictHydro"],
315 self.algorithm_steps_dict["Density"],
316 self.algorithm_steps_dict["Force"],
317 self.algorithm_steps_dict["Kick2"],
318 self.algorithm_steps_dict["ReduceGlobalQuantities"],
319 self.algorithm_steps_dict["Kick1"],
320 ]
321
322 # make a copy of the steps, so if you modify them later (e.g. by
323 # adding dependency checks), it won't modify the steps in the dict.
324 self._algorithm_steps_algorithm_steps = [copy.deepcopy(i) for i in steps]
325 return
326
328 """!
329 Define the algorithm steps to be taken during initialization here.
330 This follows the same logic as _setup_algorithm_steps. See documentation
331 there for more info.
332
333 Make sure get_algorithm_steps_dict(self) has been called before.
334 """
335 initsteps = [
336 self.algorithm_steps_dict["FlagBoundaryParticles"],
337 self.algorithm_steps_dict["Density"],
338 self.algorithm_steps_dict["Force"],
339 # TODO MLADEN: verify this later.
340 # self.algorithm_steps_dict["Kick1"],
341 self.algorithm_steps_dict["Kick2"],
342 ]
343
344 # return deep copies of algorithms steps. Otherwise,
345 # modifications down the line (like the dependency checks)
346 # will cause trouble.
347 self._initialisation_steps_initialisation_steps = [copy.deepcopy(i) for i in initsteps]
348
349 return
350
351 @property
352 def alpha_av(self):
353 """
354
355 Viscosity parameters of the Minimal SPH model. Default SWIFT values are
356 alpha_av=0.8 and beta_av=3.
357
358 """
359 return self._alpha_av
360
361 @alpha_av.setter
362 def alpha_av(self, alpha_av):
363 if alpha_av <= 0:
364 raise ValueError("Value for the alpha_av viscosiy parameter not valid.")
365 self._alpha_av = alpha_av
366
367 @property
368 def beta_av(self):
369 return self._alpha_av
370
371 @beta_av.setter
372 def beta_av(self, beta_av):
373 if beta_av <= 0:
374 raise ValueError("Value for the beta_av viscosiy parameter not valid.")
375 self._beta_av = beta_av
376
378 """
379
380 Transform namespace into cpp format. Could be used to append namespace to
381 constants in kernels (not used currently).
382
383 """
384 namespace = "::".join(self._swift_project_namespace) + "::"
385
386 print(namespace)
387
388 @property
390 return (
391 super(SPHLeapfrogFixedSearchRadius, self).readme_descriptor
392 + """
393
394 Time integrator: SPH Leapfrog
395
396 - search radius: fixed
397 - CFL factor: """
399 + """
400 - dt_initial: """
402 + """
403
404 """
405 )
cardinality is a string (you can thus use symbols as well as long as they will be defined at compile ...
Definition Double.py:6
Specialisation of dastgen2.attributes.DoubleArray which relies on Peano's tarch.
_add_dependency_checks(self)
Add dependency (particle consistency) checks.
Definition Particle.py:482
__init__(self, name, dimensions_hydro=2, cfl_factor=0.1, initial_time_step_size=1e-4, constant_time_step_size=True, swift_project_namespace="SPH", particles_per_cell=0, min_h=0.3, max_h=0.3, SPHParticle.ParticleKernelRealisation particle_interaction_kernel_realisation=SPHParticle.ParticleKernelRealisation.USE_OUTER_GUARDS)
Initialise the particle.
_setup_initialisation_steps(self)
Define the algorithm steps to be taken during initialization here.
_setup_algorithm_steps(self)
Set up the internal list of algorithm steps for this particle.
set_parameters(self)
This function translates "global" particle parameters which are constant throughout the simulation (l...
_setup_initialisation_steps(self)
Define the algorithm steps to be taken during initialization here.
_setup_algorithm_steps(self)
Set up the internal list of algorithm steps for this particle.
set_parameters(self)
This function translates "global" particle parameters which are constant throughout the simulation (l...