Peano
Loading...
Searching...
No Matches
RusanovWithOptimisticExplicitEuler.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 re
5
6import peano4.api
7import peano4.api.datamodel
8import peano4.api.dastgen2
9import toolbox.blockstructured.api
10import toolbox.loadbalancing.api.actionsets
11import exahype2.api
12from exahype2.api import Target
13import exahype2.api.dslhype as dslhype
14import dastgen2
15
16from ..Solver import Solver
17from ..blockstructured.BoundaryConditions import get_access_to_persistent_face_data
18from ..blockstructured.create_dynamic_mesh_refinement_action_set import (
19 create_dynamic_mesh_refinement_action_set,
20)
21from .actionsets import ApplyDirichletBoundaryConditions
22from .actionsets import ProjectPatchOntoFaces
23from .actionsets import RollOverUpdatedFaceWithSommerfeldBoundaryConditions
24from .actionsets import UpdateCell
25from .actionsets import ComputeGlobalCFLCondition
26from .actionsets import UpdateCellInterior
27from .actionsets import UpdateCellShell
28from .actionsets import InsertPatchInterior
29from .actionsets import DiscardPatchInterior
30from .actionsets import SpawnTaskToUpdateCellInterior
31from .actionsets import ReadjustTimeStepSizeForOptimisticGlobalTimestepping
32
33
34
66
67
68
69def timeStepWithRusanov2DOnPatchInterior(patchData: "::exahype2::CellData<double, double>&") -> "void":
70 N = Integer("patchData.numberOfCells")
71 t = DataBlock([[0, N]], "patchData.t", "double")
72 dt = DataBlock([[0, N]], "patchData.dt", "double")
73 h = DataBlock([[0, N]], "patchData.cellSize", "tarch::la::Vector<2, double>")
74 x = DataBlock([[0, N]], "patchData.cellCentre", "tarch::la::Vector<2, double>")
75 maxEigenvalues = DataBlock([[0, N]], "patchData.maxEigenvalue", "double")
76
77 input = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV], [0, N]], "patchData.QIn")
78 output = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [H, NV-H], [H, NV-H], [0, N]], "patchData.QOut")
79
80 output[:, H:NV-H, H:NV-H] = input[:, H:NV-H, H:NV-H]
81
82 if volumetricSource:
83 source = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H, NV-H], [0, N]], eval("volumetricSource", input[:, H:NV-H, H:NV-H], x, t))
84 output[0:NumberOfUnknowns, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H] + dt * source
85
86 maxEigenvaluesX = DataBlock([[0, 1], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, H-1:NV-H+1, H:NV-H], x, 0))
87 maxEigenvaluesY = DataBlock([[0, 1], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, H:NV-H, H-1:NV-H+1], x, 1))
88
89 lambdaLeft = max(maxEigenvaluesX[0:1, H-1:NV-H-1, H:NV-H], maxEigenvaluesX[0:1, H:NV-H, H:NV-H])
90 lambdaRight = max(maxEigenvaluesX[0:1, H:NV-H, H:NV-H], maxEigenvaluesX[0:1, H+1:NV-H+1, H:NV-H])
91 lambdaBottom = max(maxEigenvaluesY[0:1, H:NV-H, H-1:NV-H-1], maxEigenvaluesY[0:1, H:NV-H, H:NV-H])
92 lambdaTop = max(maxEigenvaluesY[0:1, H:NV-H, H:NV-H], maxEigenvaluesY[0:1, H:NV-H, H+1:NV-H+1])
93
94 fluxLeft = -1.0 * lambdaLeft * (input[0:NumberOfUnknowns, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, H-1:NV-H-1, H:NV-H])
95 fluxRight = -1.0 * lambdaRight * (input[0:NumberOfUnknowns, H+1:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns, H:NV-H, H:NV-H])
96 output[0:NumberOfUnknowns, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H] + dt / h(0) * (fluxLeft - fluxRight) * NV * 0.5
97
98 fluxBottom = -1.0 * lambdaBottom * (input[0:NumberOfUnknowns, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, H:NV-H, H-1:NV-H-1])
99 fluxTop = -1.0 * lambdaTop * (input[0:NumberOfUnknowns, H:NV-H, H+1:NV-H+1] - input[0:NumberOfUnknowns, H:NV-H, H:NV-H])
100 output[0:NumberOfUnknowns, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H] + dt / h(1) * (fluxBottom - fluxTop) * NV * 0.5
101
102 if flux:
103 flux_x = DataBlock([[0, NumberOfUnknowns], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("flux", input[:, H-1:NV-H+1, H:NV-H], x, 0))
104 flux_y = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("flux", input[:, H:NV-H, H-1:NV-H+1], x, 1))
105 output[0:NumberOfUnknowns, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_x[:, H-1:NV-H-1, H:NV-H] - flux_x[:, H+1:NV-H+1, H:NV-H])
106 output[0:NumberOfUnknowns, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_y[:, H:NV-H, H-1:NV-H-1] - flux_y[:, H:NV-H, H+1:NV-H+1])
107
108 if nonconservativeProduct:
109 average_x = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H-1:NV-H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H+1, H:NV-H])
110 delta_x = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H-1:NV-H, H:NV-H]
111 nonconservative_x = DataBlock([[0, NumberOfUnknowns], [H, NV-H+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_x, delta_x, x, h, t, dt, 0))
112 output[0:NumberOfUnknowns, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_x[:, H:NV-H, H:NV-H] - nonconservative_x[:, H+1:NV-H+1, H:NV-H])
113
114 average_y = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H+1])
115 delta_y = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H-1:NV-H]
116 nonconservative_y = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_y, delta_y, x, h, t, dt, 1))
117 output[0:NumberOfUnknowns, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_y[:, H:NV-H, H:NV-H] - nonconservative_y[:, H:NV-H, H+1:NV-H+1])
118
119 if maxEigenvalue:
120 maxEigenvaluesX2 = DataBlock([[0, 1], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, H:NV-H, H:NV-H], x, 0))
121 maxEigenvaluesY2 = DataBlock([[0, 1], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, H:NV-H, H:NV-H], x, 1))
122 maxEigenvaluesPerVolume = max(maxEigenvaluesX2, maxEigenvaluesY2)
123 maxEigenvalues = max(maxEigenvaluesPerVolume)
124
125
126def timeStepWithRusanov2DOnShell(patchData: "::exahype2::CellData<double, double>&") -> "void":
127 N = Integer("patchData.numberOfCells")
128 t = DataBlock([[0, N]], "patchData.t", "double")
129 dt = DataBlock([[0, N]], "patchData.dt", "double")
130 h = DataBlock([[0, N]], "patchData.cellSize", "tarch::la::Vector<2, double>")
131 x = DataBlock([[0, N]], "patchData.cellCentre", "tarch::la::Vector<2, double>")
132 maxEigenvalues = DataBlock([[0, N]], "patchData.maxEigenvalue", "double")
133
134 input = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [-H, NV + H], [-H, NV + H], [0, N]], "patchData.QIn")
135 output = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV], [0, N]], "patchData.QOut")
136
137 # Shell part YLow: [0,NV) x [0,H)
138 output[:, 0:NV, 0:H] = input[:, 0:NV, 0:H]
139
140 if volumetricSource:
141 sourceYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, H], [0, N]], eval("volumetricSource", input[:, 0:NV, 0:H], x, t))
142 output[0:NumberOfUnknowns, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:H] + dt * sourceYLow
143
144 maxEigenvaluesXYLow = DataBlock([[0, 1], [0-1, NV+1], [0, H], [0, N]], eval("maxEigenvalue", input[:, 0-1:NV+1, 0:H], x, 0))
145 maxEigenvaluesYYLow = DataBlock([[0, 1], [0, NV], [0-1, H+1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0-1:H+1], x, 1))
146
147 lambdaLeftYLow = max(maxEigenvaluesXYLow[0:1, 0-1:NV-1, 0:H], maxEigenvaluesXYLow[0:1, 0:NV, 0:H])
148 lambdaRightYLow = max(maxEigenvaluesXYLow[0:1, 0:NV, 0:H], maxEigenvaluesXYLow[0:1, 0+1:NV+1, 0:H])
149 lambdaBottomYLow = max(maxEigenvaluesYYLow[0:1, 0:NV, 0-1:H-1], maxEigenvaluesYYLow[0:1, 0:NV, 0:H])
150 lambdaTopYLow = max(maxEigenvaluesYYLow[0:1, 0:NV, 0:H], maxEigenvaluesYYLow[0:1, 0:NV, 0+1:H+1])
151
152 fluxLeftYLow = -1.0 * lambdaLeftYLow * (input[0:NumberOfUnknowns, 0:NV, 0:H] - input[0:NumberOfUnknowns, 0-1:NV-1, 0:H])
153 fluxRightYLow = -1.0 * lambdaRightYLow * (input[0:NumberOfUnknowns, 0+1:NV+1, 0:H] - input[0:NumberOfUnknowns, 0:NV, 0:H])
154 output[0:NumberOfUnknowns, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:H] + dt / h(0) * (fluxLeftYLow - fluxRightYLow) * NV * 0.5
155
156 fluxBottomYLow = -1.0 * lambdaBottomYLow * (input[0:NumberOfUnknowns, 0:NV, 0:H] - input[0:NumberOfUnknowns, 0:NV, 0-1:H-1])
157 fluxTopYLow = -1.0 * lambdaTopYLow * (input[0:NumberOfUnknowns, 0:NV, 0+1:H+1] - input[0:NumberOfUnknowns, 0:NV, 0:H])
158 output[0:NumberOfUnknowns, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:H] + dt / h(1) * (fluxBottomYLow - fluxTopYLow) * NV * 0.5
159
160 if flux:
161 flux_xYLow = DataBlock([[0, NumberOfUnknowns], [0-1, NV+1], [0, H], [0, N]], eval("flux", input[:, 0-1:NV+1, 0:H], x, 0))
162 flux_yYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0-1, H+1], [0, N]], eval("flux", input[:, 0:NV, 0-1:H+1], x, 1))
163 output[0:NumberOfUnknowns, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:H] + 0.5 * NV * dt / h(0) * (flux_xYLow[:, 0-1:NV-1, 0:H] - flux_xYLow[:, 0+1:NV+1, 0:H])
164 output[0:NumberOfUnknowns, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:H] + 0.5 * NV * dt / h(1) * (flux_yYLow[:, 0:NV, 0-1:H-1] - flux_yYLow[:, 0:NV, 0+1:H+1])
165
166 if nonconservativeProduct:
167 average_xYLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:H])
168 delta_xYLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:H]
169 nonconservative_xYLow = DataBlock([[0, NumberOfUnknowns], [0, NV+1], [0, H], [0, N]], eval("nonconservativeProduct", average_xYLow, delta_xYLow, x, h, t, dt, 0))
170 output[0:NumberOfUnknowns, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:H] + 0.5 * NV * dt / h(0) * (-nonconservative_xYLow[:, 0:NV, 0:H] - nonconservative_xYLow[:, 0+1:NV+1, 0:H])
171
172 average_yYLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H+1])
173 delta_yYLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:H]
174 nonconservative_yYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, H+1], [0, N]], eval("nonconservativeProduct", average_yYLow, delta_yYLow, x, h, t, dt, 1))
175 output[0:NumberOfUnknowns, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:H] + 0.5 * NV * dt / h(1) * (-nonconservative_yYLow[:, 0:NV, 0:H] - nonconservative_yYLow[:, 0:NV, 0+1:H+1])
176
177 # Shell part YHigh: [0,NV) x [NV-H,NV)
178 output[:, 0:NV, NV-H:NV] = input[:, 0:NV, NV-H:NV]
179
180 if volumetricSource:
181 sourceYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H, NV], [0, N]], eval("volumetricSource", input[:, 0:NV, NV-H:NV], x, t))
182 output[0:NumberOfUnknowns, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV] + dt * sourceYHigh
183
184 maxEigenvaluesXYHigh = DataBlock([[0, 1], [0-1, NV+1], [NV-H, NV], [0, N]], eval("maxEigenvalue", input[:, 0-1:NV+1, NV-H:NV], x, 0))
185 maxEigenvaluesYYHigh = DataBlock([[0, 1], [0, NV], [NV-H-1, NV+1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, NV-H-1:NV+1], x, 1))
186
187 lambdaLeftYHigh = max(maxEigenvaluesXYHigh[0:1, 0-1:NV-1, NV-H:NV], maxEigenvaluesXYHigh[0:1, 0:NV, NV-H:NV])
188 lambdaRightYHigh = max(maxEigenvaluesXYHigh[0:1, 0:NV, NV-H:NV], maxEigenvaluesXYHigh[0:1, 0+1:NV+1, NV-H:NV])
189 lambdaBottomYHigh = max(maxEigenvaluesYYHigh[0:1, 0:NV, NV-H-1:NV-1], maxEigenvaluesYYHigh[0:1, 0:NV, NV-H:NV])
190 lambdaTopYHigh = max(maxEigenvaluesYYHigh[0:1, 0:NV, NV-H:NV], maxEigenvaluesYYHigh[0:1, 0:NV, NV-H+1:NV+1])
191
192 fluxLeftYHigh = -1.0 * lambdaLeftYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H:NV] - input[0:NumberOfUnknowns, 0-1:NV-1, NV-H:NV])
193 fluxRightYHigh = -1.0 * lambdaRightYHigh * (input[0:NumberOfUnknowns, 0+1:NV+1, NV-H:NV] - input[0:NumberOfUnknowns, 0:NV, NV-H:NV])
194 output[0:NumberOfUnknowns, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV] + dt / h(0) * (fluxLeftYHigh - fluxRightYHigh) * NV * 0.5
195
196 fluxBottomYHigh = -1.0 * lambdaBottomYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H:NV] - input[0:NumberOfUnknowns, 0:NV, NV-H-1:NV-1])
197 fluxTopYHigh = -1.0 * lambdaTopYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H+1:NV+1] - input[0:NumberOfUnknowns, 0:NV, NV-H:NV])
198 output[0:NumberOfUnknowns, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV] + dt / h(1) * (fluxBottomYHigh - fluxTopYHigh) * NV * 0.5
199
200 if flux:
201 flux_xYHigh = DataBlock([[0, NumberOfUnknowns], [0-1, NV+1], [NV-H, NV], [0, N]], eval("flux", input[:, 0-1:NV+1, NV-H:NV], x, 0))
202 flux_yYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H-1, NV+1], [0, N]], eval("flux", input[:, 0:NV, NV-H-1:NV+1], x, 1))
203 output[0:NumberOfUnknowns, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(0) * (flux_xYHigh[:, 0-1:NV-1, NV-H:NV] - flux_xYHigh[:, 0+1:NV+1, NV-H:NV])
204 output[0:NumberOfUnknowns, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(1) * (flux_yYHigh[:, 0:NV, NV-H-1:NV-1] - flux_yYHigh[:, 0:NV, NV-H+1:NV+1])
205
206 if nonconservativeProduct:
207 average_xYHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, NV-H:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, NV-H:NV])
208 delta_xYHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, NV-H:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, NV-H:NV]
209 nonconservative_xYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV+1], [NV-H, NV], [0, N]], eval("nonconservativeProduct", average_xYHigh, delta_xYHigh, x, h, t, dt, 0))
210 output[0:NumberOfUnknowns, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(0) * (-nonconservative_xYHigh[:, 0:NV, NV-H:NV] - nonconservative_xYHigh[:, 0+1:NV+1, NV-H:NV])
211
212 average_yYHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H-1:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV+1])
213 delta_yYHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H-1:NV]
214 nonconservative_yYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H, NV+1], [0, N]], eval("nonconservativeProduct", average_yYHigh, delta_yYHigh, x, h, t, dt, 1))
215 output[0:NumberOfUnknowns, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(1) * (-nonconservative_yYHigh[:, 0:NV, NV-H:NV] - nonconservative_yYHigh[:, 0:NV, NV-H+1:NV+1])
216
217 # Shell part XLow: [0,H) x [H,NV-H)
218 output[:, 0:H, H:NV-H] = input[:, 0:H, H:NV-H]
219
220 if volumetricSource:
221 sourceXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H, NV-H], [0, N]], eval("volumetricSource", input[:, 0:H, H:NV-H], x, t))
222 output[0:NumberOfUnknowns, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H] + dt * sourceXLow
223
224 maxEigenvaluesXXLow = DataBlock([[0, 1], [0-1, H+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, 0-1:H+1, H:NV-H], x, 0))
225 maxEigenvaluesYXLow = DataBlock([[0, 1], [0, H], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, 0:H, H-1:NV-H+1], x, 1))
226
227 lambdaLeftXLow = max(maxEigenvaluesXXLow[0:1, 0-1:H-1, H:NV-H], maxEigenvaluesXXLow[0:1, 0:H, H:NV-H])
228 lambdaRightXLow = max(maxEigenvaluesXXLow[0:1, 0:H, H:NV-H], maxEigenvaluesXXLow[0:1, 0+1:H+1, H:NV-H])
229 lambdaBottomXLow = max(maxEigenvaluesYXLow[0:1, 0:H, H-1:NV-H-1], maxEigenvaluesYXLow[0:1, 0:H, H:NV-H])
230 lambdaTopXLow = max(maxEigenvaluesYXLow[0:1, 0:H, H:NV-H], maxEigenvaluesYXLow[0:1, 0:H, H+1:NV-H+1])
231
232 fluxLeftXLow = -1.0 * lambdaLeftXLow * (input[0:NumberOfUnknowns, 0:H, H:NV-H] - input[0:NumberOfUnknowns, 0-1:H-1, H:NV-H])
233 fluxRightXLow = -1.0 * lambdaRightXLow * (input[0:NumberOfUnknowns, 0+1:H+1, H:NV-H] - input[0:NumberOfUnknowns, 0:H, H:NV-H])
234 output[0:NumberOfUnknowns, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H] + dt / h(0) * (fluxLeftXLow - fluxRightXLow) * NV * 0.5
235
236 fluxBottomXLow = -1.0 * lambdaBottomXLow * (input[0:NumberOfUnknowns, 0:H, H:NV-H] - input[0:NumberOfUnknowns, 0:H, H-1:NV-H-1])
237 fluxTopXLow = -1.0 * lambdaTopXLow * (input[0:NumberOfUnknowns, 0:H, H+1:NV-H+1] - input[0:NumberOfUnknowns, 0:H, H:NV-H])
238 output[0:NumberOfUnknowns, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H] + dt / h(1) * (fluxBottomXLow - fluxTopXLow) * NV * 0.5
239
240 if flux:
241 flux_xXLow = DataBlock([[0, NumberOfUnknowns], [0-1, H+1], [H, NV-H], [0, N]], eval("flux", input[:, 0-1:H+1, H:NV-H], x, 0))
242 flux_yXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H-1, NV-H+1], [0, N]], eval("flux", input[:, 0:H, H-1:NV-H+1], x, 1))
243 output[0:NumberOfUnknowns, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_xXLow[:, 0-1:H-1, H:NV-H] - flux_xXLow[:, 0+1:H+1, H:NV-H])
244 output[0:NumberOfUnknowns, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_yXLow[:, 0:H, H-1:NV-H-1] - flux_yXLow[:, 0:H, H+1:NV-H+1])
245
246 if nonconservativeProduct:
247 average_xXLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H+1, H:NV-H])
248 delta_xXLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:H, H:NV-H]
249 nonconservative_xXLow = DataBlock([[0, NumberOfUnknowns], [0, H+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_xXLow, delta_xXLow, x, h, t, dt, 0))
250 output[0:NumberOfUnknowns, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_xXLow[:, 0:H, H:NV-H] - nonconservative_xXLow[:, 0+1:H+1, H:NV-H])
251
252 average_yXLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H+1])
253 delta_yXLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H-1:NV-H]
254 nonconservative_yXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_yXLow, delta_yXLow, x, h, t, dt, 1))
255 output[0:NumberOfUnknowns, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_yXLow[:, 0:H, H:NV-H] - nonconservative_yXLow[:, 0:H, H+1:NV-H+1])
256
257 # Shell part XHigh: [NV-H,NV) x [H,NV-H)
258 output[:, NV-H:NV, H:NV-H] = input[:, NV-H:NV, H:NV-H]
259
260 if volumetricSource:
261 sourceXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H, NV-H], [0, N]], eval("volumetricSource", input[:, NV-H:NV, H:NV-H], x, t))
262 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] + dt * sourceXHigh
263
264 maxEigenvaluesXXHigh = DataBlock([[0, 1], [NV-H-1, NV+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, NV-H-1:NV+1, H:NV-H], x, 0))
265 maxEigenvaluesYXHigh = DataBlock([[0, 1], [NV-H, NV], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, NV-H:NV, H-1:NV-H+1], x, 1))
266
267 lambdaLeftXHigh = max(maxEigenvaluesXXHigh[0:1, NV-H-1:NV-1, H:NV-H], maxEigenvaluesXXHigh[0:1, NV-H:NV, H:NV-H])
268 lambdaRightXHigh = max(maxEigenvaluesXXHigh[0:1, NV-H:NV, H:NV-H], maxEigenvaluesXXHigh[0:1, NV-H+1:NV+1, H:NV-H])
269 lambdaBottomXHigh = max(maxEigenvaluesYXHigh[0:1, NV-H:NV, H-1:NV-H-1], maxEigenvaluesYXHigh[0:1, NV-H:NV, H:NV-H])
270 lambdaTopXHigh = max(maxEigenvaluesYXHigh[0:1, NV-H:NV, H:NV-H], maxEigenvaluesYXHigh[0:1, NV-H:NV, H+1:NV-H+1])
271
272 fluxLeftXHigh = -1.0 * lambdaLeftXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H:NV-H] - input[0:NumberOfUnknowns, NV-H-1:NV-1, H:NV-H])
273 fluxRightXHigh = -1.0 * lambdaRightXHigh * (input[0:NumberOfUnknowns, NV-H+1:NV+1, H:NV-H] - input[0:NumberOfUnknowns, NV-H:NV, H:NV-H])
274 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] + dt / h(0) * (fluxLeftXHigh - fluxRightXHigh) * NV * 0.5
275
276 fluxBottomXHigh = -1.0 * lambdaBottomXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H:NV-H] - input[0:NumberOfUnknowns, NV-H:NV, H-1:NV-H-1])
277 fluxTopXHigh = -1.0 * lambdaTopXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H+1:NV-H+1] - input[0:NumberOfUnknowns, NV-H:NV, H:NV-H])
278 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] + dt / h(1) * (fluxBottomXHigh - fluxTopXHigh) * NV * 0.5
279
280 if flux:
281 flux_xXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H-1, NV+1], [H, NV-H], [0, N]], eval("flux", input[:, NV-H-1:NV+1, H:NV-H], x, 0))
282 flux_yXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H-1, NV-H+1], [0, N]], eval("flux", input[:, NV-H:NV, H-1:NV-H+1], x, 1))
283 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_xXHigh[:, NV-H-1:NV-1, H:NV-H] - flux_xXHigh[:, NV-H+1:NV+1, H:NV-H])
284 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_yXHigh[:, NV-H:NV, H-1:NV-H-1] - flux_yXHigh[:, NV-H:NV, H+1:NV-H+1])
285
286 if nonconservativeProduct:
287 average_xXHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H-1:NV, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV+1, H:NV-H])
288 delta_xXHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H-1:NV, H:NV-H]
289 nonconservative_xXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_xXHigh, delta_xXHigh, x, h, t, dt, 0))
290 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_xXHigh[:, NV-H:NV, H:NV-H] - nonconservative_xXHigh[:, NV-H+1:NV+1, H:NV-H])
291
292 average_yXHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H+1])
293 delta_yXHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H-1:NV-H]
294 nonconservative_yXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_yXHigh, delta_yXHigh, x, h, t, dt, 1))
295 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_yXHigh[:, NV-H:NV, H:NV-H] - nonconservative_yXHigh[:, NV-H:NV, H+1:NV-H+1])
296
297 if maxEigenvalue:
298 # Reduce per shell part, as max() over one DataBlock overwrites its target
299 maxEigenvaluesX2YLow = DataBlock([[0, 1], [0, NV], [0, H], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:H], x, 0))
300 maxEigenvaluesY2YLow = DataBlock([[0, 1], [0, NV], [0, H], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:H], x, 1))
301 maxEigenvaluesPerVolumeYLow = max(maxEigenvaluesX2YLow, maxEigenvaluesY2YLow)
302 maxEigenvalueYLow = DataBlock([[0, N]], 0.0)
303 maxEigenvalueYLow = max(maxEigenvaluesPerVolumeYLow)
304 maxEigenvaluesX2YHigh = DataBlock([[0, 1], [0, NV], [0, NV-(NV-H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, NV-H:NV], x, 0))
305 maxEigenvaluesY2YHigh = DataBlock([[0, 1], [0, NV], [0, NV-(NV-H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, NV-H:NV], x, 1))
306 maxEigenvaluesPerVolumeYHigh = max(maxEigenvaluesX2YHigh, maxEigenvaluesY2YHigh)
307 maxEigenvalueYHigh = DataBlock([[0, N]], 0.0)
308 maxEigenvalueYHigh = max(maxEigenvaluesPerVolumeYHigh)
309 maxEigenvaluesX2XLow = DataBlock([[0, 1], [0, H], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:H, H:NV-H], x, 0))
310 maxEigenvaluesY2XLow = DataBlock([[0, 1], [0, H], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:H, H:NV-H], x, 1))
311 maxEigenvaluesPerVolumeXLow = max(maxEigenvaluesX2XLow, maxEigenvaluesY2XLow)
312 maxEigenvalueXLow = DataBlock([[0, N]], 0.0)
313 maxEigenvalueXLow = max(maxEigenvaluesPerVolumeXLow)
314 maxEigenvaluesX2XHigh = DataBlock([[0, 1], [0, NV-(NV-H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, NV-H:NV, H:NV-H], x, 0))
315 maxEigenvaluesY2XHigh = DataBlock([[0, 1], [0, NV-(NV-H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, NV-H:NV, H:NV-H], x, 1))
316 maxEigenvaluesPerVolumeXHigh = max(maxEigenvaluesX2XHigh, maxEigenvaluesY2XHigh)
317 maxEigenvalueXHigh = DataBlock([[0, N]], 0.0)
318 maxEigenvalueXHigh = max(maxEigenvaluesPerVolumeXHigh)
319 maxEigenvalues = max(maxEigenvalueYLow, maxEigenvalueYHigh)
320 maxEigenvalues = max(maxEigenvalues, maxEigenvalueXLow)
321 maxEigenvalues = max(maxEigenvalues, maxEigenvalueXHigh)
322
323
324def timeStepWithRusanov3DOnPatchInterior(patchData: "::exahype2::CellData<double, double>&") -> "void":
325 N = Integer("patchData.numberOfCells")
326 t = DataBlock([[0, N]], "patchData.t", "double")
327 dt = DataBlock([[0, N]], "patchData.dt", "double")
328 h = DataBlock([[0, N]], "patchData.cellSize", "tarch::la::Vector<3, double>")
329 x = DataBlock([[0, N]], "patchData.cellCentre", "tarch::la::Vector<3, double>")
330 maxEigenvalues = DataBlock([[0, N]], "patchData.maxEigenvalue", "double")
331
332 input = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV], [0, NV], [0, N]], "patchData.QIn")
333 output = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [H, NV-H], [H, NV-H], [H, NV-H], [0, N]], "patchData.QOut")
334
335 output[:, H:NV-H, H:NV-H, H:NV-H] = input[:, H:NV-H, H:NV-H, H:NV-H]
336
337 if volumetricSource:
338 source = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H, NV-H], [H, NV-H], [0, N]], eval("volumetricSource", input[:, H:NV-H, H:NV-H, H:NV-H], x, t))
339 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + dt * source
340
341 maxEigenvaluesX = DataBlock([[0, 1], [H-1, NV-H+1], [H, NV-H], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, H-1:NV-H+1, H:NV-H, H:NV-H], x, 0))
342 maxEigenvaluesY = DataBlock([[0, 1], [H, NV-H], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, H:NV-H, H-1:NV-H+1, H:NV-H], x, 1))
343 maxEigenvaluesZ = DataBlock([[0, 1], [H, NV-H], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, H:NV-H, H:NV-H, H-1:NV-H+1], x, 2))
344
345 lambdaLeft = max(maxEigenvaluesX[0:1, H-1:NV-H-1, H:NV-H, H:NV-H], maxEigenvaluesX[0:1, H:NV-H, H:NV-H, H:NV-H])
346 lambdaRight = max(maxEigenvaluesX[0:1, H:NV-H, H:NV-H, H:NV-H], maxEigenvaluesX[0:1, H+1:NV-H+1, H:NV-H, H:NV-H])
347 lambdaBottom = max(maxEigenvaluesY[0:1, H:NV-H, H-1:NV-H-1, H:NV-H], maxEigenvaluesY[0:1, H:NV-H, H:NV-H, H:NV-H])
348 lambdaTop = max(maxEigenvaluesY[0:1, H:NV-H, H:NV-H, H:NV-H], maxEigenvaluesY[0:1, H:NV-H, H+1:NV-H+1, H:NV-H])
349 lambdaNear = max(maxEigenvaluesZ[0:1, H:NV-H, H:NV-H, H-1:NV-H-1], maxEigenvaluesZ[0:1, H:NV-H, H:NV-H, H:NV-H])
350 lambdaFar = max(maxEigenvaluesZ[0:1, H:NV-H, H:NV-H, H:NV-H], maxEigenvaluesZ[0:1, H:NV-H, H:NV-H, H+1:NV-H+1])
351
352 fluxLeft = -1.0 * lambdaLeft * (input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, H-1:NV-H-1, H:NV-H, H:NV-H])
353 fluxRight = -1.0 * lambdaRight * (input[0:NumberOfUnknowns, H+1:NV-H+1, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H])
354 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + dt / h(0) * (fluxLeft - fluxRight) * NV * 0.5
355
356 fluxBottom = -1.0 * lambdaBottom * (input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, H:NV-H, H-1:NV-H-1, H:NV-H])
357 fluxTop = -1.0 * lambdaTop * (input[0:NumberOfUnknowns, H:NV-H, H+1:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H])
358 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + dt / h(1) * (fluxBottom - fluxTop) * NV * 0.5
359
360 fluxNear = -1.0 * lambdaNear * (input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H-1:NV-H-1])
361 fluxFar = -1.0 * lambdaFar * (input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H+1:NV-H+1] - input[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H])
362 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + dt / h(2) * (fluxNear - fluxFar) * NV * 0.5
363
364 if flux:
365 flux_x = DataBlock([[0, NumberOfUnknowns], [H-1, NV-H+1], [H, NV-H], [H, NV-H], [0, N]], eval("flux", input[:, H-1:NV-H+1, H:NV-H, H:NV-H], x, 0))
366 flux_y = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("flux", input[:, H:NV-H, H-1:NV-H+1, H:NV-H], x, 1))
367 flux_z = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("flux", input[:, H:NV-H, H:NV-H, H-1:NV-H+1], x, 2))
368 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_x[:, H-1:NV-H-1, H:NV-H, H:NV-H] - flux_x[:, H+1:NV-H+1, H:NV-H, H:NV-H])
369 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_y[:, H:NV-H, H-1:NV-H-1, H:NV-H] - flux_y[:, H:NV-H, H+1:NV-H+1, H:NV-H])
370 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(2) * (flux_z[:, H:NV-H, H:NV-H, H-1:NV-H-1] - flux_z[:, H:NV-H, H:NV-H, H+1:NV-H+1])
371
372 if nonconservativeProduct:
373 average_x = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H-1:NV-H, H:NV-H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H+1, H:NV-H, H:NV-H])
374 delta_x = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H+1, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H-1:NV-H, H:NV-H, H:NV-H]
375 nonconservative_x = DataBlock([[0, NumberOfUnknowns], [H, NV-H+1], [H, NV-H], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_x, delta_x, x, h, t, dt, 0))
376 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_x[:, H:NV-H, H:NV-H, H:NV-H] - nonconservative_x[:, H+1:NV-H+1, H:NV-H, H:NV-H])
377
378 average_y = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H-1:NV-H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H+1, H:NV-H])
379 delta_y = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H-1:NV-H, H:NV-H]
380 nonconservative_y = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H, NV-H+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_y, delta_y, x, h, t, dt, 1))
381 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_y[:, H:NV-H, H:NV-H, H:NV-H] - nonconservative_y[:, H:NV-H, H+1:NV-H+1, H:NV-H])
382
383 average_z = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H, H:NV-H+1])
384 delta_z = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, H:NV-H, H:NV-H, H-1:NV-H]
385 nonconservative_z = DataBlock([[0, NumberOfUnknowns], [H, NV-H], [H, NV-H], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_z, delta_z, x, h, t, dt, 2))
386 output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, H:NV-H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(2) * (-nonconservative_z[:, H:NV-H, H:NV-H, H:NV-H] - nonconservative_z[:, H:NV-H, H:NV-H, H+1:NV-H+1])
387
388 if maxEigenvalue:
389 maxEigenvaluesX2 = DataBlock([[0, 1], [0, NV-H-(H)], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, H:NV-H, H:NV-H, H:NV-H], x, 0))
390 maxEigenvaluesY2 = DataBlock([[0, 1], [0, NV-H-(H)], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, H:NV-H, H:NV-H, H:NV-H], x, 1))
391 maxEigenvaluesZ2 = DataBlock([[0, 1], [0, NV-H-(H)], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, H:NV-H, H:NV-H, H:NV-H], x, 2))
392 maxEigenvaluesXY = max(maxEigenvaluesX2, maxEigenvaluesY2)
393 maxEigenvaluesPerVolume = max(maxEigenvaluesXY, maxEigenvaluesZ2)
394 maxEigenvalues = max(maxEigenvaluesPerVolume)
395
396
397def timeStepWithRusanov3DOnShell(patchData: "::exahype2::CellData<double, double>&") -> "void":
398 N = Integer("patchData.numberOfCells")
399 t = DataBlock([[0, N]], "patchData.t", "double")
400 dt = DataBlock([[0, N]], "patchData.dt", "double")
401 h = DataBlock([[0, N]], "patchData.cellSize", "tarch::la::Vector<3, double>")
402 x = DataBlock([[0, N]], "patchData.cellCentre", "tarch::la::Vector<3, double>")
403 maxEigenvalues = DataBlock([[0, N]], "patchData.maxEigenvalue", "double")
404
405 input = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [-H, NV + H], [-H, NV + H], [-H, NV + H], [0, N]], "patchData.QIn")
406 output = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV], [0, NV], [0, N]], "patchData.QOut")
407
408 # Shell part ZLow: [0,NV) x [0,NV) x [0,H)
409 output[:, 0:NV, 0:NV, 0:H] = input[:, 0:NV, 0:NV, 0:H]
410
411 if volumetricSource:
412 sourceZLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [0, H], [0, N]], eval("volumetricSource", input[:, 0:NV, 0:NV, 0:H], x, t))
413 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + dt * sourceZLow
414
415 maxEigenvaluesXZLow = DataBlock([[0, 1], [0-1, NV+1], [0, NV], [0, H], [0, N]], eval("maxEigenvalue", input[:, 0-1:NV+1, 0:NV, 0:H], x, 0))
416 maxEigenvaluesYZLow = DataBlock([[0, 1], [0, NV], [0-1, NV+1], [0, H], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0-1:NV+1, 0:H], x, 1))
417 maxEigenvaluesZZLow = DataBlock([[0, 1], [0, NV], [0, NV], [0-1, H+1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0:NV, 0-1:H+1], x, 2))
418
419 lambdaLeftZLow = max(maxEigenvaluesXZLow[0:1, 0-1:NV-1, 0:NV, 0:H], maxEigenvaluesXZLow[0:1, 0:NV, 0:NV, 0:H])
420 lambdaRightZLow = max(maxEigenvaluesXZLow[0:1, 0:NV, 0:NV, 0:H], maxEigenvaluesXZLow[0:1, 0+1:NV+1, 0:NV, 0:H])
421 lambdaBottomZLow = max(maxEigenvaluesYZLow[0:1, 0:NV, 0-1:NV-1, 0:H], maxEigenvaluesYZLow[0:1, 0:NV, 0:NV, 0:H])
422 lambdaTopZLow = max(maxEigenvaluesYZLow[0:1, 0:NV, 0:NV, 0:H], maxEigenvaluesYZLow[0:1, 0:NV, 0+1:NV+1, 0:H])
423 lambdaNearZLow = max(maxEigenvaluesZZLow[0:1, 0:NV, 0:NV, 0-1:H-1], maxEigenvaluesZZLow[0:1, 0:NV, 0:NV, 0:H])
424 lambdaFarZLow = max(maxEigenvaluesZZLow[0:1, 0:NV, 0:NV, 0:H], maxEigenvaluesZZLow[0:1, 0:NV, 0:NV, 0+1:H+1])
425
426 fluxLeftZLow = -1.0 * lambdaLeftZLow * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] - input[0:NumberOfUnknowns, 0-1:NV-1, 0:NV, 0:H])
427 fluxRightZLow = -1.0 * lambdaRightZLow * (input[0:NumberOfUnknowns, 0+1:NV+1, 0:NV, 0:H] - input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H])
428 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + dt / h(0) * (fluxLeftZLow - fluxRightZLow) * NV * 0.5
429
430 fluxBottomZLow = -1.0 * lambdaBottomZLow * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] - input[0:NumberOfUnknowns, 0:NV, 0-1:NV-1, 0:H])
431 fluxTopZLow = -1.0 * lambdaTopZLow * (input[0:NumberOfUnknowns, 0:NV, 0+1:NV+1, 0:H] - input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H])
432 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + dt / h(1) * (fluxBottomZLow - fluxTopZLow) * NV * 0.5
433
434 fluxNearZLow = -1.0 * lambdaNearZLow * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] - input[0:NumberOfUnknowns, 0:NV, 0:NV, 0-1:H-1])
435 fluxFarZLow = -1.0 * lambdaFarZLow * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 0+1:H+1] - input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H])
436 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + dt / h(2) * (fluxNearZLow - fluxFarZLow) * NV * 0.5
437
438 if flux:
439 flux_xZLow = DataBlock([[0, NumberOfUnknowns], [0-1, NV+1], [0, NV], [0, H], [0, N]], eval("flux", input[:, 0-1:NV+1, 0:NV, 0:H], x, 0))
440 flux_yZLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0-1, NV+1], [0, H], [0, N]], eval("flux", input[:, 0:NV, 0-1:NV+1, 0:H], x, 1))
441 flux_zZLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [0-1, H+1], [0, N]], eval("flux", input[:, 0:NV, 0:NV, 0-1:H+1], x, 2))
442 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + 0.5 * NV * dt / h(0) * (flux_xZLow[:, 0-1:NV-1, 0:NV, 0:H] - flux_xZLow[:, 0+1:NV+1, 0:NV, 0:H])
443 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + 0.5 * NV * dt / h(1) * (flux_yZLow[:, 0:NV, 0-1:NV-1, 0:H] - flux_yZLow[:, 0:NV, 0+1:NV+1, 0:H])
444 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + 0.5 * NV * dt / h(2) * (flux_zZLow[:, 0:NV, 0:NV, 0-1:H-1] - flux_zZLow[:, 0:NV, 0:NV, 0+1:H+1])
445
446 if nonconservativeProduct:
447 average_xZLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:NV, 0:H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:NV, 0:H])
448 delta_xZLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:NV, 0:H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:NV, 0:H]
449 nonconservative_xZLow = DataBlock([[0, NumberOfUnknowns], [0, NV+1], [0, NV], [0, H], [0, N]], eval("nonconservativeProduct", average_xZLow, delta_xZLow, x, h, t, dt, 0))
450 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + 0.5 * NV * dt / h(0) * (-nonconservative_xZLow[:, 0:NV, 0:NV, 0:H] - nonconservative_xZLow[:, 0+1:NV+1, 0:NV, 0:H])
451
452 average_yZLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:NV, 0:H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV+1, 0:H])
453 delta_yZLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV+1, 0:H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:NV, 0:H]
454 nonconservative_yZLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV+1], [0, H], [0, N]], eval("nonconservativeProduct", average_yZLow, delta_yZLow, x, h, t, dt, 1))
455 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + 0.5 * NV * dt / h(1) * (-nonconservative_yZLow[:, 0:NV, 0:NV, 0:H] - nonconservative_yZLow[:, 0:NV, 0+1:NV+1, 0:H])
456
457 average_zZLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, 0-1:H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, 0:H+1])
458 delta_zZLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, 0:H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, 0-1:H]
459 nonconservative_zZLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [0, H+1], [0, N]], eval("nonconservativeProduct", average_zZLow, delta_zZLow, x, h, t, dt, 2))
460 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:H] + 0.5 * NV * dt / h(2) * (-nonconservative_zZLow[:, 0:NV, 0:NV, 0:H] - nonconservative_zZLow[:, 0:NV, 0:NV, 0+1:H+1])
461
462 # Shell part ZHigh: [0,NV) x [0,NV) x [NV-H,NV)
463 output[:, 0:NV, 0:NV, NV-H:NV] = input[:, 0:NV, 0:NV, NV-H:NV]
464
465 if volumetricSource:
466 sourceZHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [NV-H, NV], [0, N]], eval("volumetricSource", input[:, 0:NV, 0:NV, NV-H:NV], x, t))
467 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + dt * sourceZHigh
468
469 maxEigenvaluesXZHigh = DataBlock([[0, 1], [0-1, NV+1], [0, NV], [NV-H, NV], [0, N]], eval("maxEigenvalue", input[:, 0-1:NV+1, 0:NV, NV-H:NV], x, 0))
470 maxEigenvaluesYZHigh = DataBlock([[0, 1], [0, NV], [0-1, NV+1], [NV-H, NV], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0-1:NV+1, NV-H:NV], x, 1))
471 maxEigenvaluesZZHigh = DataBlock([[0, 1], [0, NV], [0, NV], [NV-H-1, NV+1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0:NV, NV-H-1:NV+1], x, 2))
472
473 lambdaLeftZHigh = max(maxEigenvaluesXZHigh[0:1, 0-1:NV-1, 0:NV, NV-H:NV], maxEigenvaluesXZHigh[0:1, 0:NV, 0:NV, NV-H:NV])
474 lambdaRightZHigh = max(maxEigenvaluesXZHigh[0:1, 0:NV, 0:NV, NV-H:NV], maxEigenvaluesXZHigh[0:1, 0+1:NV+1, 0:NV, NV-H:NV])
475 lambdaBottomZHigh = max(maxEigenvaluesYZHigh[0:1, 0:NV, 0-1:NV-1, NV-H:NV], maxEigenvaluesYZHigh[0:1, 0:NV, 0:NV, NV-H:NV])
476 lambdaTopZHigh = max(maxEigenvaluesYZHigh[0:1, 0:NV, 0:NV, NV-H:NV], maxEigenvaluesYZHigh[0:1, 0:NV, 0+1:NV+1, NV-H:NV])
477 lambdaNearZHigh = max(maxEigenvaluesZZHigh[0:1, 0:NV, 0:NV, NV-H-1:NV-1], maxEigenvaluesZZHigh[0:1, 0:NV, 0:NV, NV-H:NV])
478 lambdaFarZHigh = max(maxEigenvaluesZZHigh[0:1, 0:NV, 0:NV, NV-H:NV], maxEigenvaluesZZHigh[0:1, 0:NV, 0:NV, NV-H+1:NV+1])
479
480 fluxLeftZHigh = -1.0 * lambdaLeftZHigh * (input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] - input[0:NumberOfUnknowns, 0-1:NV-1, 0:NV, NV-H:NV])
481 fluxRightZHigh = -1.0 * lambdaRightZHigh * (input[0:NumberOfUnknowns, 0+1:NV+1, 0:NV, NV-H:NV] - input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV])
482 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + dt / h(0) * (fluxLeftZHigh - fluxRightZHigh) * NV * 0.5
483
484 fluxBottomZHigh = -1.0 * lambdaBottomZHigh * (input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] - input[0:NumberOfUnknowns, 0:NV, 0-1:NV-1, NV-H:NV])
485 fluxTopZHigh = -1.0 * lambdaTopZHigh * (input[0:NumberOfUnknowns, 0:NV, 0+1:NV+1, NV-H:NV] - input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV])
486 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + dt / h(1) * (fluxBottomZHigh - fluxTopZHigh) * NV * 0.5
487
488 fluxNearZHigh = -1.0 * lambdaNearZHigh * (input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] - input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H-1:NV-1])
489 fluxFarZHigh = -1.0 * lambdaFarZHigh * (input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H+1:NV+1] - input[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV])
490 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + dt / h(2) * (fluxNearZHigh - fluxFarZHigh) * NV * 0.5
491
492 if flux:
493 flux_xZHigh = DataBlock([[0, NumberOfUnknowns], [0-1, NV+1], [0, NV], [NV-H, NV], [0, N]], eval("flux", input[:, 0-1:NV+1, 0:NV, NV-H:NV], x, 0))
494 flux_yZHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [0-1, NV+1], [NV-H, NV], [0, N]], eval("flux", input[:, 0:NV, 0-1:NV+1, NV-H:NV], x, 1))
495 flux_zZHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [NV-H-1, NV+1], [0, N]], eval("flux", input[:, 0:NV, 0:NV, NV-H-1:NV+1], x, 2))
496 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(0) * (flux_xZHigh[:, 0-1:NV-1, 0:NV, NV-H:NV] - flux_xZHigh[:, 0+1:NV+1, 0:NV, NV-H:NV])
497 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(1) * (flux_yZHigh[:, 0:NV, 0-1:NV-1, NV-H:NV] - flux_yZHigh[:, 0:NV, 0+1:NV+1, NV-H:NV])
498 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(2) * (flux_zZHigh[:, 0:NV, 0:NV, NV-H-1:NV-1] - flux_zZHigh[:, 0:NV, 0:NV, NV-H+1:NV+1])
499
500 if nonconservativeProduct:
501 average_xZHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:NV, NV-H:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:NV, NV-H:NV])
502 delta_xZHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:NV, NV-H:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:NV, NV-H:NV]
503 nonconservative_xZHigh = DataBlock([[0, NumberOfUnknowns], [0, NV+1], [0, NV], [NV-H, NV], [0, N]], eval("nonconservativeProduct", average_xZHigh, delta_xZHigh, x, h, t, dt, 0))
504 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(0) * (-nonconservative_xZHigh[:, 0:NV, 0:NV, NV-H:NV] - nonconservative_xZHigh[:, 0+1:NV+1, 0:NV, NV-H:NV])
505
506 average_yZHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:NV, NV-H:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV+1, NV-H:NV])
507 delta_yZHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV+1, NV-H:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:NV, NV-H:NV]
508 nonconservative_yZHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV+1], [NV-H, NV], [0, N]], eval("nonconservativeProduct", average_yZHigh, delta_yZHigh, x, h, t, dt, 1))
509 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(1) * (-nonconservative_yZHigh[:, 0:NV, 0:NV, NV-H:NV] - nonconservative_yZHigh[:, 0:NV, 0+1:NV+1, NV-H:NV])
510
511 average_zZHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, NV-H-1:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, NV-H:NV+1])
512 delta_zZHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, NV-H:NV+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, NV-H-1:NV]
513 nonconservative_zZHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [NV-H, NV+1], [0, N]], eval("nonconservativeProduct", average_zZHigh, delta_zZHigh, x, h, t, dt, 2))
514 output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, NV-H:NV] + 0.5 * NV * dt / h(2) * (-nonconservative_zZHigh[:, 0:NV, 0:NV, NV-H:NV] - nonconservative_zZHigh[:, 0:NV, 0:NV, NV-H+1:NV+1])
515
516 # Shell part YLow: [0,NV) x [0,H) x [H,NV-H)
517 output[:, 0:NV, 0:H, H:NV-H] = input[:, 0:NV, 0:H, H:NV-H]
518
519 if volumetricSource:
520 sourceYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, H], [H, NV-H], [0, N]], eval("volumetricSource", input[:, 0:NV, 0:H, H:NV-H], x, t))
521 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + dt * sourceYLow
522
523 maxEigenvaluesXYLow = DataBlock([[0, 1], [0-1, NV+1], [0, H], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, 0-1:NV+1, 0:H, H:NV-H], x, 0))
524 maxEigenvaluesYYLow = DataBlock([[0, 1], [0, NV], [0-1, H+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0-1:H+1, H:NV-H], x, 1))
525 maxEigenvaluesZYLow = DataBlock([[0, 1], [0, NV], [0, H], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0:H, H-1:NV-H+1], x, 2))
526
527 lambdaLeftYLow = max(maxEigenvaluesXYLow[0:1, 0-1:NV-1, 0:H, H:NV-H], maxEigenvaluesXYLow[0:1, 0:NV, 0:H, H:NV-H])
528 lambdaRightYLow = max(maxEigenvaluesXYLow[0:1, 0:NV, 0:H, H:NV-H], maxEigenvaluesXYLow[0:1, 0+1:NV+1, 0:H, H:NV-H])
529 lambdaBottomYLow = max(maxEigenvaluesYYLow[0:1, 0:NV, 0-1:H-1, H:NV-H], maxEigenvaluesYYLow[0:1, 0:NV, 0:H, H:NV-H])
530 lambdaTopYLow = max(maxEigenvaluesYYLow[0:1, 0:NV, 0:H, H:NV-H], maxEigenvaluesYYLow[0:1, 0:NV, 0+1:H+1, H:NV-H])
531 lambdaNearYLow = max(maxEigenvaluesZYLow[0:1, 0:NV, 0:H, H-1:NV-H-1], maxEigenvaluesZYLow[0:1, 0:NV, 0:H, H:NV-H])
532 lambdaFarYLow = max(maxEigenvaluesZYLow[0:1, 0:NV, 0:H, H:NV-H], maxEigenvaluesZYLow[0:1, 0:NV, 0:H, H+1:NV-H+1])
533
534 fluxLeftYLow = -1.0 * lambdaLeftYLow * (input[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] - input[0:NumberOfUnknowns, 0-1:NV-1, 0:H, H:NV-H])
535 fluxRightYLow = -1.0 * lambdaRightYLow * (input[0:NumberOfUnknowns, 0+1:NV+1, 0:H, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H])
536 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + dt / h(0) * (fluxLeftYLow - fluxRightYLow) * NV * 0.5
537
538 fluxBottomYLow = -1.0 * lambdaBottomYLow * (input[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, 0-1:H-1, H:NV-H])
539 fluxTopYLow = -1.0 * lambdaTopYLow * (input[0:NumberOfUnknowns, 0:NV, 0+1:H+1, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H])
540 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + dt / h(1) * (fluxBottomYLow - fluxTopYLow) * NV * 0.5
541
542 fluxNearYLow = -1.0 * lambdaNearYLow * (input[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, 0:H, H-1:NV-H-1])
543 fluxFarYLow = -1.0 * lambdaFarYLow * (input[0:NumberOfUnknowns, 0:NV, 0:H, H+1:NV-H+1] - input[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H])
544 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + dt / h(2) * (fluxNearYLow - fluxFarYLow) * NV * 0.5
545
546 if flux:
547 flux_xYLow = DataBlock([[0, NumberOfUnknowns], [0-1, NV+1], [0, H], [H, NV-H], [0, N]], eval("flux", input[:, 0-1:NV+1, 0:H, H:NV-H], x, 0))
548 flux_yYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0-1, H+1], [H, NV-H], [0, N]], eval("flux", input[:, 0:NV, 0-1:H+1, H:NV-H], x, 1))
549 flux_zYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, H], [H-1, NV-H+1], [0, N]], eval("flux", input[:, 0:NV, 0:H, H-1:NV-H+1], x, 2))
550 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_xYLow[:, 0-1:NV-1, 0:H, H:NV-H] - flux_xYLow[:, 0+1:NV+1, 0:H, H:NV-H])
551 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_yYLow[:, 0:NV, 0-1:H-1, H:NV-H] - flux_yYLow[:, 0:NV, 0+1:H+1, H:NV-H])
552 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + 0.5 * NV * dt / h(2) * (flux_zYLow[:, 0:NV, 0:H, H-1:NV-H-1] - flux_zYLow[:, 0:NV, 0:H, H+1:NV-H+1])
553
554 if nonconservativeProduct:
555 average_xYLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:H, H:NV-H])
556 delta_xYLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, 0:H, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, 0:H, H:NV-H]
557 nonconservative_xYLow = DataBlock([[0, NumberOfUnknowns], [0, NV+1], [0, H], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_xYLow, delta_xYLow, x, h, t, dt, 0))
558 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_xYLow[:, 0:NV, 0:H, H:NV-H] - nonconservative_xYLow[:, 0+1:NV+1, 0:H, H:NV-H])
559
560 average_yYLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H+1, H:NV-H])
561 delta_yYLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0-1:H, H:NV-H]
562 nonconservative_yYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, H+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_yYLow, delta_yYLow, x, h, t, dt, 1))
563 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_yYLow[:, 0:NV, 0:H, H:NV-H] - nonconservative_yYLow[:, 0:NV, 0+1:H+1, H:NV-H])
564
565 average_zYLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H, H:NV-H+1])
566 delta_zYLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:H, H-1:NV-H]
567 nonconservative_zYLow = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, H], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_zYLow, delta_zYLow, x, h, t, dt, 2))
568 output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, 0:H, H:NV-H] + 0.5 * NV * dt / h(2) * (-nonconservative_zYLow[:, 0:NV, 0:H, H:NV-H] - nonconservative_zYLow[:, 0:NV, 0:H, H+1:NV-H+1])
569
570 # Shell part YHigh: [0,NV) x [NV-H,NV) x [H,NV-H)
571 output[:, 0:NV, NV-H:NV, H:NV-H] = input[:, 0:NV, NV-H:NV, H:NV-H]
572
573 if volumetricSource:
574 sourceYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H, NV], [H, NV-H], [0, N]], eval("volumetricSource", input[:, 0:NV, NV-H:NV, H:NV-H], x, t))
575 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + dt * sourceYHigh
576
577 maxEigenvaluesXYHigh = DataBlock([[0, 1], [0-1, NV+1], [NV-H, NV], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, 0-1:NV+1, NV-H:NV, H:NV-H], x, 0))
578 maxEigenvaluesYYHigh = DataBlock([[0, 1], [0, NV], [NV-H-1, NV+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, 0:NV, NV-H-1:NV+1, H:NV-H], x, 1))
579 maxEigenvaluesZYHigh = DataBlock([[0, 1], [0, NV], [NV-H, NV], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, NV-H:NV, H-1:NV-H+1], x, 2))
580
581 lambdaLeftYHigh = max(maxEigenvaluesXYHigh[0:1, 0-1:NV-1, NV-H:NV, H:NV-H], maxEigenvaluesXYHigh[0:1, 0:NV, NV-H:NV, H:NV-H])
582 lambdaRightYHigh = max(maxEigenvaluesXYHigh[0:1, 0:NV, NV-H:NV, H:NV-H], maxEigenvaluesXYHigh[0:1, 0+1:NV+1, NV-H:NV, H:NV-H])
583 lambdaBottomYHigh = max(maxEigenvaluesYYHigh[0:1, 0:NV, NV-H-1:NV-1, H:NV-H], maxEigenvaluesYYHigh[0:1, 0:NV, NV-H:NV, H:NV-H])
584 lambdaTopYHigh = max(maxEigenvaluesYYHigh[0:1, 0:NV, NV-H:NV, H:NV-H], maxEigenvaluesYYHigh[0:1, 0:NV, NV-H+1:NV+1, H:NV-H])
585 lambdaNearYHigh = max(maxEigenvaluesZYHigh[0:1, 0:NV, NV-H:NV, H-1:NV-H-1], maxEigenvaluesZYHigh[0:1, 0:NV, NV-H:NV, H:NV-H])
586 lambdaFarYHigh = max(maxEigenvaluesZYHigh[0:1, 0:NV, NV-H:NV, H:NV-H], maxEigenvaluesZYHigh[0:1, 0:NV, NV-H:NV, H+1:NV-H+1])
587
588 fluxLeftYHigh = -1.0 * lambdaLeftYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] - input[0:NumberOfUnknowns, 0-1:NV-1, NV-H:NV, H:NV-H])
589 fluxRightYHigh = -1.0 * lambdaRightYHigh * (input[0:NumberOfUnknowns, 0+1:NV+1, NV-H:NV, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H])
590 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + dt / h(0) * (fluxLeftYHigh - fluxRightYHigh) * NV * 0.5
591
592 fluxBottomYHigh = -1.0 * lambdaBottomYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, NV-H-1:NV-1, H:NV-H])
593 fluxTopYHigh = -1.0 * lambdaTopYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H+1:NV+1, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H])
594 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + dt / h(1) * (fluxBottomYHigh - fluxTopYHigh) * NV * 0.5
595
596 fluxNearYHigh = -1.0 * lambdaNearYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] - input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H-1:NV-H-1])
597 fluxFarYHigh = -1.0 * lambdaFarYHigh * (input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H+1:NV-H+1] - input[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H])
598 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + dt / h(2) * (fluxNearYHigh - fluxFarYHigh) * NV * 0.5
599
600 if flux:
601 flux_xYHigh = DataBlock([[0, NumberOfUnknowns], [0-1, NV+1], [NV-H, NV], [H, NV-H], [0, N]], eval("flux", input[:, 0-1:NV+1, NV-H:NV, H:NV-H], x, 0))
602 flux_yYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H-1, NV+1], [H, NV-H], [0, N]], eval("flux", input[:, 0:NV, NV-H-1:NV+1, H:NV-H], x, 1))
603 flux_zYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H, NV], [H-1, NV-H+1], [0, N]], eval("flux", input[:, 0:NV, NV-H:NV, H-1:NV-H+1], x, 2))
604 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_xYHigh[:, 0-1:NV-1, NV-H:NV, H:NV-H] - flux_xYHigh[:, 0+1:NV+1, NV-H:NV, H:NV-H])
605 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_yYHigh[:, 0:NV, NV-H-1:NV-1, H:NV-H] - flux_yYHigh[:, 0:NV, NV-H+1:NV+1, H:NV-H])
606 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(2) * (flux_zYHigh[:, 0:NV, NV-H:NV, H-1:NV-H-1] - flux_zYHigh[:, 0:NV, NV-H:NV, H+1:NV-H+1])
607
608 if nonconservativeProduct:
609 average_xYHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, NV-H:NV, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, NV-H:NV, H:NV-H])
610 delta_xYHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+1, NV-H:NV, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:NV, NV-H:NV, H:NV-H]
611 nonconservative_xYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV+1], [NV-H, NV], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_xYHigh, delta_xYHigh, x, h, t, dt, 0))
612 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_xYHigh[:, 0:NV, NV-H:NV, H:NV-H] - nonconservative_xYHigh[:, 0+1:NV+1, NV-H:NV, H:NV-H])
613
614 average_yYHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H-1:NV, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV+1, H:NV-H])
615 delta_yYHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H-1:NV, H:NV-H]
616 nonconservative_yYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H, NV+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_yYHigh, delta_yYHigh, x, h, t, dt, 1))
617 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_yYHigh[:, 0:NV, NV-H:NV, H:NV-H] - nonconservative_yYHigh[:, 0:NV, NV-H+1:NV+1, H:NV-H])
618
619 average_zYHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV, H:NV-H+1])
620 delta_zYHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, NV-H:NV, H-1:NV-H]
621 nonconservative_zYHigh = DataBlock([[0, NumberOfUnknowns], [0, NV], [NV-H, NV], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_zYHigh, delta_zYHigh, x, h, t, dt, 2))
622 output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] = output[0:NumberOfUnknowns, 0:NV, NV-H:NV, H:NV-H] + 0.5 * NV * dt / h(2) * (-nonconservative_zYHigh[:, 0:NV, NV-H:NV, H:NV-H] - nonconservative_zYHigh[:, 0:NV, NV-H:NV, H+1:NV-H+1])
623
624 # Shell part XLow: [0,H) x [H,NV-H) x [H,NV-H)
625 output[:, 0:H, H:NV-H, H:NV-H] = input[:, 0:H, H:NV-H, H:NV-H]
626
627 if volumetricSource:
628 sourceXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H, NV-H], [H, NV-H], [0, N]], eval("volumetricSource", input[:, 0:H, H:NV-H, H:NV-H], x, t))
629 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + dt * sourceXLow
630
631 maxEigenvaluesXXLow = DataBlock([[0, 1], [0-1, H+1], [H, NV-H], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, 0-1:H+1, H:NV-H, H:NV-H], x, 0))
632 maxEigenvaluesYXLow = DataBlock([[0, 1], [0, H], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, 0:H, H-1:NV-H+1, H:NV-H], x, 1))
633 maxEigenvaluesZXLow = DataBlock([[0, 1], [0, H], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, 0:H, H:NV-H, H-1:NV-H+1], x, 2))
634
635 lambdaLeftXLow = max(maxEigenvaluesXXLow[0:1, 0-1:H-1, H:NV-H, H:NV-H], maxEigenvaluesXXLow[0:1, 0:H, H:NV-H, H:NV-H])
636 lambdaRightXLow = max(maxEigenvaluesXXLow[0:1, 0:H, H:NV-H, H:NV-H], maxEigenvaluesXXLow[0:1, 0+1:H+1, H:NV-H, H:NV-H])
637 lambdaBottomXLow = max(maxEigenvaluesYXLow[0:1, 0:H, H-1:NV-H-1, H:NV-H], maxEigenvaluesYXLow[0:1, 0:H, H:NV-H, H:NV-H])
638 lambdaTopXLow = max(maxEigenvaluesYXLow[0:1, 0:H, H:NV-H, H:NV-H], maxEigenvaluesYXLow[0:1, 0:H, H+1:NV-H+1, H:NV-H])
639 lambdaNearXLow = max(maxEigenvaluesZXLow[0:1, 0:H, H:NV-H, H-1:NV-H-1], maxEigenvaluesZXLow[0:1, 0:H, H:NV-H, H:NV-H])
640 lambdaFarXLow = max(maxEigenvaluesZXLow[0:1, 0:H, H:NV-H, H:NV-H], maxEigenvaluesZXLow[0:1, 0:H, H:NV-H, H+1:NV-H+1])
641
642 fluxLeftXLow = -1.0 * lambdaLeftXLow * (input[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, 0-1:H-1, H:NV-H, H:NV-H])
643 fluxRightXLow = -1.0 * lambdaRightXLow * (input[0:NumberOfUnknowns, 0+1:H+1, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H])
644 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + dt / h(0) * (fluxLeftXLow - fluxRightXLow) * NV * 0.5
645
646 fluxBottomXLow = -1.0 * lambdaBottomXLow * (input[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, 0:H, H-1:NV-H-1, H:NV-H])
647 fluxTopXLow = -1.0 * lambdaTopXLow * (input[0:NumberOfUnknowns, 0:H, H+1:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H])
648 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + dt / h(1) * (fluxBottomXLow - fluxTopXLow) * NV * 0.5
649
650 fluxNearXLow = -1.0 * lambdaNearXLow * (input[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, 0:H, H:NV-H, H-1:NV-H-1])
651 fluxFarXLow = -1.0 * lambdaFarXLow * (input[0:NumberOfUnknowns, 0:H, H:NV-H, H+1:NV-H+1] - input[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H])
652 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + dt / h(2) * (fluxNearXLow - fluxFarXLow) * NV * 0.5
653
654 if flux:
655 flux_xXLow = DataBlock([[0, NumberOfUnknowns], [0-1, H+1], [H, NV-H], [H, NV-H], [0, N]], eval("flux", input[:, 0-1:H+1, H:NV-H, H:NV-H], x, 0))
656 flux_yXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("flux", input[:, 0:H, H-1:NV-H+1, H:NV-H], x, 1))
657 flux_zXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("flux", input[:, 0:H, H:NV-H, H-1:NV-H+1], x, 2))
658 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_xXLow[:, 0-1:H-1, H:NV-H, H:NV-H] - flux_xXLow[:, 0+1:H+1, H:NV-H, H:NV-H])
659 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_yXLow[:, 0:H, H-1:NV-H-1, H:NV-H] - flux_yXLow[:, 0:H, H+1:NV-H+1, H:NV-H])
660 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(2) * (flux_zXLow[:, 0:H, H:NV-H, H-1:NV-H-1] - flux_zXLow[:, 0:H, H:NV-H, H+1:NV-H+1])
661
662 if nonconservativeProduct:
663 average_xXLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:H, H:NV-H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H+1, H:NV-H, H:NV-H])
664 delta_xXLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H+1, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0-1:H, H:NV-H, H:NV-H]
665 nonconservative_xXLow = DataBlock([[0, NumberOfUnknowns], [0, H+1], [H, NV-H], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_xXLow, delta_xXLow, x, h, t, dt, 0))
666 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_xXLow[:, 0:H, H:NV-H, H:NV-H] - nonconservative_xXLow[:, 0+1:H+1, H:NV-H, H:NV-H])
667
668 average_yXLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H-1:NV-H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H+1, H:NV-H])
669 delta_yXLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H-1:NV-H, H:NV-H]
670 nonconservative_yXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H, NV-H+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_yXLow, delta_yXLow, x, h, t, dt, 1))
671 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_yXLow[:, 0:H, H:NV-H, H:NV-H] - nonconservative_yXLow[:, 0:H, H+1:NV-H+1, H:NV-H])
672
673 average_zXLow = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H, H:NV-H+1])
674 delta_zXLow = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:H, H:NV-H, H-1:NV-H]
675 nonconservative_zXLow = DataBlock([[0, NumberOfUnknowns], [0, H], [H, NV-H], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_zXLow, delta_zXLow, x, h, t, dt, 2))
676 output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, 0:H, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(2) * (-nonconservative_zXLow[:, 0:H, H:NV-H, H:NV-H] - nonconservative_zXLow[:, 0:H, H:NV-H, H+1:NV-H+1])
677
678 # Shell part XHigh: [NV-H,NV) x [H,NV-H) x [H,NV-H)
679 output[:, NV-H:NV, H:NV-H, H:NV-H] = input[:, NV-H:NV, H:NV-H, H:NV-H]
680
681 if volumetricSource:
682 sourceXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H, NV-H], [H, NV-H], [0, N]], eval("volumetricSource", input[:, NV-H:NV, H:NV-H, H:NV-H], x, t))
683 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + dt * sourceXHigh
684
685 maxEigenvaluesXXHigh = DataBlock([[0, 1], [NV-H-1, NV+1], [H, NV-H], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, NV-H-1:NV+1, H:NV-H, H:NV-H], x, 0))
686 maxEigenvaluesYXHigh = DataBlock([[0, 1], [NV-H, NV], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("maxEigenvalue", input[:, NV-H:NV, H-1:NV-H+1, H:NV-H], x, 1))
687 maxEigenvaluesZXHigh = DataBlock([[0, 1], [NV-H, NV], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("maxEigenvalue", input[:, NV-H:NV, H:NV-H, H-1:NV-H+1], x, 2))
688
689 lambdaLeftXHigh = max(maxEigenvaluesXXHigh[0:1, NV-H-1:NV-1, H:NV-H, H:NV-H], maxEigenvaluesXXHigh[0:1, NV-H:NV, H:NV-H, H:NV-H])
690 lambdaRightXHigh = max(maxEigenvaluesXXHigh[0:1, NV-H:NV, H:NV-H, H:NV-H], maxEigenvaluesXXHigh[0:1, NV-H+1:NV+1, H:NV-H, H:NV-H])
691 lambdaBottomXHigh = max(maxEigenvaluesYXHigh[0:1, NV-H:NV, H-1:NV-H-1, H:NV-H], maxEigenvaluesYXHigh[0:1, NV-H:NV, H:NV-H, H:NV-H])
692 lambdaTopXHigh = max(maxEigenvaluesYXHigh[0:1, NV-H:NV, H:NV-H, H:NV-H], maxEigenvaluesYXHigh[0:1, NV-H:NV, H+1:NV-H+1, H:NV-H])
693 lambdaNearXHigh = max(maxEigenvaluesZXHigh[0:1, NV-H:NV, H:NV-H, H-1:NV-H-1], maxEigenvaluesZXHigh[0:1, NV-H:NV, H:NV-H, H:NV-H])
694 lambdaFarXHigh = max(maxEigenvaluesZXHigh[0:1, NV-H:NV, H:NV-H, H:NV-H], maxEigenvaluesZXHigh[0:1, NV-H:NV, H:NV-H, H+1:NV-H+1])
695
696 fluxLeftXHigh = -1.0 * lambdaLeftXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, NV-H-1:NV-1, H:NV-H, H:NV-H])
697 fluxRightXHigh = -1.0 * lambdaRightXHigh * (input[0:NumberOfUnknowns, NV-H+1:NV+1, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H])
698 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + dt / h(0) * (fluxLeftXHigh - fluxRightXHigh) * NV * 0.5
699
700 fluxBottomXHigh = -1.0 * lambdaBottomXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, NV-H:NV, H-1:NV-H-1, H:NV-H])
701 fluxTopXHigh = -1.0 * lambdaTopXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H+1:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H])
702 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + dt / h(1) * (fluxBottomXHigh - fluxTopXHigh) * NV * 0.5
703
704 fluxNearXHigh = -1.0 * lambdaNearXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H-1:NV-H-1])
705 fluxFarXHigh = -1.0 * lambdaFarXHigh * (input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H+1:NV-H+1] - input[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H])
706 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + dt / h(2) * (fluxNearXHigh - fluxFarXHigh) * NV * 0.5
707
708 if flux:
709 flux_xXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H-1, NV+1], [H, NV-H], [H, NV-H], [0, N]], eval("flux", input[:, NV-H-1:NV+1, H:NV-H, H:NV-H], x, 0))
710 flux_yXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H-1, NV-H+1], [H, NV-H], [0, N]], eval("flux", input[:, NV-H:NV, H-1:NV-H+1, H:NV-H], x, 1))
711 flux_zXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H, NV-H], [H-1, NV-H+1], [0, N]], eval("flux", input[:, NV-H:NV, H:NV-H, H-1:NV-H+1], x, 2))
712 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (flux_xXHigh[:, NV-H-1:NV-1, H:NV-H, H:NV-H] - flux_xXHigh[:, NV-H+1:NV+1, H:NV-H, H:NV-H])
713 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (flux_yXHigh[:, NV-H:NV, H-1:NV-H-1, H:NV-H] - flux_yXHigh[:, NV-H:NV, H+1:NV-H+1, H:NV-H])
714 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(2) * (flux_zXHigh[:, NV-H:NV, H:NV-H, H-1:NV-H-1] - flux_zXHigh[:, NV-H:NV, H:NV-H, H+1:NV-H+1])
715
716 if nonconservativeProduct:
717 average_xXHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H-1:NV, H:NV-H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV+1, H:NV-H, H:NV-H])
718 delta_xXHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV+1, H:NV-H, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H-1:NV, H:NV-H, H:NV-H]
719 nonconservative_xXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV+1], [H, NV-H], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_xXHigh, delta_xXHigh, x, h, t, dt, 0))
720 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(0) * (-nonconservative_xXHigh[:, NV-H:NV, H:NV-H, H:NV-H] - nonconservative_xXHigh[:, NV-H+1:NV+1, H:NV-H, H:NV-H])
721
722 average_yXHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H-1:NV-H, H:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H+1, H:NV-H])
723 delta_yXHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H+1, H:NV-H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H-1:NV-H, H:NV-H]
724 nonconservative_yXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H, NV-H+1], [H, NV-H], [0, N]], eval("nonconservativeProduct", average_yXHigh, delta_yXHigh, x, h, t, dt, 1))
725 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(1) * (-nonconservative_yXHigh[:, NV-H:NV, H:NV-H, H:NV-H] - nonconservative_yXHigh[:, NV-H:NV, H+1:NV-H+1, H:NV-H])
726
727 average_zXHigh = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H, H-1:NV-H] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H, H:NV-H+1])
728 delta_zXHigh = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H, H:NV-H+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, NV-H:NV, H:NV-H, H-1:NV-H]
729 nonconservative_zXHigh = DataBlock([[0, NumberOfUnknowns], [NV-H, NV], [H, NV-H], [H, NV-H+1], [0, N]], eval("nonconservativeProduct", average_zXHigh, delta_zXHigh, x, h, t, dt, 2))
730 output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] = output[0:NumberOfUnknowns, NV-H:NV, H:NV-H, H:NV-H] + 0.5 * NV * dt / h(2) * (-nonconservative_zXHigh[:, NV-H:NV, H:NV-H, H:NV-H] - nonconservative_zXHigh[:, NV-H:NV, H:NV-H, H+1:NV-H+1])
731
732 if maxEigenvalue:
733 # Reduce per shell part, as max() over one DataBlock overwrites its target
734 maxEigenvaluesX2ZLow = DataBlock([[0, 1], [0, NV], [0, NV], [0, H], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:NV, 0:H], x, 0))
735 maxEigenvaluesY2ZLow = DataBlock([[0, 1], [0, NV], [0, NV], [0, H], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:NV, 0:H], x, 1))
736 maxEigenvaluesZ2ZLow = DataBlock([[0, 1], [0, NV], [0, NV], [0, H], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:NV, 0:H], x, 2))
737 maxEigenvaluesXYZLow = max(maxEigenvaluesX2ZLow, maxEigenvaluesY2ZLow)
738 maxEigenvaluesPerVolumeZLow = max(maxEigenvaluesXYZLow, maxEigenvaluesZ2ZLow)
739 maxEigenvalueZLow = DataBlock([[0, N]], 0.0)
740 maxEigenvalueZLow = max(maxEigenvaluesPerVolumeZLow)
741 maxEigenvaluesX2ZHigh = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV-(NV-H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:NV, NV-H:NV], x, 0))
742 maxEigenvaluesY2ZHigh = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV-(NV-H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:NV, NV-H:NV], x, 1))
743 maxEigenvaluesZ2ZHigh = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV-(NV-H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:NV, NV-H:NV], x, 2))
744 maxEigenvaluesXYZHigh = max(maxEigenvaluesX2ZHigh, maxEigenvaluesY2ZHigh)
745 maxEigenvaluesPerVolumeZHigh = max(maxEigenvaluesXYZHigh, maxEigenvaluesZ2ZHigh)
746 maxEigenvalueZHigh = DataBlock([[0, N]], 0.0)
747 maxEigenvalueZHigh = max(maxEigenvaluesPerVolumeZHigh)
748 maxEigenvaluesX2YLow = DataBlock([[0, 1], [0, NV], [0, H], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:H, H:NV-H], x, 0))
749 maxEigenvaluesY2YLow = DataBlock([[0, 1], [0, NV], [0, H], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:H, H:NV-H], x, 1))
750 maxEigenvaluesZ2YLow = DataBlock([[0, 1], [0, NV], [0, H], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, 0:H, H:NV-H], x, 2))
751 maxEigenvaluesXYYLow = max(maxEigenvaluesX2YLow, maxEigenvaluesY2YLow)
752 maxEigenvaluesPerVolumeYLow = max(maxEigenvaluesXYYLow, maxEigenvaluesZ2YLow)
753 maxEigenvalueYLow = DataBlock([[0, N]], 0.0)
754 maxEigenvalueYLow = max(maxEigenvaluesPerVolumeYLow)
755 maxEigenvaluesX2YHigh = DataBlock([[0, 1], [0, NV], [0, NV-(NV-H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, NV-H:NV, H:NV-H], x, 0))
756 maxEigenvaluesY2YHigh = DataBlock([[0, 1], [0, NV], [0, NV-(NV-H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, NV-H:NV, H:NV-H], x, 1))
757 maxEigenvaluesZ2YHigh = DataBlock([[0, 1], [0, NV], [0, NV-(NV-H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:NV, NV-H:NV, H:NV-H], x, 2))
758 maxEigenvaluesXYYHigh = max(maxEigenvaluesX2YHigh, maxEigenvaluesY2YHigh)
759 maxEigenvaluesPerVolumeYHigh = max(maxEigenvaluesXYYHigh, maxEigenvaluesZ2YHigh)
760 maxEigenvalueYHigh = DataBlock([[0, N]], 0.0)
761 maxEigenvalueYHigh = max(maxEigenvaluesPerVolumeYHigh)
762 maxEigenvaluesX2XLow = DataBlock([[0, 1], [0, H], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:H, H:NV-H, H:NV-H], x, 0))
763 maxEigenvaluesY2XLow = DataBlock([[0, 1], [0, H], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:H, H:NV-H, H:NV-H], x, 1))
764 maxEigenvaluesZ2XLow = DataBlock([[0, 1], [0, H], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, 0:H, H:NV-H, H:NV-H], x, 2))
765 maxEigenvaluesXYXLow = max(maxEigenvaluesX2XLow, maxEigenvaluesY2XLow)
766 maxEigenvaluesPerVolumeXLow = max(maxEigenvaluesXYXLow, maxEigenvaluesZ2XLow)
767 maxEigenvalueXLow = DataBlock([[0, N]], 0.0)
768 maxEigenvalueXLow = max(maxEigenvaluesPerVolumeXLow)
769 maxEigenvaluesX2XHigh = DataBlock([[0, 1], [0, NV-(NV-H)], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, NV-H:NV, H:NV-H, H:NV-H], x, 0))
770 maxEigenvaluesY2XHigh = DataBlock([[0, 1], [0, NV-(NV-H)], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, NV-H:NV, H:NV-H, H:NV-H], x, 1))
771 maxEigenvaluesZ2XHigh = DataBlock([[0, 1], [0, NV-(NV-H)], [0, NV-H-(H)], [0, NV-H-(H)], [0, N]], eval("maxEigenvalue", output[:, NV-H:NV, H:NV-H, H:NV-H], x, 2))
772 maxEigenvaluesXYXHigh = max(maxEigenvaluesX2XHigh, maxEigenvaluesY2XHigh)
773 maxEigenvaluesPerVolumeXHigh = max(maxEigenvaluesXYXHigh, maxEigenvaluesZ2XHigh)
774 maxEigenvalueXHigh = DataBlock([[0, N]], 0.0)
775 maxEigenvalueXHigh = max(maxEigenvaluesPerVolumeXHigh)
776 maxEigenvalues = max(maxEigenvalueZLow, maxEigenvalueZHigh)
777 maxEigenvalues = max(maxEigenvalues, maxEigenvalueYLow)
778 maxEigenvalues = max(maxEigenvalues, maxEigenvalueYHigh)
779 maxEigenvalues = max(maxEigenvalues, maxEigenvalueXLow)
780 maxEigenvalues = max(maxEigenvalues, maxEigenvalueXHigh)
781
782class RusanovWithOptimisticExplicitEuler(Solver.NumericalScheme):
783 """!
784
785 Finite Volumes with a Rusanov Riemann solver and optimistic explicit Euler time
786 stepping, realised over a block-structured mesh.
787
788 ## Design idea
789
790 The finite volume solver operates over a block-structured mesh, where each
791 cell of the mesh holds a patch of finite volumes. This solver's fundamental
792 idea is that it splits up the patch into three classes of volumes. Assume
793 we have a NxNxN patch of volumes, where N is the patch_size, and this patch
794 is surrounded by a halo of size H. So we have in total a (N+2H)x(N+2H)x(N+2H)
795 patch that feeds into the time stepping as reconstructed patch.
796 We now distinguish three classes of volumes:
797
798 - The halo cells
799 - The shell cells, which is the outer layer of the NxNxN arrangement.
800 - The patch interior cells, which are the remaining volumes in the interior
801 of the patch, so effectively a (N-2H)x(N-2H)x(N-2H) patch.
802
803 The sketch below shows the three classes for a 2d patch with N=10 and H=1.
804 It is generated by
805 src/exahype2/api/fv/rusanov-with-optimistic-explicit-euler-illustrations.py.
806
807 @image html rusanov-with-optimistic-explicit-euler-volume-classes.png
808
809 For the new solver, we now define two new update kernels:
810
811 - timeStepWithRusanov2DOnPatchInterior which is timeStepWithRusanov2D but acts only on
812 the patch interior.
813 - timeStepWithRusanov2DOnShell which works only on the shell.
814
815 The 3D variants follow equivalently.
816
817 ## Time stepping
818
819 Every time step either runs the traditional branch (Update cell and
820 Discard patch interior) or the optimistic branch (Update cell shell and
821 Insert patch interior). The traditional branch runs if the solver's time
822 step size is zero, or if the estimated time step size, i.e. the second
823 entry of the eigenvalue data, exceeds it. The optimistic branch runs
824 otherwise. Project on faces, Compute global CFL condition and Spawn
825 patch interior update follow either branch. If at least one axis of the
826 domain is not periodic and the solver uses Dirichlet boundary
827 conditions, Apply Dirichlet BCs sets the ghost volumes along the domain
828 boundary before Reconstruct reads them (@see
829 exahype2.api.fv.actionsets.ApplyDirichletBoundaryConditions). With
830 Sommerfeld boundary conditions, the task does not exist, and the roll
831 over within Project on faces advances the ghost volumes instead (@see
832 exahype2.api.fv.actionsets.RollOverUpdatedFaceWithSommerfeldBoundaryConditions).
833 If all axes are periodic, there is no boundary treatment at all. The
834 task graph below is
835 generated by
836 src/exahype2/api/fv/rusanov-with-optimistic-explicit-euler-illustrations.py
837 from this class for a domain without periodic boundaries.
838
839 @image html rusanov-with-optimistic-explicit-euler-time-step-task-graph.png
840
841 The estimate is the admissible time step size of the current time step,
842 i.e. the one the solver could have chosen if it had not been enforced to
843 use a smaller one, times damp_time_step_size. Reconstruct
844 memorises the admissible time step size per cell (@see
845 exahype2.api.fv.actionsets.ReadjustTimeStepSizeForOptimisticGlobalTimestepping),
846 and Spawn patch interior update damps it (@see
847 exahype2.api.fv.actionsets.SpawnTaskToUpdateCellInterior). As the
848 admissible time step size is never damped itself, the estimates do not
849 shrink from time step to time step:
850
851 - Time step 0 runs with time step size 0, as the admissible time step
852 size is not known yet. Update cell runs but does not change the
853 solution, and Compute global CFL condition yields the first time step
854 size dt0. There is no admissible time step size to estimate from, so
855 the spawn does not update the patch interior.
856 - Time step 1 runs with dt0 down the traditional branch. The spawn
857 updates the patch interior with damp_time_step_size*dt0.
858 - Time step 2 finds this estimate to be smaller than the admissible time
859 step size dt1. It enforces the estimate as time step size and runs
860 the optimistic branch.
861
862 Compute global CFL condition does not only compute the next time step
863 size, it also makes it the solver's time step size and time stamp
864 straightaway (@see exahype2.api.fv.actionsets.ComputeGlobalCFLCondition).
865 It therefore depends on the spawn, which still sees the time stamp and
866 time step size of the current time step and starts the patch interior
867 update from their sum. The order in which the tasks are added to the
868 task graph does not matter.
869
870 If spawn_optimistic_calculations_as_task is set, the spawn runs the
871 patch interior update as a task. Insert patch interior waits for this
872 task before it copies the patch interior into the patch, and Discard
873 patch interior waits for it, too, so no task is still running when the
874 next one is spawned. The task numbers are held per cell.
875
876 """
877
878 def __init__(self, solver: Solver, patch_size: int, overlap: int = 1, cfl_constant: float = 0.9, damp_time_step_size: float = 0.95, spawn_optimistic_calculations_as_task: bool = True, min_volume_size: float = None, max_volume_size: float = None):
879 """!
880
881 The constructor memorises the attributes. However, it does not do this
882 within this class, but it adds the important properties straight to the
883 solver, so anybody else working against this solver knows about the
884 properties introduced by the numerical scheme.
885
886 @see @ref __create_data_structures()
887
888 patch_size: int
889 Number of Finite Volumes per axis within one patch (i.e. within
890 one mesh cell).
891
892 overlap: int
893 Number of Finite Volumes the halo layer reconstructed around a
894 patch - from the neighbouring patches' faces - is wide.
895
896 cfl_constant: float
897 Safety factor the admissible time step size (@see
898 exahype2.api.fv.actionsets.ComputeGlobalCFLCondition,
899 which this is passed straight through to) is scaled down by,
900 relative to the theoretical CFL=1 stability limit
901 h / (Dimensions * lambda_max). 0.9 is a
902 conservative default for an explicit-Euler Finite Volume scheme;
903 this codebase has no prior precedent value to match (grep for
904 "cfl" outside this class finds only exaparty's unrelated SPH
905 solvers, which default to 0.1).
906
907 damp_time_step_size: float
908 Factor, smaller or equal to one, which the current time step size
909 is multiplied with to estimate the time step size of the patch
910 interior update (@see
911 exahype2.api.fv.actionsets.SpawnTaskToUpdateCellInterior).
912
913 spawn_optimistic_calculations_as_task: bool
914 If True, Spawn patch interior update spawns the patch interior
915 update as a task, and Insert patch interior and Discard patch
916 interior wait for this task. The task number is held per cell
917 (@see __create_data_structures()). If False, Spawn patch interior
918 update runs the patch interior update straightaway.
919
920 min_volume_size: float
921 The time step asks for a refinement of a cell only if its volumes
922 are at least of this size (@ref exahype2.api.blockstructured.DynamicMeshRefinement).
923 None (the default) does not bound the refinement.
924
925 max_volume_size: float
926 The time step asks for an erase, i.e. coarsening, of a cell only
927 if the volumes of the coarser cell would be strictly smaller than
928 this size. None (the default) does not bound the coarsening.
929
930 """
931 if damp_time_step_size > 1.0:
932 raise Exception(f"damp_time_step_size has to be smaller or equal to 1, but is {damp_time_step_size}")
933 super().__init__(solver)
934 assert (
935 min_volume_size is None
936 or max_volume_size is None
937 or min_volume_size <= max_volume_size
938 ), f"min_volume_size={min_volume_size} has to be smaller than or equal to max_volume_size={max_volume_size}"
939 self.min_volume_size = min_volume_size
940 self.max_volume_size = max_volume_size
941 self.patch_size = patch_size
942 self.overlap = overlap
943 # Set in construct_output(): the generated kernels live in the project
944 # namespace, as their names are solver-specific
945 self.kernel_namespace = None
946 self.cfl_constant = cfl_constant
947 self.damp_time_step_size = damp_time_step_size
948 self.spawn_optimistic_calculations_as_task = spawn_optimistic_calculations_as_task
949 self.__create_data_structures()
950
951 @property
952 def kernel_function_name(self) -> str:
953 """!
954 Name of the generated time step kernel. It carries the solver name,
955 as do the kernels/<SolverName>_rusanov*.h files holding it, so two
956 solvers' kernels do not clash.
957 """
958 return self._solver_solver.name + "_timeStepWithRusanov"
959
960 def __create_data_structures(self):
961 """!
962
963 Create the mesh data structures (a Finite Volume patch per cell,
964 plus the face-based overlap patches used to reconstruct its halo)
965 this numerical scheme requires.
966
967 Ported from the legacy exahype2.solvers.fv.FV.create_data_structures():
968 one @ref peano4.api.datamodel.Patch of patch_size^d Finite Volumes
969 per mesh cell (_patch), plus a halo-augmented copy of it
970 (_patch_with_halo, (patch_size+2*overlap)^d, @see
971 toolbox.blockstructured.api.ReconstructPatchWithHalo) and two overlap
972 patches of width 2*overlap attached to the mesh faces:
973
974 - _patch_overlap_current_time_step: the persistent halo data
975 toolbox.blockstructured.api.ReconstructPatchWithHalo reads from to
976 rebuild _patch_with_halo (@see construct_output(), which hands it
977 _patch/_patch_overlap_current_time_step/_patch_with_halo directly
978 rather than the whole solver).
979 - _patch_overlap_update: the projection of the just-updated _patch
980 accumulated by the current sweep (@see
981 exahype2.api.fv.actionsets.ProjectPatchOntoFaces), rolled over
982 into _patch_overlap_current_time_step by
983 exahype2.api.fv.actionsets.RollOverUpdatedFace once the update is
984 complete.
985 - _max_eigenvalue_and_time_step_sizes: four doubles per cell, named
986 after the solver followed by MaxEigenvalueAndTimeStepSizes (e.g.
987 EulerFVSmoothMaxEigenvalueAndTimeStepSizes).
988 exahype2.api.fv.RusanovWithExplicitEuler calls the corresponding
989 data _max_eigenvalue. exahype2.api.fv.actionsets.UpdateCell and
990 exahype2.api.fv.actionsets.ComputeGlobalCFLCondition are shared
991 with that scheme and hence are handed this data explicitly.
992 It is modelled via DaStGen2 as an array of cardinality four. The
993 array is a peano4.api.dastgen2.Peano4SmartPointerDoubleArray,
994 i.e. the values reside on the heap. The four entries are used as
995 follows:
996 - Entry 0 holds the maximum eigenvalue of the patch.
997 exahype2.api.fv.actionsets.UpdateCell writes the one of the
998 whole patch into it, while
999 exahype2.api.fv.actionsets.UpdateCellShell writes the one of
1000 the shell only. exahype2.api.fv.actionsets.InsertPatchInterior
1001 eventually replaces it with the maximum of entries 0 and 2, so
1002 it then holds the one of the whole patch, too.
1003 exahype2.api.fv.actionsets.ComputeGlobalCFLCondition reduces
1004 entry 0 of all cells into the next time step size.
1005 - Entry 1 holds the estimated time step size of the patch
1006 interior update, i.e. entry 3 times damp_time_step_size, as
1007 written by
1008 exahype2.api.fv.actionsets.SpawnTaskToUpdateCellInterior.
1009 - Entry 2 holds the maximum eigenvalue of the patch interior, as
1010 written by exahype2.api.fv.actionsets.UpdateCellInterior when
1011 exahype2.api.fv.actionsets.SpawnTaskToUpdateCellInterior
1012 updates the patch interior.
1013 - Entry 3 holds the time step size the solver could have chosen
1014 in the current time step if it had not been enforced to use a
1015 smaller one, i.e. the admissible time step size, as written by
1016 exahype2.api.fv.actionsets.ReadjustTimeStepSizeForOptimisticGlobalTimestepping.
1017 - _patch_interior: the patch interior, i.e. the patch without its
1018 shell, so (patch_size-2*overlap)^d Finite Volumes. It matches the
1019 region timeStepWithRusanov2DOnPatchInterior() and
1020 timeStepWithRusanov3DOnPatchInterior() update.
1021 - _optimistic_cell_update_task_number: only if
1022 spawn_optimistic_calculations_as_task is set, otherwise None. One
1023 integer per cell, named after the solver followed by
1024 OptimisticCellUpdateTaskNumber (e.g.
1025 EulerFVSmoothOptimisticCellUpdateTaskNumber), which holds the
1026 number of the task updating the patch interior, or
1027 tarch::multicore::TaskEnumerator::NoNumber if there is none. The
1028 class furthermore holds a static tarch::multicore::TaskEnumerator
1029 (getTaskEnumerator()) handing out the task numbers.
1030
1031 All six are assigned onto self._solver (the shared
1032 exahype2.api.Solver instance both this numerical scheme and its
1033 surrounding solver's Solver.GridConstruction/Solver.Plotter
1034 counterparts hold a reference to as self._solver), not onto self,
1035 since the remaining FV action sets that need them - @see
1036 exahype2.api.fv.actionsets.ApplyDirichletBoundaryConditions,
1037 exahype2.api.fv.actionsets.UpdateCell,
1038 exahype2.api.fv.actionsets.ComputeGlobalCFLCondition - are handed that
1039 shared solver and look the patches up as solver._patch /
1040 solver._patch_overlap_*.
1041
1042 Storage is hardwired to smart pointers (peano4.api.datamodel's
1043 PatchToDoubleArrayWithSmartPointer), matching the legacy default
1044 (Storage.SmartPointers for both cell and face data) - the legacy
1045 Storage enum's other two backends (CallStack, Heap) are not
1046 exposed here.
1047
1048 Two pieces of the legacy method are deliberately not ported yet:
1049 - The cell/face label creation
1050 (exahype2.grid.create_cell_label()/create_face_label()) - that
1051 module has not been ported to the new API.
1052 - The load/store/compute guards
1053 (FV._provide_cell_data_to_compute_kernels_default_guard() and
1054 friends) - they are keyed off the legacy SolverState state
1055 machine, which the task-graph architecture has no equivalent of:
1056 a task only ever runs within its own phase's task graph in the
1057 first place, so "are we still constructing the grid" does not
1058 need to be re-checked at runtime the same way.
1059
1060 """
1061 total_unknowns = self._solver_solver.unknowns + self._solver_solver.auxiliary_variables
1062 identifier = self._solver_solver._unknown_identifier()
1063
1064 self._solver_solver._patch = peano4.api.datamodel.Patch(
1065 (self.patch_size, self.patch_size, self.patch_size),
1066 total_unknowns,
1067 identifier,
1068 )
1069 self._solver_solver._patch_with_halo = peano4.api.datamodel.Patch(
1070 (self.patch_size+2*self.overlap, self.patch_size+2*self.overlap, self.patch_size+2*self.overlap),
1071 total_unknowns,
1072 identifier + "WithHalo",
1073 )
1074 self._solver_solver._patch_interior = peano4.api.datamodel.Patch(
1075 (self.patch_size-2*self.overlap, self.patch_size-2*self.overlap, self.patch_size-2*self.overlap),
1076 total_unknowns,
1077 identifier + "PatchInterior",
1078 )
1079 self._solver_solver._patch_overlap_current_time_step = peano4.api.datamodel.Patch(
1080 (2 * self.overlap, self.patch_size, self.patch_size),
1081 total_unknowns,
1082 identifier + "CurrentTimeStep",
1083 )
1084 self._solver_solver._patch_overlap_update = peano4.api.datamodel.Patch(
1085 (2 * self.overlap, self.patch_size, self.patch_size),
1086 total_unknowns,
1087 identifier + "Update",
1088 )
1089 self._solver_solver._max_eigenvalue_and_time_step_sizes = peano4.api.datamodel.DaStGen2(self._solver_solver.name + "MaxEigenvalueAndTimeStepSizes")
1090 self._solver_solver._max_eigenvalue_and_time_step_sizes.data.add_attribute(
1091 peano4.api.dastgen2.Peano4SmartPointerDoubleArray("maxEigenvalue", 4)
1092 )
1093
1094 self._solver_solver._optimistic_cell_update_task_number = None
1095 if self.spawn_optimistic_calculations_as_task:
1096 self._solver_solver._optimistic_cell_update_task_number = peano4.api.datamodel.DaStGen2(
1097 self._solver_solver.name + "OptimisticCellUpdateTaskNumber"
1098 )
1099 self._solver_solver._optimistic_cell_update_task_number.data.add_attribute(
1100 dastgen2.attributes.Integer(name="taskNumber")
1101 )
1102 self._solver_solver._optimistic_cell_update_task_number.data.add_attribute(
1103 dastgen2.attributes.UserDefinedType(
1104 name="taskEnumerator",
1105 type="tarch::multicore::TaskEnumerator",
1106 include='#include "tarch/multicore/TaskEnumerator.h"',
1107 qualifier=dastgen2.attributes.Attribute.Qualifier.STATIC,
1108 user_type_has_toString_method=False,
1109 )
1110 )
1111
1112
1113 def __add_data_to_repository(self, peano_project):
1114 """!
1115
1116 Consult @ref page_exahype_programmers_guide_new_solver for some information
1117 how this routine fits to the typical workflow when we introduce a new
1118 solver.
1119
1120 """
1121 peano4.api.Logger().log("""
1122Patch data
1123----------
1124{}
1125Patch overlap data
1126------------------
1127{}
1128{}""".format( str(self._solver_solver._patch), str(self._solver_solver._patch_overlap_current_time_step), str(self._solver_solver._patch_overlap_update) ) )
1129
1130 peano_project.datamodel.add_cell(self._solver_solver._patch)
1131 peano_project.datamodel.add_cell(self._solver_solver._patch_with_halo)
1132 peano_project.datamodel.add_cell(self._solver_solver._patch_interior)
1133 peano_project.datamodel.add_face(self._solver_solver._patch_overlap_current_time_step)
1134 peano_project.datamodel.add_face(self._solver_solver._patch_overlap_update)
1135 peano_project.datamodel.add_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes)
1136 if self._solver_solver._optimistic_cell_update_task_number is not None:
1137 peano_project.datamodel.add_cell(self._solver_solver._optimistic_cell_update_task_number)
1138 self._solver_solver._optimistic_cell_update_task_number.grid_selector = self._solver_solver.grid_selector
1139
1140 # The task graph compiler only uses, loads and stores data where the
1141 # solver's grid selector holds, too (@see peano4.api.datamodel.DoF).
1142 for patch in (
1143 self._solver_solver._patch,
1144 self._solver_solver._patch_with_halo,
1145 self._solver_solver._patch_interior,
1146 self._solver_solver._patch_overlap_current_time_step,
1147 self._solver_solver._patch_overlap_update,
1148 self._solver_solver._max_eigenvalue_and_time_step_sizes,
1149 ):
1150 patch.grid_selector = self._solver_solver.grid_selector
1151
1152 self._solver_solver._patch.generator =peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
1153 self._solver_solver._patch, "double"
1154 )
1155 self._solver_solver._patch_with_halo.generator = peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
1156 self._solver_solver._patch_with_halo, "double"
1157 )
1158 self._solver_solver._patch_interior.generator = peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
1159 self._solver_solver._patch_interior, "double"
1160 )
1161 self._solver_solver._patch_overlap_current_time_step.generator = peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
1162 self._solver_solver._patch_overlap_current_time_step, "double"
1163 )
1164 self._solver_solver._patch_overlap_update.generator = peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
1165 self._solver_solver._patch_overlap_update, "double"
1166 )
1167
1168 self._solver_solver._patch_overlap_current_time_step.generator.merge_method_definition = (
1169 toolbox.blockstructured.api.get_face_merge_implementation(
1170 self._solver_solver._patch_overlap_current_time_step
1171 )
1172 )
1173 self._solver_solver._patch_overlap_current_time_step.generator.includes += """
1174#include "peano4/utils/Loop.h"
1175"""
1176
1177 # Only faces selected by the solver's grid selector hold valid data,
1178 # so only those are exchanged.
1179 self._solver_solver._patch_overlap_current_time_step.generator.send_condition = self._solver_solver.grid_selector.createFaceGuard()
1180 self._solver_solver._patch_overlap_current_time_step.generator.receive_and_merge_condition = self._solver_solver.grid_selector.createFaceGuard()
1181
1182 self._solver_solver._patch_overlap_update.generator.send_condition = "false"
1183 self._solver_solver._patch_overlap_update.generator.receive_and_merge_condition = "false"
1184
1185 def create_default_mesh_adaptation(self) -> Solver.MeshAdaptation:
1186 """!
1187
1188 The Finite Volume schemes adapt the mesh and balance the load
1189 through exahype2.api.fv.AdaptMeshAndBalanceLoad by default
1190
1191 """
1192 from .AdaptMeshAndBalanceLoad import AdaptMeshAndBalanceLoad
1193
1194 return AdaptMeshAndBalanceLoad(self._solver_solver)
1195
1196 def construct_output(self,
1197 exahype2_project,
1198 peano_project,
1199 time_step):
1200 """!
1201
1202 """
1203 peano4.api.Logger().log("Create time stepping task graph for Finite Volume scheme over block-structured meshes with Rusanov flux")
1204
1205 self.__add_data_to_repository(peano_project)
1206
1207 self.kernel_namespace = peano_project.namespace
1208 # Construct the call first, as it validates the PDE terms' variants
1209 self._solver_solver.compute_kernel_call = self.__construct_kernel_call(Target.Serial, False)
1210 self._solver_solver.compute_kernel_call_on_patch_interior = self.__construct_kernel_call(Target.Serial, False, self.kernel_function_name + "OnPatchInterior")
1211 self._solver_solver.compute_kernel_call_on_shell = self.__construct_kernel_call(Target.Serial, False, self.kernel_function_name + "OnShell")
1212 self.__create_time_step_functions(peano_project.directory)
1213
1214 reconstruction_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Reconstruct", entity=peano4.api.taskgraph.Task.GridEntity.CELL)
1215 update_cell_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Update cell", entity=peano4.api.taskgraph.Task.GridEntity.CELL)
1216 update_cell_shell_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Update cell shell", entity=peano4.api.taskgraph.Task.GridEntity.CELL)
1217 insert_patch_interior_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Insert patch interior", entity=peano4.api.taskgraph.Task.GridEntity.CELL, may_spawn_background_tasks=False)
1218 discard_patch_interior_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Discard patch interior", entity=peano4.api.taskgraph.Task.GridEntity.CELL, may_spawn_background_tasks=False)
1219 project_on_faces_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Project on faces", entity=peano4.api.taskgraph.Task.GridEntity.CELL)
1220 compute_cfl_condition_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Compute global CFL condition", entity=peano4.api.taskgraph.Task.GridEntity.CELL, may_spawn_background_tasks=False)
1221 spawn_patch_interior_task = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Spawn patch interior update", entity=peano4.api.taskgraph.Task.GridEntity.CELL, may_spawn_background_tasks=False)
1222
1223 time_step.add_task(reconstruction_task)
1224 time_step.add_task(update_cell_task)
1225 time_step.add_task(update_cell_shell_task)
1226 time_step.add_task(insert_patch_interior_task)
1227 time_step.add_task(discard_patch_interior_task)
1228 time_step.add_task(project_on_faces_task)
1229 time_step.add_task(compute_cfl_condition_task)
1230 time_step.add_task(spawn_patch_interior_task)
1231
1232 # The load balancing is switched off throughout the time step, and
1233 # switched on again only at its very end, so it may then decide on
1234 # splits. Peano carries these splits out in the subsequent mesh sweep.
1235 reconstruction_task.add_action_set(
1236 toolbox.loadbalancing.api.actionsets.ToggleLoadBalancing(
1237 enable_load_balancing=False, plug_into_prepare_traversal=True
1238 )
1239 )
1240
1241 # The very last task evaluates the adaptivity criterion if the
1242 # solver has one (create_dynamic_mesh_refinement_action_set()), i.e. it
1243 # collects refine and erase requests for the mesh adaptation phase,
1244 # and it switches the load balancing on again and informs it that
1245 # the time step has finished.
1246 dynamic_mesh_refinement_task = peano4.api.taskgraph.Task(
1247 name=f"{self._solver.name}_Dynamic mesh refinement",
1248 entity=peano4.api.taskgraph.Task.GridEntity.CELL,
1249 may_spawn_background_tasks=False,
1250 )
1251 dynamic_mesh_refinement_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1252 dynamic_mesh_refinement = create_dynamic_mesh_refinement_action_set(
1253 self._solver_solver,
1254 self.patch_size,
1255 guard=self._solver_solver.grid_selector.createCellGuard(),
1256 min_volume_size=self.min_volume_size,
1257 max_volume_size=self.max_volume_size,
1258 )
1259 if dynamic_mesh_refinement is not None:
1260 dynamic_mesh_refinement_task.add_action_set(dynamic_mesh_refinement)
1261 dynamic_mesh_refinement_task.add_action_set(
1262 toolbox.loadbalancing.api.actionsets.LoadBalancing(enable_load_balancing=True)
1263 )
1264 dynamic_mesh_refinement_task.add_dependency(compute_cfl_condition_task)
1265 dynamic_mesh_refinement_task.add_dependency(project_on_faces_task)
1266 # No other task depends on the discard, so we have to make it a
1267 # predecessor explicitly, as the dynamic mesh refinement has to be
1268 # the very last task
1269 dynamic_mesh_refinement_task.add_dependency(discard_patch_interior_task)
1270 time_step.add_task(dynamic_mesh_refinement_task)
1271
1272
1273 reconstruction_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1274 reconstruction_task.use_cell(self._solver_solver._patch_with_halo, peano4.api.taskgraph.Task.DataAccess.WRITE)
1275
1277 reconstruction_task.use_face(
1278 self._solver_solver._patch_overlap_current_time_step,
1279 get_access_to_persistent_face_data(self._solver_solver, exahype2_project, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD),
1280 )
1281 reconstruction_task.add_action_set(
1282 toolbox.blockstructured.api.ReconstructPatchWithHalo(
1283 patch=self._solver_solver._patch,
1284 patch_overlap=self._solver_solver._patch_overlap_current_time_step,
1285 reconstructed_patch=self._solver_solver._patch_with_halo,
1286 guard=self._solver_solver.grid_selector.createCellGuard(),
1287 add_assertions_to_halo_exchange=False,
1288 additional_includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
1289 )
1290 )
1291
1294 ApplyDirichletBoundaryConditions.add_task_to_time_step(self._solver_solver, exahype2_project, time_step, reconstruction_task)
1295
1297 reconstruction_task.use_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1298 reconstruction_task.add_action_set(
1299 ReadjustTimeStepSizeForOptimisticGlobalTimestepping(
1300 self._solver_solver,
1301 guard=self._solver_solver.grid_selector.createCellGuard(),
1302 )
1303 )
1304
1305
1321 max_eigenvalue = f"fineGridCell{self._solver._max_eigenvalue_and_time_step_sizes.name}"
1322 # Cell data holding the number of the task which updates the patch
1323 # interior, or None if the patch interior is updated straightaway
1324 task_number = self._solver_solver._optimistic_cell_update_task_number
1325 time_step_size = f"repositories::{self._solver.get_name_of_global_instance()}.getMinTimeStepSize()"
1326 cell_guard = self._solver_solver.grid_selector.createCellGuard()
1327
1332 traditional_condition = (
1333 f"(tarch::la::equals({time_step_size}, 0.0) or tarch::la::greater({max_eigenvalue}.getMaxEigenvalue(1), {time_step_size}))"
1334 )
1335 traditional_guard = f"({cell_guard}) and {traditional_condition}"
1336 optimistic_guard = f"({cell_guard}) and not {traditional_condition}"
1337
1338
1339
1340 update_cell_task.add_dependency(reconstruction_task)
1341 update_cell_task.use_cell(self._solver_solver._patch_with_halo, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD)
1342 update_cell_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.WRITE)
1343 update_cell_task.use_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1344 update_cell_task.add_action_set(
1345 UpdateCell(
1346 self._solver_solver,
1347 guard=traditional_guard,
1348 max_eigenvalue=self._solver_solver._max_eigenvalue_and_time_step_sizes,
1349 )
1350 )
1351 update_cell_task.tags.add(Solver.SolutionUpdateTaskTag)
1352
1353
1354 discard_patch_interior_task.add_dependency(update_cell_task)
1355 if task_number is not None:
1356 discard_patch_interior_task.use_cell(task_number, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD)
1357 discard_patch_interior_task.use_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1358 discard_patch_interior_task.add_action_set(
1359 DiscardPatchInterior(
1360 self._solver_solver,
1361 guard=traditional_guard,
1362 task_number=task_number,
1363 )
1364 )
1365
1366
1369 update_cell_shell_task.add_dependency(reconstruction_task)
1370 update_cell_shell_task.use_cell(self._solver_solver._patch_with_halo, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD)
1371 update_cell_shell_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1372 update_cell_shell_task.use_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1373 update_cell_shell_task.add_action_set(
1374 UpdateCellShell(
1375 self._solver_solver,
1376 guard=optimistic_guard,
1377 )
1378 )
1379
1380
1381 insert_patch_interior_task.add_dependency(update_cell_shell_task)
1382 if task_number is not None:
1383 insert_patch_interior_task.use_cell(task_number, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD)
1384 insert_patch_interior_task.use_cell(self._solver_solver._patch_interior, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD)
1385 insert_patch_interior_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1386 insert_patch_interior_task.use_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1387 insert_patch_interior_task.add_action_set(
1388 InsertPatchInterior(
1389 self._solver_solver,
1390 self.patch_size,
1391 self.overlap,
1392 guard=optimistic_guard,
1393 task_number=task_number,
1394 )
1395 )
1396 insert_patch_interior_task.tags.add(Solver.SolutionUpdateTaskTag)
1397
1398
1400 compute_cfl_condition_task.add_dependency(update_cell_task)
1401 compute_cfl_condition_task.add_dependency(insert_patch_interior_task)
1402 # The spawn has to see the solver's time stamp and time step size
1403 # before the CFL condition advances them
1404 compute_cfl_condition_task.add_dependency(spawn_patch_interior_task)
1405 compute_cfl_condition_task.use_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1406 compute_cfl_condition_task.add_action_set(
1407 ComputeGlobalCFLCondition(
1408 self._solver_solver,
1409 self.cfl_constant,
1410 self.patch_size,
1411 guard=cell_guard,
1412 max_eigenvalue=self._solver_solver._max_eigenvalue_and_time_step_sizes,
1413 )
1414 )
1415
1416
1417 project_on_faces_task.add_dependency(update_cell_task)
1418 project_on_faces_task.add_dependency(insert_patch_interior_task)
1419 project_on_faces_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1420 project_on_faces_task.use_face(self._solver_solver._patch_overlap_update, peano4.api.taskgraph.Task.DataAccess.TEMPORARY)
1421 project_on_faces_task.use_face(
1422 self._solver_solver._patch_overlap_current_time_step,
1423 get_access_to_persistent_face_data(self._solver_solver, exahype2_project, peano4.api.taskgraph.Task.DataAccess.WRITE),
1424 )
1425 project_on_faces_task.add_action_set(
1426 ProjectPatchOntoFaces(
1427 self._solver_solver._patch,
1428 self._solver_solver._patch_overlap_update,
1429 guard=cell_guard,
1430 includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
1431 )
1432 )
1433
1435 project_on_faces_task.add_action_set(
1436 RollOverUpdatedFaceWithSommerfeldBoundaryConditions.create_roll_over(
1437 self._solver_solver,
1438 exahype2_project,
1439 guard=self._solver_solver.grid_selector.createFaceGuard(),
1440 includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
1441 time_step_size=f"repositories::{self._solver.get_name_of_global_instance()}.getMinTimeStepSize()",
1442 )
1443 )
1444
1445
1448 spawn_patch_interior_task.add_dependency(update_cell_task)
1449 spawn_patch_interior_task.add_dependency(insert_patch_interior_task)
1450 spawn_patch_interior_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1451 spawn_patch_interior_task.use_cell(self._solver_solver._max_eigenvalue_and_time_step_sizes, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
1452 spawn_patch_interior_task.use_cell(self._solver_solver._patch_interior, peano4.api.taskgraph.Task.DataAccess.WRITE)
1453 if task_number is not None:
1454 spawn_patch_interior_task.use_cell(task_number, peano4.api.taskgraph.Task.DataAccess.WRITE)
1455 spawn_patch_interior_task.add_action_set(
1456 SpawnTaskToUpdateCellInterior(
1457 self._solver_solver,
1458 self.damp_time_step_size,
1459 guard=cell_guard,
1460 task_number=task_number,
1461 )
1462 )
1463
1464 if peano4.api.Logger().is_visualisation:
1465 time_step.visualise(title="Time stepping task graph for {}".format(self._solver_solver.name))
1466 if peano4.api.Logger().is_text:
1467 time_step.log_as_text(title="Time stepping task graph for {}".format(self._solver_solver.name))
1468
1469 def __create_time_step_functions(self, output_directory):
1470 """!
1471
1472 Calling this function is delayed until all PDE term information is provided.
1473 In this way we can avoid generating code that will not be used. The following
1474 information is required to translate a Python function using DSLHyPE:
1475 - template_parameters: the variables that are used within the Python function
1476 to represent external information
1477 - functor_parameters: the PDE terms that can be used by the solver
1478 - function: the Python function itself
1479
1480 TODO: Currently the MLIR types are passed manually here, this should ideally
1481 be inferred
1482
1483 """
1484
1485 template_parameters = [
1486 dslhype.SyntaxTree.Argument("NV", dslhype.SyntaxTree.TInteger()),
1487 dslhype.SyntaxTree.Argument("H", dslhype.SyntaxTree.TInteger()),
1488 dslhype.SyntaxTree.Argument("NumberOfUnknowns", dslhype.SyntaxTree.TInteger()),
1489 dslhype.SyntaxTree.Argument("NumberOfAuxiliaryVariables", dslhype.SyntaxTree.TInteger())
1490 ]
1491
1492 functor_parameters = [
1493 dslhype.SyntaxTree.Argument("flux", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::FluxFunctor&"), "(!llvm.ptr, !llvm.ptr, !llvm.ptr, f64, f64, i32, !llvm.ptr) -> ()"),
1494 dslhype.SyntaxTree.Argument("nonconservativeProduct", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::NonconservativeProductFunctor&"), "(!llvm.ptr, !llvm.ptr, !llvm.ptr, !llvm.ptr, f64, f64, i32, !llvm.ptr) -> ()"),
1495 dslhype.SyntaxTree.Argument("volumetricSource", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::VolumetricSourceFunctor&"), "(!llvm.ptr, !llvm.ptr, f64, !llvm.ptr) -> ()"),
1496 dslhype.SyntaxTree.Argument("maxEigenvalue", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::MaxEigenvalueFunctor&"), "(!llvm.ptr, !llvm.ptr, !llvm.ptr, f64, f64, i32, !llvm.ptr) -> ()")
1497 ]
1498
1499 time_step_dsl_function_2d = exahype2.api.fv.timeStepWithRusanov2D
1500 time_step_dsl_function_3d = exahype2.api.fv.timeStepWithRusanov3D
1501 dslhype.translate_dslhype_function(
1502 function = time_step_dsl_function_2d,
1503 function_namespace = self.kernel_namespace,
1504 template_parameters=template_parameters,
1505 functor_parameters=functor_parameters,
1506 pde_terms=self._solver_solver.pde_terms,
1507 output_file_name=self._solver_solver.name + "_rusanov2d",
1508 target=dslhype.Target.Serial,
1509 stateless=False,
1510 function_name=self.kernel_function_name,
1511 output_directory=output_directory,
1512 )
1513 dslhype.translate_dslhype_function(
1514 function = time_step_dsl_function_3d,
1515 function_namespace = self.kernel_namespace,
1516 template_parameters=template_parameters,
1517 functor_parameters=functor_parameters,
1518 pde_terms=self._solver_solver.pde_terms,
1519 output_file_name=self._solver_solver.name + "_rusanov3d",
1520 target=dslhype.Target.Serial,
1521 stateless=False,
1522 function_name=self.kernel_function_name,
1523 output_directory=output_directory,
1524 )
1525 # Kernels on the patch interior and on the shell
1526 for function, file_suffix, function_suffix in [
1527 (timeStepWithRusanov2DOnPatchInterior, "_rusanov2d_patch_interior", "OnPatchInterior"),
1528 (timeStepWithRusanov3DOnPatchInterior, "_rusanov3d_patch_interior", "OnPatchInterior"),
1529 (timeStepWithRusanov2DOnShell, "_rusanov2d_shell", "OnShell"),
1530 (timeStepWithRusanov3DOnShell, "_rusanov3d_shell", "OnShell"),
1531 ]:
1532 dslhype.translate_dslhype_function(
1533 function = function,
1534 function_namespace = self.kernel_namespace,
1535 template_parameters=template_parameters,
1536 functor_parameters=functor_parameters,
1537 pde_terms=self._solver_solver.pde_terms,
1538 output_file_name=self._solver_solver.name + file_suffix,
1539 target=dslhype.Target.Serial,
1540 stateless=False,
1541 function_name=self.kernel_function_name + function_suffix,
1542 output_directory=output_directory,
1543 )
1544 file = open(os.path.join(output_directory, "kernels", f"{self._solver.name}_rusanov.h"), "w")
1545 file.write(f"""#if Dimensions == 2
1546 #include "{self._solver.name}_rusanov2d.h"
1547 #include "{self._solver.name}_rusanov2d_patch_interior.h"
1548 #include "{self._solver.name}_rusanov2d_shell.h"
1549#else
1550 #include "{self._solver.name}_rusanov3d.h"
1551 #include "{self._solver.name}_rusanov3d_patch_interior.h"
1552 #include "{self._solver.name}_rusanov3d_shell.h"
1553#endif""")
1554 file.close()
1555
1556 def __construct_kernel_call(self, target: Target, stateless: bool, function_name: str = None):
1557 """!
1558
1559 This function constructs the calls to the compute kernels based on the
1560 application PDE terms and the version chosen.
1561
1562 @param function_name: str
1563 Name of the kernel to call. Defaults to kernel_function_name, i.e.
1564 the kernel updating the whole patch.
1565
1566 """
1567 available_terms = ["flux", "nonconservativeProduct", "volumetricSource", "maxEigenvalue"]
1568 if function_name is None:
1569 function_name = self.kernel_function_name
1570
1571 if stateless:
1572 return f"""::{"::".join(self.kernel_namespace)}::{function_name}<
1573 {self.patch_size},
1574 {self.overlap},
1575 {self._solver.unknowns},
1576 {self._solver.auxiliary_variables}
1577 >(
1578 patchData,
1579 {self._solver.name}
1580 );"""
1581 else:
1582 # Variants per term the kernel takes: {name: {stateless: PDETerm}}
1583 variants = {}
1584 for application_term in self._solver_solver.pde_terms:
1585 if application_term.NAME in available_terms:
1586 variants.setdefault(application_term.NAME, {})[application_term.stateless] = application_term
1587
1588 terms_with_both_variants = [name for name in variants if len(variants[name]) == 2]
1589 terms_without_stateless_variant = [name for name in variants if True not in variants[name]]
1590 if terms_with_both_variants and terms_without_stateless_variant:
1591 raise Exception(
1592 f"Solver {self._solver.name}: the PDE terms {terms_with_both_variants} are "
1593 f"registered both with and without state. In this case, the generated code "
1594 f"asks cellCanUseStatelessPDETerms() per cell before it invokes the kernel: "
1595 f"if it returns true, the kernel is called with the stateless variants of all "
1596 f"PDE terms, otherwise with the variants with state. However, the PDE terms "
1597 f"{terms_without_stateless_variant} have no stateless variant, so a kernel "
1598 f"call would mix stateless and stateful routines. Either add stateless "
1599 f"variants of {terms_without_stateless_variant}, or register "
1600 f"{terms_with_both_variants} only once."
1601 )
1602
1603 def kernel_call(stateless: bool) -> str:
1604 pde_term_calls = []
1605 for term in available_terms:
1606 if term not in variants:
1607 continue
1608 # A term registered only once is used in either case
1609 application_term = variants[term].get(stateless, next(iter(variants[term].values())))
1610 # Commas inside template brackets, e.g. Vector<Dimensions, double>,
1611 # do not separate parameters
1612 arguments = [
1613 parameter.split()[-1]
1614 for parameter in re.split(r",(?![^<]*>)", application_term.SIGNATURE)
1615 ]
1616 # Stateless terms are static and take a dummy tag (@see PDETerm)
1617 if application_term.stateless:
1618 arguments.append("::exahype2::Solver::Offloadable::Yes")
1619 pde_term_calls.append(f"""
1620 [&](
1621 {application_term.SIGNATURE}
1622 )->void {{
1623 repositories::{self._solver.get_name_of_global_instance()}.{application_term.name}({", ".join(arguments)});
1624 }}""")
1625 return f"""::{"::".join(self.kernel_namespace)}::{function_name}<
1626 {self.patch_size},
1627 {self.overlap},
1628 {self._solver.unknowns},
1629 {self._solver.auxiliary_variables}
1630>(
1631 patchData,
1632{",\n".join(pde_term_calls)}
1633);"""
1634
1635 if not terms_with_both_variants:
1636 return kernel_call(stateless=False)
1637 return f"""if (repositories::{self._solver.get_name_of_global_instance()}.cellCanUseStatelessPDETerms(
1638 patchData.cellCentre[0],
1639 patchData.cellSize[0],
1640 patchData.t[0],
1641 patchData.dt[0]
1642)) {{
1643{kernel_call(stateless=True)}
1644}} else {{
1645{kernel_call(stateless=False)}
1646}}"""
UpdateCell.UpdateCell
Action set implementing the cell update of the continuous Galerkin solver.
Definition UpdateCell.py:9
dastgen2.attributes.UserDefinedType.UserDefinedType
Wrapper around user-defined attribute.
Definition UserDefinedType.py:6
exahype2.api.Solver.Solver.NumericalScheme
The concrete discretisation (Finite Volumes, ADER-DG, ...) realised by a Solver.
Definition Solver.py:60
exahype2.api.fv.AdaptMeshAndBalanceLoad.AdaptMeshAndBalanceLoad
Default mesh adaptation of the Finite Volume schemes.
Definition AdaptMeshAndBalanceLoad.py:13
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler
Finite Volumes with a Rusanov Riemann solver and optimistic explicit Euler time stepping,...
Definition RusanovWithOptimisticExplicitEuler.py:782
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.create_default_mesh_adaptation
Solver.MeshAdaptation create_default_mesh_adaptation(self)
The Finite Volume schemes adapt the mesh and balance the load through exahype2.api....
Definition RusanovWithOptimisticExplicitEuler.py:1185
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.patch_size
patch_size
Optimistic branch: insert the patch interior once the shell is updated.
Definition RusanovWithOptimisticExplicitEuler.py:941
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.construct_output
construct_output(self, exahype2_project, peano_project, time_step)
Has to be implemented by subclass.
Definition RusanovWithOptimisticExplicitEuler.py:1199
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler._solver
_solver
Build the reconstruction tasks befilling the patch_with_halo, which is logically temporary.
Definition RusanovWithOptimisticExplicitEuler.py:1194
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.__init__
__init__(self, Solver solver, int patch_size, int overlap=1, float cfl_constant=0.9, float damp_time_step_size=0.95, bool spawn_optimistic_calculations_as_task=True, float min_volume_size=None, float max_volume_size=None)
The constructor memorises the attributes.
Definition RusanovWithOptimisticExplicitEuler.py:878
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.cfl_constant
cfl_constant
Both branches: reduce the cells' maximum eigenvalues into the next time step size.
Definition RusanovWithOptimisticExplicitEuler.py:946
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.__add_data_to_repository
__add_data_to_repository(self, peano_project)
Consult Introduce a new numerical scheme for some information how this routine fits to the typical wo...
Definition RusanovWithOptimisticExplicitEuler.py:1113
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.__construct_kernel_call
__construct_kernel_call(self, Target target, bool stateless, str function_name=None)
This function constructs the calls to the compute kernels based on the application PDE terms and the ...
Definition RusanovWithOptimisticExplicitEuler.py:1556
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.__create_time_step_functions
__create_time_step_functions(self, output_directory)
Calling this function is delayed until all PDE term information is provided.
Definition RusanovWithOptimisticExplicitEuler.py:1469
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.RusanovWithOptimisticExplicitEuler.__create_data_structures
__create_data_structures(self)
Create the mesh data structures (a Finite Volume patch per cell, plus the face-based overlap patches ...
Definition RusanovWithOptimisticExplicitEuler.py:960
exahype2.api.fv.actionsets.ComputeGlobalCFLCondition.ComputeGlobalCFLCondition
Reduce the cells' maximum eigenvalues and derive the next global time step size from the CFL conditio...
Definition ComputeGlobalCFLCondition.py:8
exahype2.api.fv.actionsets.ProjectPatchOntoFaces.ProjectPatchOntoFaces
Project a Finite Volume patch's data onto the overlap held on its surrounding faces,...
Definition ProjectPatchOntoFaces.py:6
exahype2.api.fv.actionsets.UpdateCellShell.UpdateCellShell
Update the shell of a Finite Volume patch.
Definition UpdateCellShell.py:8
peano4.api.Logger.Logger
Singleton logger for the Peano 4 Python API.
Definition Logger.py:37
peano4.api.dastgen2.Peano4SmartPointerDoubleArray.Peano4SmartPointerDoubleArray
Specialisation of array using Peano's tarch with a smart pointer.
Definition Peano4SmartPointerDoubleArray.py:6
peano4.api.datamodel.DaStGen2.DaStGen2
Default superclass for any data model in Peano which is stored within the grid.
Definition DaStGen2.py:188
peano4.api.datamodel.Patch.Patch
Patch implementation of a degree of freedom.
Definition Patch.py:7
peano4.api.taskgraph.Task.Task
One task within a peano4.api.taskgraph.Graph.
Definition Task.py:8
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.timeStepWithRusanov3DOnShell
"void" timeStepWithRusanov3DOnShell("::exahype2::CellData<double, double>&" patchData)
Rusanov update of the shell of a 3d patch.
Definition RusanovWithOptimisticExplicitEuler.py:397
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.timeStepWithRusanov2DOnPatchInterior
"void" timeStepWithRusanov2DOnPatchInterior("::exahype2::CellData<double, double>&" patchData)
DSLHyPE time step kernels on the patch interior and on the shell of a patch.
Definition RusanovWithOptimisticExplicitEuler.py:69
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.timeStepWithRusanov3DOnPatchInterior
"void" timeStepWithRusanov3DOnPatchInterior("::exahype2::CellData<double, double>&" patchData)
Rusanov update of the interior of a 3d patch.
Definition RusanovWithOptimisticExplicitEuler.py:324
exahype2.api.fv.RusanovWithOptimisticExplicitEuler.timeStepWithRusanov2DOnShell
"void" timeStepWithRusanov2DOnShell("::exahype2::CellData<double, double>&" patchData)
Rusanov update of the shell of a 2d patch.
Definition RusanovWithOptimisticExplicitEuler.py:126