Peano 4
Loading...
Searching...
No Matches
acoustic-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 = [10.0, 10.0, 10.0]
42offset = [0.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", "acousticwave"],
72 project_name="AcousticWaveTutorial",
73 directory=".",
74 executable="acoustic-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="AcousticWaveSolver",
85 patch_size=patch_size,
86 unknowns=dimensions + 1, # [p, u, v, (w)]
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 + 1, auxiliary_variables=0, dimensions=dimensions
99)
100p = my_pde.name_Q_entry(0, "p")
101v = my_pde.name_Q_entries(1, dimensions, "v")
102
103rho = sympy.symbols("rho")
104c = sympy.symbols("c")
105K0 = rho * c * c
106
107#
108# Define the equation system
109#
110# Flux [unknowns, dimensions]
111if dimensions == 2:
112 my_pde.F[0, 0] = K0 * v[0]
113 my_pde.F[1, 0] = 1.0 / rho * p
114 my_pde.F[2, 0] = 0
115
116 my_pde.F[0, 1] = K0 * v[1]
117 my_pde.F[1, 1] = 0
118 my_pde.F[2, 1] = 1.0 / rho * p
119else:
120 my_pde.F[0, 0] = K0 * v[0]
121 my_pde.F[1, 0] = 1.0 / rho * p
122 my_pde.F[2, 0] = 0
123 my_pde.F[3, 0] = 0
124
125 my_pde.F[0, 1] = K0 * v[1]
126 my_pde.F[1, 1] = 0
127 my_pde.F[2, 1] = 1.0 / rho * p
128 my_pde.F[3, 1] = 0
129
130 my_pde.F[0, 2] = K0 * v[2]
131 my_pde.F[1, 2] = 0
132 my_pde.F[2, 2] = 0
133 my_pde.F[3, 2] = 1.0 / rho * p
134
135# Eigenvalues [unknowns, dimensions]
136if dimensions == 2:
137 my_pde.eigenvalues[0] = -c
138 my_pde.eigenvalues[1] = 0
139 my_pde.eigenvalues[2] = c
140else:
141 my_pde.eigenvalues[0] = -c
142 my_pde.eigenvalues[1] = 0
143 my_pde.eigenvalues[2] = 0
144 my_pde.eigenvalues[3] = c
145
146#
147# Since 'my_pde' only holds the PDE for the advection equation without initial- or boundary conditions,
148# we still need to properly define initial- and boundary conditions.
149# This gives us then a complete description of a 'scenario'.
150#
151#
152# Initial conditions
153#
154my_pde.initial_values[0] = 0
155my_pde.initial_values[1] = 0
156my_pde.initial_values[2] = 0
157if dimensions == 3:
158 my_pde.initial_values[3] = 0
159
160#
161# Boundary conditions
162#
163my_pde.boundary_values[0] = 0
164my_pde.boundary_values[1] = 0
165my_pde.boundary_values[2] = 0
166if dimensions == 3:
167 my_pde.boundary_values[3] = 0
168
169#
170# Source term
171#
172sigma = sympy.symbols("sigma")
173t0 = sympy.symbols("t0")
174M0 = sympy.symbols("M0")
175t = sympy.symbols("t")
176
177max_h = sympy.symbols("MaxAdmissibleVolumeH")
178force = M0 * sympy.exp(-((t - t0) * (t - t0)) / (2.0 * sigma * sigma))
179
180if dimensions == 2:
181 point_source = sympy.sqrt((5 - my_pde.x[0]) ** 2 + (5 - my_pde.x[1]) ** 2)
182 my_pde.sources[0] = sympy.Piecewise((force, point_source <= max_h), (0, True)) # p
183else:
184 point_source = sympy.sqrt(
185 (5 - my_pde.x[0]) ** 2 + (5 - my_pde.x[1]) ** 2 + (5 - my_pde.x[2]) ** 2
186 )
187 my_pde.sources[0] = sympy.Piecewise((force, point_source <= max_h), (0, True)) # p
188
189my_pde.sources[1] = 0
190my_pde.sources[2] = 0
191if dimensions == 3:
192 my_pde.sources[3] = 0
193
194my_pde.substitute_expression(rho, 2.7)
195my_pde.substitute_expression(c, 6.0)
196my_pde.substitute_expression(sigma, 0.1149)
197my_pde.substitute_expression(t0, 0.7)
198my_pde.substitute_expression(M0, 1000.0)
199
200#
201# Specify which implementation our solvers uses.
202# Here we want to set the implementation we get from our symbolically defined PDE,
203# i.e., we get the C++ implementation which is generated by ExaHyPE's 'symhype' package.
204#
205my_solver.set_implementation(
206 initial_conditions=my_pde.implementation_of_initial_conditions(),
207 boundary_conditions=my_pde.implementation_of_boundary_conditions(),
208 flux=my_pde.implementation_of_flux(),
209 eigenvalues=my_pde.implementation_of_max_eigenvalue(),
210 source_term=my_pde.implementation_of_sources(),
211)
212
213#
214# To see which variables (unknowns + auxiliary variables) we can visualise,
215# let's add a plot description for the variables to our solver.
216#
217my_solver.plot_description = my_pde.unknown_identifier_for_plotter()
218
219#
220# Add the solver to our project
221#
222my_project.add_solver(my_solver)
223
224#
225# Configure some global parameters
226#
227my_project.set_global_simulation_parameters(
228 dimensions=dimensions,
229 size=size[0:dimensions],
230 offset=offset[0:dimensions],
231 min_end_time=end_time,
232 max_end_time=end_time,
233 first_plot_time_stamp=0.0,
234 time_in_between_plots=time_in_between_two_snapshots,
235 periodic_BC=[False, False, False],
236)
237
238#
239# This defines where the output files should go.
240# If you omit this, output files are automatically put into the application's folder.
241#
242my_project.set_output_path("solution")
243if not os.path.exists("solution"):
244 os.makedirs("solution")
245
246if visualise:
247 output_patch_file = "solution/solution-AcousticWaveSolver.peano-patch-file"
248 if os.path.isfile(output_patch_file):
249 subprocess.call(
250 [
251 "python3",
252 "../../../python/peano4/visualisation/render.py",
253 str(output_patch_file),
254 ]
255 )
256 sys.exit(0)
257
258#
259# If you only target a sequential execution, you can omit the load balancing.
260# However, for a parallel build and execution you need to set load balancing.
261# This requires the 'LoadBalancing' toolbox to be built.
262# The type of load balancer can greatly impact the speedup and overall performance.
263# For an overview of available load balancer refer to the documentation.
264#
265my_project.set_load_balancing(
266 "toolbox::loadbalancing::strategies::SpreadOutHierarchically",
267 "new ::exahype2::LoadBalancingConfiguration()",
268)
269
270#
271# We need to set the location of our core libraries ('Peano4').
272# This helps us to resolve any dependencies.
273# Additionally, we specify the build mode which you can also change to a different mode.
274#
275my_project.set_Peano4_installation("../../../", mode=compile_mode)
276
277#
278# We generate and grab the underlying core project of 'Peano4'.
279# This gives us access to some functions we want to use to finalise and build this project.
280#
281my_project = my_project.generate_Peano4_project(verbose=True)
282
283#
284# Optional: Set an floating-point exception handler to the Peano 4 project.
285#
286my_project.set_fenv_handler("FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW")
287
288#
289# Finally, we want to build our project.
290# First, all of the necessary glue code is generated in the application folder,
291# then 'make' is invoked automatically which compiles the generated code and links against our core libraries
292# and toolboxes which have been built before.
293# You can also always invoke 'make' yourself to compile, or cleanup with 'make clean'.
294#
295my_project.build(make=use_make, make_clean_first=True, throw_away_data_after_build=True)
296
297print("\nPDE is: ")
298print(my_pde.__str__())
299print("Executable is acoustic-wave-" + str(dimensions) + "d")
300print("Clean object files via 'make clean'")
301print("(Re)compile the generated code via 'make -j'")
302print("Remove all generated code via 'make distclean'")
303print("Regenerate all code by running 'python3 acoustic-wave.py' again")
ExaHyPE 2 project.
Definition Project.py:18
Helper class to model a hyperbolic PDE in first-order conservative formulation.