Peano 4
Loading...
Searching...
No Matches
planet-orbit.py
Go to the documentation of this file.
1import os
2import argparse
3import numpy as np
4
5import peano4
6import swift2
7import dastgen2
8
9import shutil
10
11
14
15parser = argparse.ArgumentParser(description="SPH benchmarking script")
16parser.add_argument(
17 "-j", "--parallel-builds", dest="j", type=int, default=-1, help="Parallel builds"
18)
19parser.add_argument(
20 "-v",
21 "--verbose",
22 dest="verbose",
23 action="store_true",
24 default=False,
25 help="Verbose",
26)
27parser.add_argument(
28 "-c",
29 "--clean",
30 dest="clean_temporary_directories",
31 action="store_true",
32 help="Remove temporary directories before script starts",
33)
34parser.add_argument(
35 "-d",
36 "--dimensions",
37 dest="dimensions",
38 type=int,
39 default=2,
40 help="Spatial dimensions (default 2)",
41)
42parser.add_argument(
43 "-et",
44 "--end-time",
45 dest="end_time",
46 type=float,
47 default=1.0,
48 help="End of simulation",
49)
50parser.add_argument(
51 "-dt",
52 "--timestep-size",
53 dest="timestep_size",
54 type=float,
55 default=1e-3,
56 help="Global timestep size (initial value)",
57)
58parser.add_argument(
59 "-plot",
60 "--plot-delta",
61 dest="plot_delta",
62 type=float,
63 default=0.01,
64 help="Plot delta (0 switches off)",
65)
66parser.add_argument(
67 "-cs",
68 "--cell-size",
69 dest="cell_size",
70 type=float,
71 default=0.3,
72 help="Cell size (default 0.3)",
73)
74parser.add_argument(
75 "-pd", "--peano-dir", dest="peanodir", default="../../..", help="Peano4 directory"
76)
77parser.add_argument(
78 "-asserts",
79 "--asserts",
80 dest="asserts",
81 action="store_true",
82 default=False,
83 help="Switch on assertions",
84)
85parser.add_argument(
86 "-integrator",
87 "--time-integrator",
88 dest="time_integrator",
89 choices=["leapfrog", "euler", "euler-dynamic", "Euler", "Euler-dynamic"],
90 required=True,
91 help="Type of time integrator to use (default: leapfrog)",
92)
93args = parser.parse_args()
94
95
96
99
100dimensions = args.dimensions
101
102
103#! [tutorials swift2 planet-orbit Create project]
105 namespace=["tests", "swift2", "planetorbit"],
106 project_name="planetorbit",
107 executable="orbit-" + str(args.cell_size),
108)
109#! [tutorials swift2 planet-orbit Create project]
110
111
112"""
113Declare the particle type:
114We need to specify the time integration scheme that will be used to evolve
115this particle. Current choices include:
116 - Leapfrog (default)
117 - Explicit Euler with fixed search radius
118 - Explicit Euler with variable search radius
119
120"""
121
122if args.time_integrator == "leapfrog":
123 #! [tutorials swift2 planet-orbit Instantiate particle with leapfrog ODE integrator]
125 name="Planet",
126 cfl_factor=0.2,
127 initial_time_step_size=args.timestep_size,
128 particles_per_cell=0,
129 min_h=args.cell_size,
130 max_h=0.3,
131 )
132 particle_set = project.add_particle_species(particle)
133 #! [tutorials swift2 planet-orbit Instantiate particle with leapfrog ODE integrator]
134elif args.time_integrator == "euler" or "Euler":
136 name="Planet",
137 cfl_factor=0.2,
138 initial_time_step_size=args.timestep_size,
139 particles_per_cell=0,
140 min_h=args.cell_size,
141 max_h=0.3,
142 )
143 particle_set = project.add_particle_species(particle)
144elif args.time_integrator == "euler-dynamic" or "Euler-dynamic":
146 name="Planet",
147 cfl_factor=0.2,
148 initial_time_step_size=args.timestep_size,
149 particles_per_cell=0,
150 min_h=args.cell_size,
151 max_h=0.3,
152 )
153 particle_set = project.add_particle_species(particle)
154else:
155 print("Time integrator not supported. Please check the list of options.")
156
157
158#! [tutorials swift2 planet-orbit Declare additional attributes]
159energy_kin_attr = dastgen2.attributes.Double("energyKin")
160energy_pot_attr = dastgen2.attributes.Double("energyPot")
161energy_tot_attr = dastgen2.attributes.Double("energyTot")
162particle.data.add_attribute(energy_kin_attr)
163particle.data.add_attribute(energy_pot_attr)
164particle.data.add_attribute(energy_tot_attr)
165#! [tutorials swift2 planet-orbit Declare additional attributes]
166
167
168#! [tutorials swift2 planet-orbit Configure project]
169if args.asserts:
170 build_mode = peano4.output.CompileMode.Asserts
171else:
172 build_mode = peano4.output.CompileMode.Release
173
174offset = [0, 0, 0]
175domain_size = [1, 1, 1]
176periodic_boundary_conditions = [False, False, False] # Periodic BC
177
178project.set_global_simulation_parameters(
179 dimensions=dimensions,
180 offset=offset,
181 domain_size=domain_size,
182 min_end_time=args.end_time,
183 max_end_time=0.0,
184 first_plot_time_stamp=0.0,
185 time_in_between_plots=args.plot_delta,
186 periodic_BC=periodic_boundary_conditions,
187 plotter_precision=8, # plotter precision
188)
189
190project.set_load_balancing(
191 "toolbox::loadbalancing::strategies::SpreadOutHierarchically",
192 "new toolbox::loadbalancing::DefaultConfiguration()",
193)
194
195project.set_Peano4_installation(args.peanodir, build_mode)
196#! [tutorials swift2 planet-orbit Configure project]
197
198
199
200#! [tutorials swift2 planet-orbit Define constants]
201# Orbital parameters
202N_orbits = 8
203# Number of orbits
204G_NEWTON = 6.67384e-11
205# Newton's constant [m^3/kg/s^2]
206M_SUN = 1.9885e30
207# Sun mass [kg]
208M_EARTH = 5.97219e24
209# Earth mass [kg]
210R_MAX = 152097701000.0
211# [m]
212R_MIN = 147098074000.0
213# [m]
214V_MAX = 30287.0
215# [m/s]
216
217# Initial positions in internal units
218LENGTH_UNIT = 4.0 * R_MAX
219SUN_X_COORD = str(0.5)
220SUN_Y_COORD = str(0.5)
221EARTH_X_INI = str(float(SUN_X_COORD) + R_MAX / LENGTH_UNIT)
222EARTH_Y_INI = SUN_Y_COORD
223
224# Derived quantities
225e = (R_MAX - R_MIN) / (R_MAX + R_MIN)
226# Eccentricity [dimensionless]
227b = np.sqrt(R_MAX * R_MIN)
228# Semi-minor axis [m]
229p = b * np.sqrt(1.0 - e * e)
230# Semi-lactus rectum [m]
231a = p / (1.0 - e * e)
232# Semi-major axis [m]
233T = np.sqrt(4.0 * np.pi * np.pi * a * a * a / (G_NEWTON * (M_SUN + M_EARTH)))
234# Period [s] */
235
236# Time-step size
237dt = 1e-3 * T
238# [s]
239dt /= T
240N = N_orbits / dt
241GLOBAL_TIME_STEP_SIZE = str(dt)
242
243# Initial velocity in internal units
244EARTH_V_INI_X = str(0.0)
245EARTH_V_INI_Y = str(V_MAX / LENGTH_UNIT * T)
246#! [tutorials swift2 planet-orbit Define constants]
247
248
249
250
251#! [tutorials swift2 planet-orbit Particle force]
252
289particle_force = """
290::swift2::kernels::forAllLocalParticles(
291 marker,
292 localParticles,
293 [=](
294 const peano4::datamanagement::CellMarker& marker,
295 globaldata::Planet& localParticle
296 ) -> void {
297 if ( ::swift2::kernels::localParticleCanBeUpdatedInCellKernel(
298 marker,
299 localParticle
300 )) {
301 #if Dimensions==2
302 tarch::la::Vector<Dimensions,double> pos_sun = {SUN_X_COORD,SUN_Y_COORD};
303 #else
304 tarch::la::Vector<Dimensions,double> pos_sun = {SUN_X_COORD,SUN_Y_COORD,0.0};
305 #endif
306
307 /* G_NEWTON in internal units */
308 double const L3 = LENGTH_UNIT * LENGTH_UNIT * LENGTH_UNIT;
309 double const T2 = PERIOD * PERIOD;
310 double const GN = G_NEWTON * M_SUN * T2 / L3;
311 //double GN = 1.0;
312
313 tarch::la::Vector<Dimensions,double> earthSunDistance = localParticle.getX() - pos_sun;
314 const double r = tarch::la::norm2(earthSunDistance);
315 if (tarch::la::greater(r,0.0)) {
316 const double r_inv = 1.0 / r;
317 localParticle.setA( - GN * earthSunDistance * (r_inv * r_inv * r_inv) );
318
319 /* Update Energy */
320 const double vel = tarch::la::norm2( localParticle.getV() );
321 const double E_kin = 0.5 * vel * vel;
322 const double E_pot = - GN * r_inv;
323
324 /* Update total energy (overall minus sign) */
325 localParticle.setEnergyKin( E_kin );
326 localParticle.setEnergyPot(-E_pot );
327 localParticle.setEnergyTot( - ( E_kin + E_pot ) );
328 } // end of radius check
329 } // end of "may I update" check
330 } // end of functor
331);
332"""
333
334particle.cell_kernel = particle_force
335#! [tutorials swift2 planet-orbit Particle force]
336
337
338#! [tutorials swift2 planet-orbit Initial conditions]
339
349search_radius = args.cell_size / 100.0
350initialisation_snippet = (
351 """
352 /*
353 * In this problem we will initialize just a single particle (Earth)
354 * This will orbit around a point where the hypothetical Sun is (no particle needed here)
355 */
356
357 /* Internal search radius used by Peano.
358 * Should always be larger than actual SPH interaction radius (~2h).
359 * But no really used here. */
360 particle->setSearchRadius("""
361 + str(search_radius)
362 + """);
363 particle->setV( {0, 0, 0} );
364 particle->setA( {0, 0, 0} );
365
366 // Orbiting particle: We can also have a free moving sun, but by default, this
367 // one rests and the earth is actually orbiting.
368 if( particle->getX(0) == EARTH_X_INI ){
369 particle->setV(0, EARTH_V_INI_X) ;
370 particle->setV(1, EARTH_V_INI_Y) ;
371 }
372
373 /*
374 * Set initial energy
375 */
376
377 /* G_NEWTON in internal units */
378 double const L3 = LENGTH_UNIT * LENGTH_UNIT * LENGTH_UNIT;
379 double const T2 = PERIOD * PERIOD;
380 double const GN = G_NEWTON * M_SUN * T2 / L3;
381
382 /* Sun-Earth initial distance */
383 const double r_ini = particle->getX()(0) - SUN_X_COORD;
384 const double r_inv = r_ini ? 1.0 / r_ini : 0.0;
385 const double vel = tarch::la::norm2( particle->getV() );
386
387 const double E_kin = 0.5 * vel * vel;
388 const double E_pot = - GN * r_inv;
389
390 particle->setEnergyKin( E_kin );
391 particle->setEnergyPot(-E_pot );
392
393 /* Overall minus sign for plotting */
394 particle->setEnergyTot( -( E_kin + E_pot ) );
395
396 logInfo( "Init()", "=====================================" );
397 logInfo( "Init()", "Initial energy:");
398 logInfo( "Init()", "Kinetic energy :" << E_kin);
399 logInfo( "Init()", "Potential energy:" << E_pot);
400 logInfo( "Init()", "Total energy :" << E_kin + E_pot);
401 logInfo( "Init()", "Search Radius :" << particle->getSearchRadius() );
402 logInfo( "Init()", "Initial velocity:" << particle->getV() );
403 logInfo( "Init()", "=====================================" );
404
405
406"""
407)
408
409#
410# Insert Earth into system
411#
412project.algorithm_step_initial_conditions.add_action_set(
414 particle_set=particle_set,
415 initialisation_call=initialisation_snippet,
416 coordinates=[(EARTH_X_INI, EARTH_Y_INI, 0)],
417 )
418)
419#! [tutorials swift2 planet-orbit Initial conditions]
420
421
422#! [tutorials swift2 planet-orbit Plot]
423
447 "particles",
448 particle_set,
449 "globaldata::{}::getSpecies().getMinTimeStamp()".format(particle.name),
450)
451particle_plotter.add_attribute_to_plot(particle.velocity, dimensions)
452particle_plotter.add_attribute_to_plot(energy_kin_attr, 1)
453particle_plotter.add_attribute_to_plot(energy_pot_attr, 1)
454particle_plotter.add_attribute_to_plot(energy_tot_attr, 1)
455project.algorithm_step_plot.add_action_set(particle_plotter)
456#! [tutorials swift2 planet-orbit Plot]
457
458
459
460
461#! [tutorials swift2 planet-orbit Validation]
462particle.add_to_reduction( f"""
463nonCriticalAssertion( vertexdata::{particle.name}Set::getNumberOfRemainingLocalParticles() == 1 );
464""")
465#! [tutorials swift2 planet-orbit Validation]
466
467
468
469
479peano4_project = project.generate_Peano4_project(verbose=args.verbose)
480#! [tutorials swift2 planet-orbit Create Peano project]
481
482#! [tutorials swift2 planet-orbit Export constants]
483
484# Export symbols and values for compilation (e.g. for SPH kernel)
485# Dummy values for non SPH
486peano4_project.output.makefile.add_CXX_flag("-D" + "HYDRO_DIMENSION=1")
487peano4_project.output.makefile.add_CXX_flag("-D" + "HYDRO_GAMMA_5_3")
488peano4_project.output.makefile.add_CXX_flag("-D" + "QUARTIC_SPLINE_KERNEL")
489
490peano4_project.constants.export_const_with_type(
491 "GLOBAL_TIME_STEP_SIZE", GLOBAL_TIME_STEP_SIZE, "double"
492)
493peano4_project.constants.export_const_with_type("G_NEWTON", str(G_NEWTON), "double")
494peano4_project.constants.export_const_with_type(
495 "LENGTH_UNIT", str(LENGTH_UNIT), "double"
496)
497peano4_project.constants.export_const_with_type("M_SUN", str(M_SUN), "double")
498peano4_project.constants.export_const_with_type("SUN_X_COORD", SUN_X_COORD, "double")
499peano4_project.constants.export_const_with_type("SUN_Y_COORD", SUN_Y_COORD, "double")
500peano4_project.constants.export_const_with_type("EARTH_X_INI", EARTH_X_INI, "double")
501peano4_project.constants.export_const_with_type("EARTH_Y_INI", EARTH_Y_INI, "double")
502peano4_project.constants.export_const_with_type(
503 "EARTH_V_INI_X", EARTH_V_INI_X, "double"
504)
505peano4_project.constants.export_const_with_type(
506 "EARTH_V_INI_Y", EARTH_V_INI_Y, "double"
507)
508peano4_project.constants.export_const_with_type("ECCENTRICITY", str(e), "double")
509peano4_project.constants.export_const_with_type("PERIOD", str(T), "double")
510peano4_project.constants.export_const_with_type("R_MIN", str(R_MIN), "double")
511peano4_project.constants.export_const_with_type("R_MAX", str(R_MAX), "double")
512#! [tutorials swift2 planet-orbit Export constants]
513
514# If we want to start a project from scratch, we may want to remove some of the generated code
515if args.clean_temporary_directories:
516 dirs = ["vertexdata/", "repositories/", "globaldata/", "observers/"]
517 for dir in dirs:
518 shutil.rmtree(dir)
519
520peano4_project.generate(
521 overwrite=peano4.output.Overwrite.Always, throw_away_data_after_generation=True
522)
523peano4_project.build(make_clean_first=True, number_of_parallel_builds=args.j)
524
525# remove output files from previous runs
526output_files = [
527 f
528 for f in os.listdir(".")
529 if f.endswith(".peano-patch-file") or f.endswith(".vtu") or f.endswith(".pvd")
530]
531for f in output_files:
532 os.remove(f)
cardinality is a string (you can thus use symbols as well as long as they will be defined at compile ...
Definition Double.py:6
Swift2 project.
Definition Project.py:23
Extension of the explicit Euler to support dynamic search radius adoption.
Simple particle with fixed search radius subject to explicit Euler.