Peano 4
Loading...
Searching...
No Matches
FirstOrderConservativePDEFormulation.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 sympy
4import numbers
5
6from .PDE import PDE
7
8
10 """
11 Helper class to model a hyperbolic PDE in first-order conservative
12 formulation.
13
14 To model your PDE, you typically run through a couple of single steps.
15 First, include this package plus sympy. Next, create a new instance
16 of this class to which you pass the number of equations in your PDE
17 system, i.e., the number of unknowns you want to develop, plus the
18 dimension.
19
20 As a consequence, the solver holds an equation pde which describes the
21 PDE's right-hand side in the representation
22
23 \partial Q(t) + div_x F(Q)
24
25 Most people prefer not to work with a vector Q but with some symbolic
26 names. Feel free to derive them from Q via constructs similar to
27
28 rho = euler.Q[0]
29 j = sympy.Array(euler.Q[1:4])
30 E = euler.Q[4]
31 """
32
33 def __init__(self, unknowns, auxiliary_variables, dimensions):
34 PDE.__init__(self, unknowns, auxiliary_variables, dimensions)
35
36 self.F = sympy.symarray("F", (unknowns, dimensions))
37 self.ncp = sympy.symarray("ncp", (unknowns, dimensions))
38 self.eigenvalues = sympy.symarray("alpha", (unknowns, dimensions))
39 self.sources = sympy.symarray("S", (unknowns))
40
41 def __str__(self):
42 result = ""
43 for i in range(0, self.unknownsunknowns):
44 result += "d_t(" + str(self.Q[i]) + "(x)) + div("
45 result += str(self.F[i])
46 if isinstance(self.ncp[i], numbers.Number):
47 result += " + "
48 result += str(self.ncp[i])
49 result += ") = "
50 if isinstance(self.sources[i], numbers.Number) or isinstance(self.sources[i], sympy.Expr):
51 result += str(self.sources[i]) + "\n"
52 else:
53 result += "0\n"
54 result += "\n"
55 for i in range(0, self.unknownsunknowns):
56 result += "\lambda_{max," + str(i) + "}(x) = "
57 result += str(self.eigenvalues[i])
58 result += "\n"
59 return result
60
61 def LaTeX(self):
62 result = "\\begin{eqnarray*}"
63 for i in range(0, self.unknownsunknowns):
64 result += "\partial _t " + str(self.Q[i]) + "(x) + div("
65 for j in range(0, len(self.F[i])):
66 if j != 0:
67 result += ", "
68 result += str(self.F[i][j])
69 if isinstance(self.ncp[i], numbers.Number):
70 result += " + "
71 result += str(self.ncp[i][j])
72 result += ") & = & "
73 if isinstance(self.sources[i], numbers.Number) or isinstance(self.sources[i], sympy.Expr):
74 result += str(self.sources[i]) + " \\\\"
75 else:
76 result += "0 \\\\"
77 for i in range(0, self.unknownsunknowns):
78 result += "\\lambda _{max," + str(i) + "}(x) & = & ["
79 for j in range(0, len(self.eigenvalues[i])):
80 if j != 0:
81 result += ","
82 result += str(self.eigenvalues[i][j])
83 result += "] \\\\"
84 result += "\\end{eqnarray*}"
85 result = result.replace("**", "^")
86 result = result.replace("sqrt", "\\sqrt")
87 return result
88
89 def substitute_expression(self, expression, new_expression):
90 """
91 Usually used to set a symbolic variable (expression) to a value
92 (new_expression). We run over all internal expressions and set
93 them accordingly.
94 """
95 for i in range(0, self.unknownsunknowns):
96 if not isinstance(self.sources[i], numbers.Number):
97 self.sources[i] = self.sources[i].subs(expression, new_expression)
98 for d in range(0, self.dimensionsdimensions):
99 if not isinstance(self.F[i, d], numbers.Number):
100 self.F[i, d] = self.F[i, d].subs(expression, new_expression)
101 if not isinstance(self.ncp[i, d], numbers.Number):
102 self.ncp[i, d] = self.ncp[i, d].subs(expression, new_expression)
103 if not isinstance(self.eigenvalues[i, d], numbers.Number):
104 self.eigenvalues[i, d] = self.eigenvalues[i, d].subs(
105 expression, new_expression
106 )
107
108 def implementation_of_flux(self, invoke_evalf_before_output=False):
109 """
110 Return implementation for flux along one coordinate axis (d) as C code.
111 invoke_evalf_before_output: boolean
112 If your expression is a symbolic expression (default) then we use evalf
113 before we pipe it into the output. If your expression is something numeric,
114 then evalf will fail (as it is not defined for scalar quantities).
115 """
116 result = "assertion(normal >= 0);\n"
117 result += "assertion(normal < Dimensions);\n"
119 is_cell_mapping=False
120 )
121 result += "switch( normal ) {\n"
122 for d in range(0, self.dimensionsdimensions):
123 result += " case " + str(d) + ":\n"
124 for i in range(0, self.unknownsunknowns):
125 if invoke_evalf_before_output:
126 result += sympy.printing.cxxcode(
127 self.F[i, d].evalf(), assign_to="F[" + str(i) + "]"
128 )
129 else:
130 result += (
131 "F["
132 + str(i)
133 + "] = "
134 + sympy.printing.cxxcode(self.F[i, d])
135 + ";"
136 )
137 result += "\n"
138 result += "assertion(!std::isnan(F[" + str(i) + "]));\n"
139 result += "nonCriticalAssertion(F[" + str(i) + "] == F[" + str(i) + "]);\n"
140 result += " break;\n"
141 result += "}\n"
142 return result
143
144 def implementation_of_ncp(self, invoke_evalf_before_output=False):
145 """
146 Return implementation for nonconservative product (ncp) along one coordinate axis (d) as C code.
147 invoke_evalf_before_output: boolean
148 If your expression is a symbolic expression (default) then we use evalf
149 before we pipe it into the output. If your expression is something numeric,
150 then evalf will fail (as it is not defined for scalar quantities).
151 """
152 result = "assertion(normal >= 0);\n"
153 result += "assertion(normal < Dimensions);\n"
155 is_cell_mapping=False, has_delta=True
156 )
157 result += "switch( normal ) {\n"
158 for d in range(0, self.dimensionsdimensions):
159 result += " case " + str(d) + ":\n"
160 for i in range(0, self.unknownsunknowns):
161 if invoke_evalf_before_output:
162 result += sympy.printing.cxxcode(
163 self.ncp[i, d].evalf(), assign_to="BTimesDeltaQ[" + str(i) + "]"
164 )
165 else:
166 result += (
167 "BTimesDeltaQ["
168 + str(i)
169 + "] = "
170 + sympy.printing.cxxcode(self.ncp[i, d])
171 + ";"
172 )
173 result += "\n"
174 result += " break;\n"
175 result += "}\n"
176 return result
177
178 def implementation_of_eigenvalues(self, invoke_evalf_before_output=False):
179 """
180 Return eigenvalues
181
182 This yields a set of eigenvalues, i.e., one per unknown. For many solvers such as
183 Rusanov, you need only one.
184
185 d: int
186 The axis along which we wanna have the eigenvalues
187
188 invoke_evalf_before_output: boolean
189 If your expression is a symbolic expression (default) then we use evalf
190 before we pipe it into the output. If your expression is something numeric,
191 then evalf will fail (as it is not defined for scalar quantities).
192 """
193 result = "assertion(normal >= 0);\n"
194 result += "assertion(normal < Dimensions);\n"
196 is_cell_mapping=False
197 )
198 result += "switch (normal) {\n"
199 for d in range(0, self.dimensionsdimensions):
200 result += " case " + str(d) + ":\n"
201 for i in range(0, self.unknownsunknowns):
202 if invoke_evalf_before_output:
203 result += sympy.printing.cxxcode(
204 self.eigenvalues[i, d].evalf(),
205 assign_to="lambda[" + str(i) + "]",
206 )
207 else:
208 result += (
209 "lambda["
210 + str(i)
211 + "] = "
212 + sympy.printing.cxxcode(self.eigenvalues[i, d])
213 + ";"
214 )
215 result += "\n"
216 result += " break;\n"
217 result += "}\n"
218 return result
219
221 self, invoke_evalf_before_output=False, use_absolute_values=True
222 ):
223 """
224 Return maximum eigenvalue
225 """
226 result = "assertion(normal >= 0);\n"
227 result += "assertion(normal < Dimensions);\n"
229 is_cell_mapping=False
230 )
231 result += "double lambda[" + str(self.unknownsunknowns) + "];\n"
232 result += "switch (normal) {\n"
233 for d in range(0, self.dimensionsdimensions):
234 result += " case " + str(d) + ":\n"
235 for i in range(0, self.unknownsunknowns):
236 if invoke_evalf_before_output:
237 result += sympy.printing.cxxcode(
238 self.eigenvalues[i, d].evalf(),
239 assign_to="lambda[" + str(i) + "]",
240 )
241 else:
242 result += (
243 "lambda["
244 + str(i)
245 + "] = "
246 + sympy.printing.cxxcode(self.eigenvalues[i, d])
247 + ";"
248 )
249 result += "\n"
250 result += "assertion(!std::isnan(lambda[" + str(i) + "]));\n"
251 result += " break;\n"
252 result += "}\n"
253 result += "double result = 0.0;\n"
254 for i in range(0, self.unknownsunknowns):
255 if use_absolute_values:
256 result += (
257 "result = std::max( result, std::abs(lambda[" + str(i) + "]) );\n"
258 )
259 else:
260 result += "result = std::max( result, lambda[" + str(i) + "] );\n"
261 result += "return result;\n"
262 return result
263
264 def implementation_of_sources(self, invoke_evalf_before_output=False):
265 """
266 Return implementation for sources as C code.
267 """
268 result = "assertion(normal >= 0);\n"
269 result += "assertion(normal < Dimensions);\n"
271 is_cell_mapping=True
272 )
273 for i in range(0, self.unknownsunknowns):
274 if invoke_evalf_before_output:
275 result += sympy.printing.cxxcode(
276 self.sources[i].evalf(), assign_to="S[" + str(i) + "]"
277 )
278 else:
279 result += (
280 "S["
281 + str(i)
282 + "] = "
283 + sympy.printing.cxxcode(self.sources[i])
284 + ";"
285 )
286 result += "\n"
287 return result
Helper class to model a hyperbolic PDE in first-order conservative formulation.
implementation_of_ncp(self, invoke_evalf_before_output=False)
Return implementation for nonconservative product (ncp) along one coordinate axis (d) as C code.
implementation_of_max_eigenvalue(self, invoke_evalf_before_output=False, use_absolute_values=True)
Return maximum eigenvalue.
implementation_of_flux(self, invoke_evalf_before_output=False)
Return implementation for flux along one coordinate axis (d) as C code.
implementation_of_sources(self, invoke_evalf_before_output=False)
Return implementation for sources as C code.
substitute_expression(self, expression, new_expression)
Usually used to set a symbolic variable (expression) to a value (new_expression).
_implementation_of_mapping_onto_named_quantities(self, is_cell_mapping=True, is_boundary=False, has_delta=False)
Return the C code that maps the quantities from Q onto properly labelled quantities.
Definition PDE.py:34