This page was generated from unit-5.5-cuda/EulerEquations.ipynb.
5.5.3 Euler equations¶
state variables:
mass density \(\rho\)
momentum \(m = \rho u\), with velocity \(u\)
energy density \(E\)
conservation of mass, momentum and energy:
\begin{eqnarray*} \frac{\partial \rho}{\partial t} & = & -\operatorname{div} \rho u \\ \frac{\partial \rho u}{\partial t} & = & -\operatorname{div} (\rho u \otimes u + p I) \\ \frac{\partial E}{\partial t} & = & -\operatorname{div} (E+p) u \end{eqnarray*}
closure equation for pressure: \(p = (\gamma-1) (E - \rho \tfrac{1}{2} |u|^2)\)
with heat capacity ration \(\gamma\) depending on gas, \(\gamma=1.4\) for atmosphere.
Compact notation:
\[\frac{\partial U}{\partial t} + \operatorname{div} F(U) = 0\]
with vector of state variables (in \(R^{d+2}\)):
\[\begin{split}U = \left( \begin{array}{c}
\rho \\ m \\ E
\end{array} \right)\end{split}\]
and flux in \(R^{(d+2) \times d}\):
\[\begin{split}F = \left( \begin{array}{c}
\rho u \\
\rho u \otimes u + p I \\
(E + p) u
\end{array} \right)\end{split}\]
[1]:
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.comp import MFOpts
from ngsolve.gpu import * # registers cuda, metal, or the host reference device
from time import sleep, time
print ("gpu backend:", backend)
gpu backend: metal
[2]:
mesh = Mesh(unit_square.GenerateMesh(maxh=0.02))
order=3
fesT = L2(mesh, order=order)**4
feshat = FacetFESpace(mesh, order=order)**4
fes = fesT*feshat
gfu = GridFunction(fes)
rho0 = 1+1*exp(-400*( (x-0.5)**2 + (y-0.5)**2))
with TaskManager():
gfu.components[0].Set( (rho0, 0, 0, 1) )
[3]:
gamma = 1.4 # Gas constant
def Flux (U):
rho = U[0]
u = U[1:3]/rho
E = U[3]
p = (gamma-1)*(E-rho/2*(u*u))
return CF ( (rho*u, rho*OuterProduct(u,u)+p*Id(2),
(E+p)*u)).Reshape((4,2))
ngsglobals.msg_level = 0
n = specialcf.normal(2)
stab = 1
def NumFlux(u, uo):
return 0.5*(Flux(u)+Flux(uo))
[4]:
truecompile = False
(u,uhat), (v,vhat) = fes.TnT()
# bfa1 = BilinearForm(fes, nonassemble=True)
term1a = InnerProduct (NumFlux(u, 2*uhat-u), v.Operator("normal")).Compile(truecompile, wait=True)*dx(element_boundary=True)
term1b = (2*stab*v*(u-uhat)).Compile(truecompile, wait=True)*dx(element_vb=BND)
term2 = (-InnerProduct(Flux(u),Grad(v))).Compile(truecompile, wait=True)*dx
# fused matrix-free operator B^T D B, the flux F(u) is evaluated pointwise in the kernel
bfa = BilinearForm (term1a+term1b+term2, mf=MFOpts(nonlinear=True)).Assemble()
embT, embhat = fes.embeddings
resT, reshat = fes.restrictions
rangeT = fes.Range(0)
rangehat = fes.Range(1)
invm1 = embT@fesT.Mass(1).Inverse()@embT.T
# traceop = fesT.TraceOperator(feshat, average=True)
traceop = 0.5*fesT.TraceOperator(feshat, average=False)
traceop = traceop.CreateDeviceMatrix()
uT, vT = fesT.TnT()
with TaskManager():
invm = BilinearForm(uT*vT*dx, diagonal=True).Assemble().mat.Inverse()
invm_host = embT@invm@embT.T
invm = invm_host.CreateDeviceMatrix()
print(rangeT, rangehat)
[0,232480) [232480,373568)
[5]:
rho0 = 0.2+1*exp(-400*( (x-0.5)**2 + (y-0.5)**2))
with TaskManager():
gfu.components[0].Set( (rho0, 0, 0, rho0) )
gfubnd = GridFunction(fes)
gfubnd.components[1].vec.data = traceop * gfu.components[0].vec
Projector(fes.GetDofs(mesh.Boundaries(".*")), True).Project(gfubnd.vec)
gf_rho = gfu.components[0][0]
gf_u = gfu.components[0][1:3] / gf_rho
gf_E = gfu.components[0][3]
gf_p = (gamma-1)*(gf_E-gf_rho/2*(gf_u*gf_u)) # pressure
gf_a = sqrt(gamma*gf_p/gf_rho) # speed of sound
gf_M = Norm(gf_u) / gf_a # Mach number
print ("Density")
scene_rho = Draw (gf_rho, mesh, deformation=True)
print ("Velocity")
scene_u = Draw(gf_u, mesh, vectors={"grid_size":100})
print ("Mach number")
scene_M = Draw(Norm(gf_u)/gf_a, mesh)
t = 0
tend = 0.3
tau = 1e-4
i = 0
bfamat = bfa.mat.CreateDeviceMatrix()
vec = gfu.vec.CreateDeviceVector(copy=True)
vecbnd = gfubnd.vec.CreateDeviceVector(copy=True)
hv = gfu.vec.CreateDeviceVector()
with TaskManager():
while t < tend:
vec[rangehat] = traceop * vec[rangeT] + vecbnd[rangehat]
hv.data = bfamat * vec
vec -= tau * invm * hv
t += tau
i += 1
if i%20 == 0:
gfu.vec.data = vec
scene_rho.Redraw()
scene_u.Redraw()
scene_M.Redraw()
Density
Velocity
Mach number
[6]:
print (bfamat.GetOperatorInfo())
SumMatrix, h = 373568, w = 373568
SumMatrix, h = 373568, w = 373568
SumMatrix, h = 373568, w = 373568
SumMatrix, h = 373568, w = 373568
SumMatrix, h = 373568, w = 373568
GPU_BTDTBMatrix<float>, h = 373568, w = 373568
GPU_BTDTBMatrix<float>, h = 373568, w = 373568
GPU_BTDTBMatrix<float>, h = 373568, w = 373568
GPU_BTDTBMatrix<float>, h = 373568, w = 373568
GPU_BTDTBMatrix<float>, h = 373568, w = 373568
GPU_BTDTBMatrix<float>, h = 373568, w = 373568
Code generation for the device¶
the volume-part of the bilinear-form
\[\begin{split}\int_T F(u) \nabla v
\qquad
\text{with}
\qquad
F = \left( \begin{array}{c}
\rho u \\
\rho u \otimes u + p I \\
(E + p) u
\end{array} \right)\end{split}\]
is represented in NGSolve as a list of operations:
[7]:
print (term2)
Compiled CF:
Step 0: trial-function diffop = Id, dim=4
Step 1: ComponentCoefficientFunction 0
input: 0
Step 2: 1
Step 3: binary operation '/'
input: 2 1
Step 4: subtensor [ first: 1, num: ( 2,), dist: ( 1,) ], dim=2
input: 0
Step 5: scalar-vector multiply, dim=2
input: 3 4
Step 6: scalar-vector multiply, dim=2
input: 1 5
Step 7: reshape, dims = 2 x 1
input: 5
Step 8: reshape, dims = 1 x 2
input: 5
Step 9: matrix-matrix multiply, dims = 2 x 2
input: 7 8
Step 10: scalar-matrix multiply, dims = 2 x 2
input: 1 9
Step 11: ComponentCoefficientFunction 3
input: 0
Step 12: scale 0.5
input: 1
Step 13: innerproduct, fix size = 2
input: 5 5
Step 14: binary operation '*'
input: 12 13
Step 15: binary operation '-'
input: 11 14
Step 16: scale 0.4
input: 15
Step 17: Identity matrix, dims = 2 x 2
Step 18: scalar-matrix multiply, dims = 2 x 2
input: 16 17
Step 19: binary operation '+', dims = 2 x 2
input: 10 18
Step 20: binary operation '+'
input: 11 16
Step 21: scalar-vector multiply, dim=2
input: 20 5
Step 22: VectorialCoefficientFunction, dim=8
input: 6 19 21
Step 23: unary operation ' ', dim=8
input: 22
Step 24: reshape, dims = 4 x 2
input: 23
Step 25: test-function diffop = grad, dims = 4 x 2
Step 26: innerproduct, fix size = 8
input: 24 25
Step 27: scale -1
input: 26
VOL
to extract the flux, we differentiate w.r.t. the test-function
[8]:
print (term2[0].coef.Diff(Grad(v)))
Compiled CF:
Step 0: trial-function diffop = Id, dim=4
Step 1: ComponentCoefficientFunction 0
input: 0
Step 2: 1
Step 3: binary operation '/'
input: 2 1
Step 4: subtensor [ first: 1, num: ( 2,), dist: ( 1,) ], dim=2
input: 0
Step 5: scalar-vector multiply, dim=2
input: 3 4
Step 6: scalar-vector multiply, dim=2
input: 1 5
Step 7: reshape, dims = 2 x 1
input: 5
Step 8: reshape, dims = 1 x 2
input: 5
Step 9: matrix-matrix multiply, dims = 2 x 2
input: 7 8
Step 10: scalar-matrix multiply, dims = 2 x 2
input: 1 9
Step 11: ComponentCoefficientFunction 3
input: 0
Step 12: scale 0.5
input: 1
Step 13: innerproduct, fix size = 2
input: 5 5
Step 14: binary operation '*'
input: 12 13
Step 15: binary operation '-'
input: 11 14
Step 16: scale 0.4
input: 15
Step 17: Identity matrix, dims = 2 x 2
Step 18: scalar-matrix multiply, dims = 2 x 2
input: 16 17
Step 19: binary operation '+', dims = 2 x 2
input: 10 18
Step 20: binary operation '+'
input: 11 16
Step 21: scalar-vector multiply, dim=2
input: 20 5
Step 22: VectorialCoefficientFunction, dim=8
input: 6 19 21
Step 23: unary operation ' ', dim=8
input: 22
Step 24: reshape, dims = 4 x 2
input: 23
Step 25: scale -1, dims = 4 x 2
input: 24
And for this program we generate low-level code. It becomes the physics part of the fused matrix-free device kernel, which for every element gathers the dofs, evaluates the trial functions and the flux at the integration points, applies the test functions and scatters the result. The same source compiles for CUDA, Metal and the host reference device; here is the physics part of the generated kernel:
[9]:
import tempfile, os
kernelfile = os.path.join(tempfile.gettempdir(), "euler_kernel.txt")
BilinearForm (term2, mf=MFOpts(nonlinear=True, write_kernel=kernelfile)).Assemble().mat.CreateDeviceMatrix()
src = open(kernelfile).read()
print (src[src.find("// THE PHYSICS"):src.find("// END PHYSICS")])
// THE PHYSICS (WIP)
int i = 0;
{ // equation 0
auto values_0 = [&](int ip, int nr) { return xvals(0+nr); };
bool constexpr has_values_0 = true;
Real comp_0_0(0x0p+0 /* 0 */);
Real comp_0_1(0x0p+0 /* 0 */);
Real comp_0_2(0x0p+0 /* 0 */);
Real comp_0_3(0x0p+0 /* 0 */);
// step 0: trial-function diffop = Id
Real var_0_0(0x0p+0 /* 0 */);
var_0_0 = 0.0;
Real var_0_1(0x0p+0 /* 0 */);
var_0_1 = 0.0;
Real var_0_2(0x0p+0 /* 0 */);
var_0_2 = 0.0;
Real var_0_3(0x0p+0 /* 0 */);
var_0_3 = 0.0;
// if (tmp_0_0->fel) {
if (has_values_0) {
var_0_0 = values_0(i,0);
var_0_1 = values_0(i,1);
var_0_2 = values_0(i,2);
var_0_3 = values_0(i,3);
} else {
var_0_0 = comp_0_0;
var_0_1 = comp_0_1;
var_0_2 = comp_0_2;
var_0_3 = comp_0_3;
}
// step 1: ComponentCoefficientFunction 0
Real var_1;
var_1 = var_0_0;
// step 2: 1
Real var_2;
var_2 = 0x1p+0 /* 1 */;
// step 3: binary operation '/'
Real var_3;
var_3 = var_2 / var_1;
// step 4: subtensor [ first: 1, num: ( 2,), dist: ( 1,) ]
Real var_4_0;
Real var_4_1;
var_4_0 = var_0_1;
var_4_1 = var_0_2;
// step 5: scalar-vector multiply
Real var_5_0;
Real var_5_1;
var_5_0 = (var_3 * var_4_0);
var_5_1 = (var_3 * var_4_1);
// step 6: scalar-vector multiply
Real var_6_0;
Real var_6_1;
var_6_0 = (var_1 * var_5_0);
var_6_1 = (var_1 * var_5_1);
// step 7: reshape
Real var_7_0_0;
Real var_7_1_0;
var_7_0_0 = (var_5_0);
var_7_1_0 = (var_5_1);
// step 8: reshape
Real var_8_0_0;
Real var_8_0_1;
var_8_0_0 = (var_5_0);
var_8_0_1 = (var_5_1);
// step 9: matrix-matrix multiply
Real var_9_0_0;
Real var_9_0_1;
Real var_9_1_0;
Real var_9_1_1;
var_9_0_0 = ((var_7_0_0 * var_8_0_0));
var_9_0_1 = ((var_7_0_0 * var_8_0_1));
var_9_1_0 = ((var_7_1_0 * var_8_0_0));
var_9_1_1 = ((var_7_1_0 * var_8_0_1));
// step 10: scalar-matrix multiply
Real var_10_0_0;
Real var_10_0_1;
Real var_10_1_0;
Real var_10_1_1;
var_10_0_0 = (var_1 * var_9_0_0);
var_10_0_1 = (var_1 * var_9_0_1);
var_10_1_0 = (var_1 * var_9_1_0);
var_10_1_1 = (var_1 * var_9_1_1);
// step 11: ComponentCoefficientFunction 3
Real var_11;
var_11 = var_0_3;
// step 12: scale 0.5
Real var_12;
var_12 = (0x1p-1 /* 0.5 */ * var_1);
// step 13: innerproduct, fix size = 2
Real var_13;
var_13 = (((var_5_0 * var_5_0)) + (var_5_1 * var_5_1));
// step 14: binary operation '*'
Real var_14;
var_14 = var_12 * var_13;
// step 15: binary operation '-'
Real var_15;
var_15 = var_11 - var_14;
// step 16: scale 0.4
Real var_16;
var_16 = (0x1.9999999999998p-2 /* 0.4 */ * var_15);
// step 17: Identity matrix
Real var_17_0_0;
Real var_17_0_1;
Real var_17_1_0;
Real var_17_1_1;
var_17_0_0 = 1.0;
var_17_0_1 = 0.0;
var_17_1_0 = 0.0;
var_17_1_1 = 1.0;
// step 18: scalar-matrix multiply
Real var_18_0_0;
Real var_18_0_1;
Real var_18_1_0;
Real var_18_1_1;
var_18_0_0 = (var_16 * var_17_0_0);
var_18_0_1 = (var_16 * var_17_0_1);
var_18_1_0 = (var_16 * var_17_1_0);
var_18_1_1 = (var_16 * var_17_1_1);
// step 19: binary operation '+'
Real var_19_0_0;
Real var_19_0_1;
Real var_19_1_0;
Real var_19_1_1;
var_19_0_0 = var_10_0_0 + var_18_0_0;
var_19_0_1 = var_10_0_1 + var_18_0_1;
var_19_1_0 = var_10_1_0 + var_18_1_0;
var_19_1_1 = var_10_1_1 + var_18_1_1;
// step 20: binary operation '+'
Real var_20;
var_20 = var_11 + var_16;
// step 21: scalar-vector multiply
Real var_21_0;
Real var_21_1;
var_21_0 = (var_20 * var_5_0);
var_21_1 = (var_20 * var_5_1);
// step 22: VectorialCoefficientFunction
Real var_22_0;
Real var_22_1;
Real var_22_2;
Real var_22_3;
Real var_22_4;
Real var_22_5;
Real var_22_6;
Real var_22_7;
var_22_0 = var_6_0;
var_22_1 = var_6_1;
var_22_2 = var_19_0_0;
var_22_3 = var_19_0_1;
var_22_4 = var_19_1_0;
var_22_5 = var_19_1_1;
var_22_6 = var_21_0;
var_22_7 = var_21_1;
// step 23: unary operation ' '
Real var_23_0;
Real var_23_1;
Real var_23_2;
Real var_23_3;
Real var_23_4;
Real var_23_5;
Real var_23_6;
Real var_23_7;
var_23_0 = (var_22_0);
var_23_1 = (var_22_1);
var_23_2 = (var_22_2);
var_23_3 = (var_22_3);
var_23_4 = (var_22_4);
var_23_5 = (var_22_5);
var_23_6 = (var_22_6);
var_23_7 = (var_22_7);
// step 24: reshape
Real var_24_0_0;
Real var_24_0_1;
Real var_24_1_0;
Real var_24_1_1;
Real var_24_2_0;
Real var_24_2_1;
Real var_24_3_0;
Real var_24_3_1;
var_24_0_0 = (var_23_0);
var_24_0_1 = (var_23_1);
var_24_1_0 = (var_23_2);
var_24_1_1 = (var_23_3);
var_24_2_0 = (var_23_4);
var_24_2_1 = (var_23_5);
var_24_3_0 = (var_23_6);
var_24_3_1 = (var_23_7);
// step 25: scale -1
Real var_25_0_0;
Real var_25_0_1;
Real var_25_1_0;
Real var_25_1_1;
Real var_25_2_0;
Real var_25_2_1;
Real var_25_3_0;
Real var_25_3_1;
var_25_0_0 = (-0x1p+0 /* -1 */ * var_24_0_0);
var_25_0_1 = (-0x1p+0 /* -1 */ * var_24_0_1);
var_25_1_0 = (-0x1p+0 /* -1 */ * var_24_1_0);
var_25_1_1 = (-0x1p+0 /* -1 */ * var_24_1_1);
var_25_2_0 = (-0x1p+0 /* -1 */ * var_24_2_0);
var_25_2_1 = (-0x1p+0 /* -1 */ * var_24_2_1);
var_25_3_0 = (-0x1p+0 /* -1 */ * var_24_3_0);
var_25_3_1 = (-0x1p+0 /* -1 */ * var_24_3_1);
yvals(0) = var_25_0_0;
yvals(1) = var_25_0_1;
yvals(2) = var_25_1_0;
yvals(3) = var_25_1_1;
yvals(4) = var_25_2_0;
yvals(5) = var_25_2_1;
yvals(6) = var_25_3_0;
yvals(7) = var_25_3_1;
}
}