Peano
Loading...
Searching...
No Matches
MUSCLHancockDSLKernel3D.py
Go to the documentation of this file.
1# The following parameters shoule be set when calling Parser().parse(...)
2
3# Parameters used to generate non stateless code (set in parser)
4"""
5flux: "const fluxFunctor&",
6nonconservativeProduct: "const nonconservativeProductFunctor&",
7sourceTerm: "const sourceFunctor&",
8maxEigenvalue: "const maxEigenvalueFunctor&",
9"""
10
11def timeStepWithMUSCLHancock(patchData: "CellData<double, double>&") -> "void":
12 N = Integer("patchData.numberOfCells")
13 t = DataBlock([[0, N]], "patchData.t", "double")
14 dt = DataBlock([[0, N]], "patchData.dt", "double")
15 h = DataBlock([[0, N]], "patchData.cellSize", "tarch::la::Vector<3, double>")
16 x = DataBlock([[0, N]], "patchData.cellCentre", "tarch::la::Vector<3, double>")
17 maxEigenvalues = DataBlock([[0, N]], "patchData.maxEigenvalue", "double")
18
19 input = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [-H, NV + H], [-H, NV + H], [-H, NV + H], [0, N]], "patchData.QIn")
20 output = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV], [0, NV], [0, N]], "patchData.QOut")
21 output = input[:, 0:NV, 0:NV, 0:NV]
22
23 if EvaluateSource:
24 source = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [0, NV], [0, N]], eval("sourceTerm", input[:, 0:NV, 0:NV, 0:NV], x, h, t, dt))
25 output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] = output[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] + dt * source
26
27 #Compute slope with van Albada limiter
28 epsilon = 0.01
29 slopeRight = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV+2, 0:NV, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV+1, 0:NV, 0:NV]
30 slopeLeft = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV+1, 0:NV, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -2:NV, 0:NV, 0:NV]
31 slopeX = ((slopeRight * slopeRight + epsilon * epsilon)*slopeLeft + (slopeLeft * slopeLeft + epsilon * epsilon)*slopeRight) \
32 / (slopeLeft * slopeLeft + slopeRight * slopeRight + 2 * epsilon * epsilon)
33
34 slopeTop = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV+2, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV+1, 0:NV]
35 slopeBottom = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV+1, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -2:NV, 0:NV]
36 slopeY = ((slopeTop * slopeTop + epsilon * epsilon)*slopeBottom + (slopeBottom * slopeBottom + epsilon * epsilon)*slopeTop) \
37 / (slopeBottom * slopeBottom + slopeTop * slopeTop + 2 * epsilon * epsilon)
38
39 slopeFar = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, 0:NV+2] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, -1:NV+1]
40 slopeNear = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, -1:NV+1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, -2:NV]
41 slopeZ = ((slopeFar * slopeFar + epsilon * epsilon)*slopeNear + (slopeNear * slopeNear + epsilon * epsilon)*slopeFar) \
42 / (slopeNear * slopeNear + slopeFar * slopeFar + 2 * epsilon * epsilon)
43
44 #Compute face values
45 leftFace = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV+1, 0:NV, 0:NV] - 0.5 * slopeX
46 rightFace = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV+1, 0:NV, 0:NV] + 0.5 * slopeX
47 bottomFace = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV+1, 0:NV] - 0.5 * slopeY
48 topFace = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV+1, 0:NV] + 0.5 * slopeY
49 nearFace = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, -1:NV+1] - 0.5 * slopeZ
50 farFace = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, -1:NV+1] + 0.5 * slopeZ
51
52 #Evolve face values by half time step
53 if EvaluateFlux:
54 fluxLeft = DataBlock([[0, NumberOfUnknowns], [-1, NV+1], [0, NV], [0, NV], [0, N]], eval("flux", leftFace, x, h, t, dt, 0))
55 fluxRight = DataBlock([[0, NumberOfUnknowns], [-1, NV+1], [0, NV], [0, NV], [0, N]], eval("flux", rightFace, x, h, t, dt, 0))
56 fluxBottom = DataBlock([[0, NumberOfUnknowns], [0, NV], [-1, NV+1], [0, NV], [0, N]], eval("flux", bottomFace, x, h, t, dt, 1))
57 fluxTop = DataBlock([[0, NumberOfUnknowns], [0, NV], [-1, NV+1], [0, NV], [0, N]], eval("flux", topFace, x, h, t, dt, 1))
58 fluxNear = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [-1, NV+1], [0, N]], eval("flux", nearFace, x, h, t, dt, 2))
59 fluxFar = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [-1, NV+1], [0, N]], eval("flux", farFace, x, h, t, dt, 2))
60
61 leftFace[0:NumberOfUnknowns, -1:NV+1, 0:NV, 0:NV] = leftFace[0:NumberOfUnknowns, -1:NV+1, 0:NV, 0:NV] \
62 - 0.5 * (dt / h(0)) * (fluxRight - fluxLeft)
63 rightFace[0:NumberOfUnknowns, -1:NV+1, 0:NV, 0:NV] = rightFace[0:NumberOfUnknowns, -1:NV+1, 0:NV, 0:NV] \
64 - 0.5 * (dt / h(0)) * (fluxRight - fluxLeft)
65 bottomFace[0:NumberOfUnknowns, 0:NV, -1:NV+1, 0:NV] = bottomFace[0:NumberOfUnknowns, 0:NV, -1:NV+1, 0:NV] \
66 - 0.5 * (dt / h(1)) * (fluxTop - fluxBottom)
67 topFace[0:NumberOfUnknowns, 0:NV, -1:NV+1, 0:NV] = topFace[0:NumberOfUnknowns, 0:NV, -1:NV+1, 0:NV] \
68 - 0.5 * (dt / h(1)) * (fluxTop - fluxBottom)
69 nearFace[0:NumberOfUnknowns, 0:NV, 0:NV, -1:NV+1] = nearFace[0:NumberOfUnknowns, 0:NV, 0:NV, -1:NV+1] \
70 - 0.5 * (dt / h(2)) * (fluxFar - fluxNear)
71 farFace[0:NumberOfUnknowns, 0:NV, 0:NV, -1:NV+1] = farFace[0:NumberOfUnknowns, 0:NV, 0:NV, -1:NV+1] \
72 - 0.5 * (dt / h(2)) * (fluxFar - fluxNear)
73
74 if EvaluateSource:
75 sourceLeftFace = DataBlock([[0, NumberOfUnknowns], [-1, NV+1], [0, NV], [0, NV], [0, N]], eval("sourceTerm", leftFace, x, h, t, dt))
76 sourceRightFace = DataBlock([[0, NumberOfUnknowns], [-1, NV+1], [0, NV], [0, NV], [0, N]], eval("sourceTerm", rightFace, x, h, t, dt))
77 sourceBottomFace = DataBlock([[0, NumberOfUnknowns], [0, NV], [-1, NV+1], [0, NV], [0, N]], eval("sourceTerm", bottomFace, x, h, t, dt))
78 sourceTopFace = DataBlock([[0, NumberOfUnknowns], [0, NV], [-1, NV+1], [0, NV], [0, N]], eval("sourceTerm", topFace, x, h, t, dt))
79 sourceNearFace = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [-1, NV+1], [0, N]], eval("sourceTerm", nearFace, x, h, t, dt))
80 sourceFarFace = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [-1, NV+1], [0, N]], eval("sourceTerm", farFace, x, h, t, dt))
81
82 leftFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] = leftFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] + 0.5 * dt * sourceLeftFace
83 rightFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] = rightFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] + 0.5 * dt * sourceRightFace
84 bottomFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] = bottomFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] + 0.5 * dt * sourceBottomFace
85 topFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] = topFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] + 0.5 * dt * sourceTopFace
86 nearFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] = nearFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] + 0.5 * dt * sourceNearFace
87 farFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] = farFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] + 0.5 * dt * sourceFarFace
88
89 #Evolve volumes
90 maxEigenvaluesLeft = DataBlock([[0, 1], [-1, NV+1], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", leftFace, x, h, t, dt, 0))
91 maxEigenvaluesRight = DataBlock([[0, 1], [-1, NV+1], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", rightFace, x, h, t, dt, 0))
92 lambdaX = max(maxEigenvaluesLeft[0:1, 0:NV+1, 0:NV, 0:NV], maxEigenvaluesRight[0:1, -1:NV, 0:NV, 0:NV])
93
94 maxEigenvaluesBottom = DataBlock([[0, 1], [0, NV], [-1, NV+1], [0, NV], [0, N]], eval("maxEigenvalue", bottomFace, x, h, t, dt, 1))
95 maxEigenvaluesTop = DataBlock([[0, 1], [0, NV], [-1, NV+1], [0, NV], [0, N]], eval("maxEigenvalue", topFace, x, h, t, dt, 1))
96 lambdaY = max(maxEigenvaluesBottom[0:1, 0:NV, 0:NV+1, 0:NV], maxEigenvaluesTop[0:1, 0:NV, -1:NV, 0:NV])
97
98 maxEigenvaluesNear = DataBlock([[0, 1], [0, NV], [0, NV], [-1, NV+1], [0, N]], eval("maxEigenvalue", nearFace, x, h, t, dt, 2))
99 maxEigenvaluesFar = DataBlock([[0, 1], [0, NV], [0, NV], [-1, NV+1], [0, N]], eval("maxEigenvalue", farFace, x, h, t, dt, 2))
100 lambdaZ = max(maxEigenvaluesNear[0:1, 0:NV, 0:NV, 0:NV+1], maxEigenvaluesFar[0:1, 0:NV, 0:NV, -1:NV])
101
102 rusanovFluxLeft = - 0.5 * lambdaX[0:1, 0:NV, 0:NV, 0:NV] * (leftFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] - rightFace[0:NumberOfUnknowns, -1:NV-1, 0:NV, 0:NV])
103 rusanovFluxRight = - 0.5 * lambdaX[0:1, 1:NV+1, 0:NV, 0:NV] * (leftFace[0:NumberOfUnknowns, 1:NV+1, 0:NV, 0:NV] - rightFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV])
104 rusanovFluxBottom = - 0.5 * lambdaY[0:1, 0:NV, 0:NV, 0:NV] * (bottomFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] - topFace[0:NumberOfUnknowns, 0:NV, -1:NV-1, 0:NV])
105 rusanovFluxTop = - 0.5 * lambdaY[0:1, 0:NV, 1:NV+1, 0:NV] * (bottomFace[0:NumberOfUnknowns, 0:NV, 1:NV+1, 0:NV] - topFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV])
106 rusanovFluxNear = - 0.5 * lambdaZ[0:1, 0:NV, 0:NV, 0:NV] * (nearFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] - farFace[0:NumberOfUnknowns, 0:NV, 0:NV, -1:NV-1])
107 rusanovFluxFar = - 0.5 * lambdaZ[0:1, 0:NV, 0:NV, 1:NV+1] * (nearFace[0:NumberOfUnknowns, 0:NV, 0:NV, 1:NV+1] - farFace[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV])
108
109 if EvaluateFlux:
110 timeAveragedfluxLeft = DataBlock([[0, NumberOfUnknowns], [-1, NV + 1], [0, NV], [0, NV], [0, N]], eval("flux", leftFace, x, h, t, dt, 0))
111 timeAveragedfluxRight = DataBlock([[0, NumberOfUnknowns], [-1, NV + 1], [0, NV], [0, NV], [0, N]], eval("flux", rightFace, x, h, t, dt, 0))
112 timeAveragedfluxBottom = DataBlock([[0, NumberOfUnknowns], [0, NV], [-1, NV + 1], [0, NV], [0, N]], eval("flux", bottomFace, x, h, t, dt, 1))
113 timeAveragedfluxTop = DataBlock([[0, NumberOfUnknowns], [0, NV], [-1, NV + 1], [0, NV], [0, N]], eval("flux", topFace, x, h, t, dt, 1))
114 timeAveragedfluxNear = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [-1, NV + 1], [0, N]], eval("flux", nearFace, x, h, t, dt, 2))
115 timeAveragedfluxFar = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [-1, NV + 1], [0, N]], eval("flux", farFace, x, h, t, dt, 2))
116
117 rusanovFluxLeft = rusanovFluxLeft + 0.5 * (timeAveragedfluxLeft[:, 0:NV, 0:NV, 0:NV] + timeAveragedfluxRight[:, -1:NV-1, 0:NV, 0:NV])
118 rusanovFluxRight = rusanovFluxRight + 0.5 * (timeAveragedfluxLeft[:, 1:NV+1, 0:NV, 0:NV] + timeAveragedfluxRight[:, 0:NV, 0:NV, 0:NV])
119 rusanovFluxBottom = rusanovFluxBottom + 0.5 * (timeAveragedfluxBottom[:, 0:NV, 0:NV, 0:NV] + timeAveragedfluxTop[:, 0:NV, -1:NV-1, 0:NV])
120 rusanovFluxTop = rusanovFluxTop + 0.5 * (timeAveragedfluxBottom[:, 0:NV, 1:NV+1, 0:NV] + timeAveragedfluxTop[:, 0:NV, 0:NV, 0:NV])
121 rusanovFluxNear = rusanovFluxNear + 0.5 * (timeAveragedfluxNear[:, 0:NV, 0:NV, 0:NV] + timeAveragedfluxFar[:, 0:NV, 0:NV, -1:NV-1])
122 rusanovFluxFar = rusanovFluxFar + 0.5 * (timeAveragedfluxNear[:, 0:NV, 0:NV, 1:NV+1] + timeAveragedfluxFar[:, 0:NV, 0:NV, 0:NV])
123
124 output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + (dt / h(0)) * NV * (rusanovFluxLeft - rusanovFluxRight)
125 output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + (dt / h(1)) * NV * (rusanovFluxBottom - rusanovFluxTop)
126 output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + (dt / h(2)) * NV * (rusanovFluxNear - rusanovFluxFar)
127
128 #if EvaluateNonconservativeProduct:
129 # average_x = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV, 0:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV + 1, 0:NV])
130 # delta_x = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV + 1, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV, 0:NV]
131 # nonconservative_x = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV + 1], [0, NV], [0, N]], eval("nonconservativeProduct", average_x, delta_x, x, h, t, dt, 0))
132 # output = output + 0.5 * NV * dt / h(0) * (-nonconservative_x[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV] - nonconservative_x[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 1:NV+1, 0:NV])
133 #
134 # average_y = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV + 1])
135 # delta_y = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV + 1] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV]
136 # nonconservative_y = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV + 1], [0, N]], eval("nonconservativeProduct", average_y, delta_y, x, h, t, dt, 1))
137 # output = output + 0.5 * NV * dt / h(1) * (-nonconservative_y[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV] - nonconservative_y[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 1:NV+1])
138
139 if EvaluateMaximumEigenvalueAfterTimeStep:
140 maxEigenvaluesX2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, h, t, dt, 0))
141 maxEigenvaluesY2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, h, t, dt, 1))
142 maxEigenvaluesZ2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, h, t, dt, 2))
143 maxEigenvaluesXY = max(maxEigenvaluesX2, maxEigenvaluesY2)
144 maxEigenvaluesPerVolume = max(maxEigenvaluesXY, maxEigenvaluesZ2)
145 maxEigenvalues = max(maxEigenvaluesPerVolume)