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