Peano 4
Loading...
Searching...
No Matches
elastic-wave.py
Go to the documentation of this file.
1# This file is part of the ExaHyPE2 project. For conditions of distribution and
2# use, please see the copyright notice at www.peano-framework.org
3import os
4import sys
5import sympy
6import subprocess
7
8#
9# We import the required modules for this project.
10# We always need the 'peano4' module as this is our core project.
11# Since we are creating an ExaHyPE 2 application, we additionally
12# need to import the 'exahype2' module.
13#
14import peano4
15import exahype2
16
17#
18# We specify the space dimensions here.
19# We support either 2- or 3-dimensional problems.
20#
21dimensions = 2
22
23#
24# The number of finite volumes per axis in one patch.
25#
26patch_size = 16
27
28#
29# The size of a finite volume/cell per axis.
30#
31volume_size = 0.1
32
33#
34# The simulation end time.
35#
36end_time = 1
37
38#
39# Choose domain size and offset.
40#
41size = [20.0, 20.0, 20.0]
42offset = [-5.0, 0.0, 0.0]
43
44#
45# If set to False, the project will not compile the project after generation (i.e., you need to invoke the Makefile yourself).
46#
47use_make = True
48
49#
50# Set to true, to visualise/convert the output data using a postprocessor.
51#
52visualise = False
53
54#
55# Choose how often a snapshot is written.
56#
57time_in_between_two_snapshots = end_time / 100
58# time_in_between_two_snapshots = 0 # Set to 0 to disable I/O
59
60#
61# Switch between 'Release', 'Debug', 'Asserts', 'Trace', 'Stats'.
62#
63compile_mode = peano4.output.CompileMode.Release
64
65#
66# We first create a new ExaHyPE 2 project.
67# For this, we specify the (nested) namespaces, the name of our main file and our executable name.
68# We amend the executable name by the space dimensions.
69#
70my_project = exahype2.Project(
71 namespace=["tutorials", "exahype2", "elasticwave"],
72 project_name="ElasticWaveTutorial",
73 directory=".",
74 executable="elastic-wave-" + str(dimensions) + "d",
75)
76
77#
78# Add the Finite Volumes solver using named arguments.
79# This is the way you can add further PDE terms.
80# This requires the 'BlockStructured' toolbox and 'ExaHyPE' to be built.
81# Rusanov is the type of flux that is used to solve the Riemann problem at boundaries between cells.
82#
84 name="ElasticWaveSolver",
85 patch_size=patch_size,
86 unknowns=dimensions + 3, # [u, v, sxx, syy, sxy]
87 auxiliary_variables=0,
88 min_volume_h=volume_size,
89 max_volume_h=volume_size,
90 time_step_relaxation=0.5,
91)
92
93#
94# We want to define our PDE symbolically.
95# For this we use the 'symhype' package (not to be confused with 'sympy') from 'ExaHyPE'.
96#
98 unknowns=dimensions + 3, auxiliary_variables=0, dimensions=dimensions
99)
100v = my_pde.name_Q_entries(0, dimensions, "v")
101sigma = my_pde.name_Q_entries(dimensions, 3, "sigma")
102
103rho = sympy.symbols("rho")
104cp = sympy.symbols("cp")
105cs = sympy.symbols("cs")
106
107# Lamé parameters
108mu = rho * cs * cs
109lamb = rho * cp * cp - 2.0 * mu
110
111#
112# Define the equation system
113#
114# Flux [unknowns, dimensions]
115my_pde.F[0, 0] = -1.0 / rho * sigma[0]
116my_pde.F[1, 0] = -1.0 / rho * sigma[2]
117my_pde.F[2, 0] = -(2.0 * mu + lamb) * v[0]
118my_pde.F[3, 0] = -lamb * v[0]
119my_pde.F[4, 0] = -mu * v[1]
120
121my_pde.F[0, 1] = -1.0 / rho * sigma[2]
122my_pde.F[1, 1] = -1.0 / rho * sigma[1]
123my_pde.F[2, 1] = -lamb * v[1]
124my_pde.F[3, 1] = -(2.0 * mu + lamb) * v[1]
125my_pde.F[4, 1] = -mu * v[0]
126
127# Eigenvalues [unknowns, dimensions]
128my_pde.eigenvalues[0] = cp
129my_pde.eigenvalues[1] = cs
130my_pde.eigenvalues[2] = 0
131my_pde.eigenvalues[3] = -cs
132my_pde.eigenvalues[4] = -cp
133
134#
135# Since 'my_pde' only holds the PDE for the advection equation without initial- or boundary conditions,
136# we still need to properly define initial- and boundary conditions.
137# This gives us then a complete description of a 'scenario'.
138#
139#
140# Initial conditions
141#
142my_pde.initial_values[0] = 0
143my_pde.initial_values[1] = 0
144my_pde.initial_values[2] = 0
145my_pde.initial_values[3] = 0
146my_pde.initial_values[4] = 0
147
148#
149# Boundary conditions
150#
151my_pde.boundary_values[0] = 0
152my_pde.boundary_values[1] = 0
153my_pde.boundary_values[2] = 0
154my_pde.boundary_values[3] = 0
155my_pde.boundary_values[4] = 0
156
157#
158# Source term
159#
160t0 = sympy.symbols("t0")
161M0 = sympy.symbols("M0")
162t = sympy.symbols("t")
163force = M0 * t / (t0 * t0) * (sympy.exp(-t / t0))
164
165max_h = sympy.symbols("MaxAdmissibleVolumeH")
166point_source = sympy.sqrt((10 - my_pde.x[0]) ** 2 + (10 - my_pde.x[1]) ** 2)
167
168my_pde.sources[0] = 0
169my_pde.sources[1] = 0
170my_pde.sources[2] = 0
171my_pde.sources[3] = 0
172my_pde.sources[4] = sympy.Piecewise((force, point_source <= max_h), (0, True))
173
174my_pde.substitute_expression(rho, 2.7)
175my_pde.substitute_expression(cp, 6.0)
176my_pde.substitute_expression(cs, 3.464)
177my_pde.substitute_expression(t0, 0.1)
178my_pde.substitute_expression(M0, 1000.0)
179
180#
181# Specify which implementation our solvers uses.
182# Here we want to set the implementation we get from our symbolically defined PDE,
183# i.e., we get the C++ implementation which is generated by ExaHyPE's 'symhype' package.
184#
185my_solver.set_implementation(
186 initial_conditions=my_pde.implementation_of_initial_conditions(),
187 boundary_conditions=my_pde.implementation_of_boundary_conditions(),
188 flux=my_pde.implementation_of_flux(),
189 eigenvalues=my_pde.implementation_of_max_eigenvalue(),
190 source_term=my_pde.implementation_of_sources(),
191)
192
193#
194# To see which variables (unknowns + auxiliary variables) we can visualise,
195# let's add a plot description for the variables to our solver.
196#
197my_solver.plot_description = my_pde.unknown_identifier_for_plotter()
198
199#
200# Add the solver to our project
201#
202my_project.add_solver(my_solver)
203
204#
205# Configure some global parameters
206#
207my_project.set_global_simulation_parameters(
208 dimensions=dimensions,
209 size=size[0:dimensions],
210 offset=offset[0:dimensions],
211 min_end_time=end_time,
212 max_end_time=end_time,
213 first_plot_time_stamp=0.0,
214 time_in_between_plots=time_in_between_two_snapshots,
215 periodic_BC=[False, False, False],
216)
217
218#
219# This defines where the output files should go.
220# If you omit this, output files are automatically put into the application's folder.
221#
222my_project.set_output_path("solution")
223if not os.path.exists("solution"):
224 os.makedirs("solution")
225
226if visualise:
227 output_patch_file = "solution/solution-ElasticWaveSolver.peano-patch-file"
228 if os.path.isfile(output_patch_file):
229 subprocess.call(
230 [
231 "python3",
232 "../../../python/peano4/visualisation/render.py",
233 str(output_patch_file),
234 ]
235 )
236 sys.exit(0)
237
238#
239# If you only target a sequential execution, you can omit the load balancing.
240# However, for a parallel build and execution you need to set load balancing.
241# This requires the 'LoadBalancing' toolbox to be built.
242# The type of load balancer can greatly impact the speedup and overall performance.
243# For an overview of available load balancer refer to the documentation.
244#
245my_project.set_load_balancing(
246 "toolbox::loadbalancing::strategies::SpreadOutHierarchically",
247 "new ::exahype2::LoadBalancingConfiguration()",
248)
249
250#
251# We need to set the location of our core libraries ('Peano4').
252# This helps us to resolve any dependencies.
253# Additionally, we specify the build mode which you can also change to a different mode.
254#
255my_project.set_Peano4_installation("../../../", mode=compile_mode)
256
257#
258# We generate and grab the underlying core project of 'Peano4'.
259# This gives us access to some functions we want to use to finalise and build this project.
260#
261my_project = my_project.generate_Peano4_project(verbose=True)
262
263#
264# Optional: Set an floating-point exception handler to the Peano 4 project.
265#
266my_project.set_fenv_handler("FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW")
267
268#
269# Finally, we want to build our project.
270# First, all of the necessary glue code is generated in the application folder,
271# then 'make' is invoked automatically which compiles the generated code and links against our core libraries
272# and toolboxes which have been built before.
273# You can also always invoke 'make' yourself to compile, or cleanup with 'make clean'.
274#
275my_project.build(make=use_make, make_clean_first=True, throw_away_data_after_build=True)
276
277print("\nPDE is: ")
278print(my_pde.__str__())
279print("Executable is elastic-wave-" + str(dimensions) + "d")
280print("Clean object files via 'make clean'")
281print("(Re)compile the generated code via 'make -j'")
282print("Remove all generated code via 'make distclean'")
283print("Regenerate all code by running 'python3 elastic-wave.py' again")
ExaHyPE 2 project.
Definition Project.py:18
Helper class to model a hyperbolic PDE in first-order conservative formulation.