The form of the equation we want to solve is
u,t = gi(u; {h(u),x}),xi ,
where the comma is used to denote partial differentiation
and for lists we use a semicolon
and where we suppress all state-space indices.
To simplify the theory we write
gi(u; (h(u),x)) = gi(u; (h,u.u,x)) =: fi(u; u,x),
where f might be an ugly expression.
So for theoretical purposes we can take the general equation to be
u,t = fi(u; (u,x)),xi.
We use the chain rule to decompose this
into advective and diffusive parts:
u,t = fi,u.u,xi + fi,(u,xj).u,{xj xi}.
Define
Ai := fi,u and
Bij := fi,(u,xj).
Adopt the shorthand
u,i = u,xi and
u,ij = u,{xi xj}.
We seek generally necessary conditions for well-posedness,
so for the purposes of studying well-posedness freeze these
matrices of coefficients. We now have the system
u,t = Ai.u,i + Bij.u,ij
This can only be well-posed if it is so for data that is
homogenous perpendicular to a direction n. So assume
u(x,t) = u(x•n⊗n,t) = u(s,t) where s := x•n (so s,x = n).
Then
u,x = u,s*s,x = u,s*n and
u,xx = u,ss*n⊗n.
So we get the 1d system
u,t = ni*Ai.u,s + ni*nj*Bij.u,ss.
(Using n and m for state-space indices this writes as
u_n,t = n_i*A_nmi*u_m,s + n_i*n_j*B_nmij*u_m,ss.)
Avoiding all indices this writes as
u,t = (n•A).u,s + (n•B•n).u,ss,
where I use "•" to indicate contraction of
spacial indices and "." to indicate contraction of
state-space indices.
For this to be well-posed we need
eigenvalues of n•A are real for all n and
eigenvalues of n•B•n are positive for all n.
So we see that for any direction n in physical space
there is a set of directions in state space
(1) in which waves propagate each with their own
speed (as defined by the eigenstructure of n•A, and
(2) in which diffusion occurs independently, each
with its own direction and rate (as defined by the
eigenstructure of n•B•n.
In the scalar case one need merely check that n•B•n
is always positive, i.e. that B is positive definite. Is
there a way to check that n•B•n has positive
eigenvalues for all n without looking at every n? How can the
eigenstructure vary with n? Can a convexity argument be used?
Scalar case.
(In this section all dot products denote contraction
over physical space.) For scalar u, in the equation
u,t = Ai*u,i + Bij*u,ij = A.u,x + B:u,xx
A is an advection velocity and B is a 3x3 WLOG symmetric matrix
whose eigenstructure reveals rates of diffusion in orthogonal
directions.
Indeed, transforming variables by x' = L.x, where L
is an orthonormal matrix of left eigenvectors with transpose LT,
B:u,xx = Bij*u,{xi xj}
= (xm'_,xi)*Bij*(xn'_,xj)*u_,{xm' xn'}
= L•B•LT:u_,{x'x'}
That is,
B:u,xx = (Dx.B.Dx)u = (Dx'.(L.B.LT).Dx')u = V:u_,{x'x'}
where the diagonal matrix V := L.B.LT should have nonnegative entries.
This raises the question:
* Is there any natural extension to systems of the notion of
orthogonal directions of diffusion in physical space?
Scalar heat equation.
Consider the scalar heat equation:
u,t = a^2*u,xx = u,yy where y = x/a.
here a is the characteristic length scale
of diffusion in one unit of time.
...
Consider a generic scalar diffusion equation with constant coefficients:
u,t = B:u,xx
Transform physical space by x' = L•x.
We have seen that
u,t = V:u_,{x'x'} where V := L.B.LT.
Let B = diag(λ1,λ2) and let L be a rotation matrix
L = [c -s] where c=cos(theta) and s=sin(theta).
[s c]
So
V = [c -s] [λ1 0] [ c s]
[s c] [0 λ2] [-s c]
= [c -s] [λ1c λ1s]
[s c] [-λ2s λ2c]
= [λ1cc+λ2ss (λ1-λ2)cs]
[(λ1-λ2)cs λ1ss+λ2cc]
Choose theta = pi/4. Then c = s = 1/sqrt(2) and
V = [(λ1+λ2)/2 (λ1-λ2)/2] = [a c]
[(λ1-λ2)/2 (λ1+λ2)/2] [c a]
where a=(λ1+λ2)/2 and c=(λ1-λ2)/2.
Note that the nonnegative definiteness conditions
λ1≥0 and λ2≥0
become the conditions
a ≥ 0 and |a| ≥ |c|.
In other words, rotating the domain of the diffusion equation
u,t = λ1*u,xx + λ2*u,yy
by 45 degrees gives a diffusion equation with mixed derivatives
u,t = a u,xx + 2c u,xy + a*u,yy
where the eigenvalues of diffusion are
λ1= a+c and
λ2= a-c.
In practice the form of diffusive terms looks like
u,t = Div(A:Grad h),
where u is a tuple of conserved quantities (e.g. mass,
momentum, and energy) which specifies the state, comma denotes
partial differentation, A is a "fourth-order" tensor (with
two indices for physical space, a (suppressed) index for the
state-space of u, and an index for the state-space of h)
comprised of closure coefficients and h is a function of u
(e.g. the temperature or the velocity). Physical collision
models give A(u) and anomalous models give A(t,x,u).
We make the definition
g(t; x; u; Grad h) = A : Grad h
In components we would write
u,t = (A_ijk*h_j,k),i = g_i,i, where
g_i := A_ijk*h_j,k.
In case the coefficients A are constant,
we can pull out the derivative and write
u,t = Div(Div(A.h)) = Div(Div(F)), where
F := A.h.
To derive the discontinuous galerkin form
we first rewrite the differential equation as a system:
--- generic --- --- constant-coefficient ---
u,t = Div g u,t = Div(v)
g = A:w v = Div F
w = Grad h F = A.h.
We then multiply by a test function and integrate
by parts to get:
--- generic ---
int(u*phi),t = oint(n.g*phi) - int(g.Grad phi)
g(u,w) = A:w
int(w*phi) = oint(n*h*phi) - int(h*Grad phi)
--- constant-coefficient ---
int(u*phi),t = oint(n.v*phi) - int(v.Grad phi)
int(v*phi) = oint(n.F*phi) - int(F.Grad phi)
F(u) = A.h(u)
In the generic case the user must supply two
callback functions defining g and h,
whereas in the constant-coefficient case the user
must supply a single callback function defining F.
Each use of a callback function involves evaluating
the arguments at gaussian quadrature points.
The question is how to evaluate the boundary fluxes.
I think that for an unstructured mesh we could just
take the average from both sides of the boundary,
but for repeated derivatives in a cartesian mesh
I believe that a singularity arises in the case
of taking the average which causes loss of an
order of accuracy for methods that would normally
give an odd order of accuracy.
As a guiding principle we seek a method which agrees
with best methods for the heat equation.
--- correspondence between generic and constant-coefficient cases ---
Remark: We can tighten the correspondence between the
generic and constant-coefficient cases when we can write
A(u) = a(u)*C,
where a is a scalar and C is tensor of integers. Then
u,t = Div(a*C:Grad h)
= Div(a*Div(C.h))
=F=
====v===
=====g====
This framework fits the (10-moment) heat flux and the contribution
of the stress tensor to the momentum equation but not the contribution
of the stress tensor to the energy equation.
--- overview --- --- generic case --- u,t = (Axx.(h,x)),x + (Axy.(h,y)),x + (Ayx.(h,x)),y + (Ayy.(h,y)),y = (Axx.(wxx)),x + (Axy.(wxy)),x + (Ayx.(wyx)),y + (Ayy.(wyy)),y = ( gxx ),x + ( gxy ),x + ( gyx ),y + ( gyy ),y --- constant-coefficient case --- u,t = ((Axx.h),x),x + ((Axy.h),y),x + ((Ayx.h),x),y + ((Ayy.h),y),y = ( Fxx ,x),x + ( Fxy ,y),x + ( Fyx ,x),y + ( Fyy ,y),y = ( vxx ),x + ( vxy ),x + ( vyx ),y + ( vyy ),y Note that wij = h,j and gij = Aij.wij whereas Fij = Aij.h and vij = Fij,j So as quantities gij=vij, but I think of gij as a callback function. Conceptually, the algorithms are: --- generic --- w = ComputeGrad(h,u) u,t = ComputeDiv(g,u,w) --- constant-coefficient --- u,t = ComputeDivDiv(F,u) --- constant-coefficient case --- In two dimensions the constant-coefficient diffusive system is u,t = (Axx.h),xx + (Axy.h),yx + (Ayx.h),xy + (Ayy.h),yy Making auxiliary definitions: --- constant-coefficient --- u,t = vx,x + vy,y vx = Fxx,x + Fxy,y vy = Fyx,x + Fyy,y Fxx = Axx.h Fxy = Axy.h Fyx = Ayx.h Fyy = Ayy.h In order to use one-sided boundary values for repeated derivatives and two-sided averages for mixed derivatives, we rewrite this as --- constant-coefficient --- u,t = vxx,x + vxy,x + vyx,y + vyy,y vxx = Fxx,x [boundary flux: use low-sided] vxy = Fxy,y [boundary flux: use averaged] vyx = Fyx,x [boundary flux: use averaged] vyy = Fyy,y [boundary flux: use low-sided] Fxx = Axx.h [boundary flux: use high-sided] Fxy = Axy.h [boundary flux: use averaged] Fyx = Ayx.h [boundary flux: use averaged] Fyy = Ayy.h [boundary flux: use high-sided] Multiplying by a test function and integrating by parts gives int(vxx*phi) = -iint(Fxx*phi,x) + oint(nx*Fxx*phi) [high-sided], int(vyx*phi) = -iint(Fyx*phi,x) + oint(nx*Fyx*phi) [averaged], int(vxy*phi) = -iint(Fxy*phi,y) + oint(ny*Fxy*phi) [averaged], int(vyy*phi) = -iint(Fyy*phi,y) + oint(ny*Fyy*phi) [high-sided], int(u*phi),t = -iint(vxx*phi,x) + oint(nx*vxx*phi) [low-sided] -iint(vxy*phi,x) + oint(nx*vxy*phi) [averaged] -iint(vyx*phi,y) + oint(ny*vyx*phi) [averaged] -iint(vyy*phi,y) + oint(ny*vyy*phi) [low-sided]. Here nx is the unit vector in the x direction and ny is the unit vector in the y direction. Following James's precedent with the heat equation, for the vxx and fxx terms and for the vyy and fyy terms I use alternating one-sided differences, whereas for the vxy and fxy terms and for the vyx and fyx terms I used a two-sided average. --- generic case --- In two dimensions the generic diffusive system is u,t = (Axx.(h,x)),x + (Axy.(h,y)),x + (Ayx.(h,x)),y + (Ayy.(h,y)),y Making auxiliary definitions, === generic === u,t = gx,x + gy,y gx = Axx.wx + Axy.wy gy = Ayx.wx + Ayy.wy wx = h,x wy = h,y In order to use one-sided boundary values for repeated derivatives and two-sided averages for mixed derivatives, we rewrite this as === generic === u,t = gxx,x + gy,y = gxx,x + gxy,x + gyx,y + gyy,y gx = gxx + gxy gy = gyx + gyy gxx = Axx.wxx [use low-sided] gxy = Axy.wxy [use averaged] gyx = Ayx.wyx [use averaged] gyy = Ayy.wyy [use low-sided] wxx = h,x [use high-sided] wyx = h,x [use averaged] wxy = h,y [use averaged] wyy = h,y [use high-sided] --- algorithm --- [wxx,wyx,wxy,wyy] = ComputeGrad(h,u) u,t = ComputeDiv(gxx,gxy,gyx,gyy,u,wxx,wyx,wxy,wyy) Multiplying by a test function and integrating by parts gives int(wxx*phi) = -iint(h*phi,x) + oint(nx*h*phi) [flux: use high-sided], int(wyx*phi) = -iint(h*phi,x) + oint(nx*h*phi) [flux: use averaged], int(wxy*phi) = -iint(h*phi,y) + oint(ny*h*phi) [flux: use averaged], int(wyy*phi) = -iint(h*phi,y) + oint(ny*h*phi) [flux: use high-sided], int(u*phi),t = -iint(gxx*phi,x) + oint(nx*gxx*phi) [flux: use low-sided] -iint(gxy*phi,x) + oint(nx*gxy*phi) [flux: use averaged] -iint(gyx*phi,y) + oint(ny*gyx*phi) [flux: use averaged] -iint(gyy*phi,y) + oint(ny*gyy*phi) [flux: use low-sided]. I'm not 100% sure that taking one-sided differences will give the desired order of accuracy in the non-constant-coefficient case. As you refine the mesh the amount of change in the coefficient values from one side of a cell to the other should decrease like an order-1 function. So my intuition is that for low-order methods you will still get the desired order of accuracy, but maybe not for higher order. If we were using an unstructured mesh I think we could get away with always taking the average. Details for implementation: int(wxx*phi) = oint(nx*h*phi) - iint(h*phi,x) IntEdgeXlow(h,u) L2ProjectDx(h,u) int(wyx*phi) = oint(nx*h*phi) - iint(h*phi,x) IntEdgeXave(h,u) L2ProjectDx(h,u) int(wxy*phi) = oint(ny*h*phi) - iint(h*phi,y) IntEdgeYave(h,u) L2ProjectDy(h,u) int(wyy*phi) = oint(ny*h*phi) - iint(h*phi,y) IntEdgeYlow(h,u) L2ProjectDy(h,u) int(u*phi),t = oint(nx*gxx*phi) - iint(gxx*phi,x) IntEdgeXhgh2(gxx,u,wxx) L2ProjectGrad2(gxx,u,wxx) + oint(nx*gxy*phi) - iint(gxy*phi,x) IntEdgeXave2(gxy,u,wxy) L2ProjectGrad2(gxy,u,wxy) + oint(ny*gyx*phi) - iint(gyx*phi,y) IntEdgeYave2(gyx,u,wyy) L2ProjectGrad2(gyx,u,wyx) + oint(ny*gyy*phi) - iint(gyy*phi,y) IntEdgeYhgh2(gyy,u,wyy) L2ProjectGrad2(gyy,u,wyy) where h and gij are function pointers, and u and wij represent stored data. An algorithm to implement this is then: * Evaluate u and h(u) at quadrature points. Could store this data to avoid repeated evaluation. * Use guassian quadrature and appropriate fluxes to evaluate integrals yielding polynomial representations of wxx, wyx, wxy, and wyy. Store this data and free previous data. * Evaluate gxx(u,wxx), gyx(u,wyx), gxy(u,wxy), and gyy(u,wyy) at quadrature points using appropriate fluxes. (In the constant-coefficient case these evaluations are not necessary; gij is a linear function of wij,i and it suffices to .) * Use guassian quadrature and appropriate fluxes to evaluate integrals yielding int(u*phi),t; add this to Lstar and free previous data. Potentially conflicting goals of implementation: * Avoid unnecessary evaluation. + this goal often conflicts with next three. * Avoid unnecessary storage. * Avoid complexity. * Avoid sacrificing generality to accomodate specific details. Data representation constraints: * For differentiation you need the polynomial representation. * For arithmetic (except addition or subtraction) you need quadrature point sample representation. Implementing mechanisms and issues: * L2Project takes a polynomial representation, evaluates an arbitrary function at quadrature points, and returns a polynomial representation. It involves no boundary evaluations since it involves no derivatives. It is used to handle source terms, not differentiated terms. * To handle differentiated terms, ConstructL uses L2ProjectGrad and RiemannSolve (which evaluates boundary states, computes an intermediate state, and integrates its contribution to each basis function) * The idea of L2Project* is to find the polynomial representation closest to the data set. * To obtain a polynomial representation of w we need something like L2ProjectD[xy] and something like IntegrateEdge.
Recall that we are trying to solve the system
u,t = (Axx.wx + Axy.wy),x + (Ayx.wx + Ayy.wy),y,
wx = h(u),x
wy = h(u),y
In the hopes of getting a hyperbolic system we
modify the last two equations:
u,t = (Axx.wx + Axy.wy),x + (Ayx.wx + Ayy.wy),y,
eps^2*wx,t + h(u),x = wx
eps^2*wy,t + h(u),y = wy
As eps goes to zero, this should approximate the diffusion
equation. (But wave speeds are on the order of 1/eps and should
go to infinity.) In particular, for the special case h(u) = u
and Aij = identity tensor this gives a hyperbolic system with
wave speed 1/eps.
What are the actual equations that we need to solve? + full 5-moment closed system The five-moment equations (Navier-Stokes) have two diffusive terms: one for heat flux and one for viscosity. rho,t + Div(rho*u) = 0, (rho*u),t + Div(rho*u*u)+Grad p = Div(vstress) + drag + (q/m)*rho*(E + u cross B) energy,t + Div(u*(energy+p)) + Div(heatFlux) = Div(vstress.u) + (q/m)rho*u.E + u.drag + Qf + Qt, where E is electric field, B is magnetic field, p = pressure = (2/3)*(energy-rho*u.u/2), vstress = viscousStress = 2*viscosity*(Sym(Grad u) - idtens*Div(u)/3), drag = dragCoefficient*rho*rho_p(u_p - u), p = other species (when used as an index), heatFlux = -heatConductivity*rho*Grad T, T = temperature = p/n, n = numberDensity = rho/m, Qf = frictional heating = drag.(u_p-u)*(m_p/(m_p+m)), and Qt = thermalEquilibration = thermalConductivity*rho*rho_p*(T_p-T), where viscosity, dragCoefficient, heatConductivity, and thermalConductivity are all parameters potentially dependent on the state. Note that all coupling to the electromagnetic field and to other species is through nondifferentiated source terms. + diffusive 5-moment closed system Neglecting all nondiffusive terms, the system reduces to: (rho*u),t = Div(vstress) energy,t + Div(heatFlux) = Div(vstress.u) where again vstress = viscousStress = 2*viscosity*(Sym(Grad u) - idtens*Div(u)/3), heatFlux = -heatConductivity*rho*Grad T, + full 10-moment closed system The ten-moment equations have one diffusive term, for heat flux. They are: rho,t + Div(rho*u) = 0, (rho*u),t + Div(Energy) = drag + (q/m)*rho*(E + u cross B), Energy,t + 3*Div Sym(u*Energy) - 2*Div(rho*u*u*u) + Div HeatFlux = (q/m)*2*Sym(rho*u*E + Energy cross B) + 2*Sym(u*drag) + R + QQf + QQt, where rho is mass density, u is fluid velocity, Energy is the energy tensor, * denotes scalar or tensor product, . denotes contracted (dot) product, Sym symmetrizes its argument tensor, E is electric field, B is magnetic field, n = numberDensity = rho/m, P = pressureTensor = Energy-rho*u*u, p = pressure = trace(P)/3, TT = temperatureTensor = P/n, T = temperature = p/n = trace(TT)/3, p = other species (when used as an index), w_p = u_p - u = drift velocity of species p, drag = dragCoefficient*rho*rho_p*w_p, idtens = identity tensor, Tinv = TT^{-1} = inverse of temperature tensor, HeatFlux = a_1*3*Sym(Grad Tinv) + a_0*3*Sym(idtens*trace(3*Sym(Grad Tinv))), R = isotropizationTensor = (p*idtens - P)/isoPeriod, QQf = frictional heating tensor = (2*m_p/(m_p+m))*drag.((a_para-a_perp)*idtens*w_p+a_perp*w_p*idtens), QQt = thermalEquilibration = (2/3)*tho*thermalConductivity*rho_p*(TT_p-TT), where isoPeriod, dragCoefficient, a_0, a_1, thermalConductivity, and a_perp are positive parameters determined by a collision integral, 1 ≥ 2*a_perp, and a_para = 1 - 2*a_perp. and where these coefficients are related to the coefficients of the 5-moment closure by viscosity = p*isoPeriod, heatConductivity = (25*a_0 + 5*a_1)/(2*T*T). + diffusive 10-moment closed system Neglecting all nondiffusive terms, the system reduces to: Energy,t = - Div HeatFlux. where again HeatFlux = a_1*3*Sym(Grad Tinv) + a_0*3*Sym(idtens*trace(3*Sym(Grad Tinv))), Tinv = TT^{-1} = inverse of temperature tensor, TT = temperatureTensor = P/n, P = pressureTensor = Energy-rho*u*u, n = numberDensity = rho/m, and where heatConductivity = (25*a_0 + 5*a_1)/(2*T*T).
We now want to write out the diffusive terms in coordinates. + 5-moment diffusive terms Recall the diffusive 5-moment system: (rho*u),t = Div(vstress) energy,t + Div(heatFlux) = Div(vstress.u) where vstress = viscousStress = 2*mu*(Sym(Grad u) - idtens*Div(u)/3), heatFlux = -k*Grad T, where mu = viscosity, k = heatConductivity*rho. In components, vstress_ij = mu*(u_j,i + u_i,j - u_k,k*(2/3)*delta_ij). In two dimensions, writing out the definition of the stress and the momentum equation in components gives vstress_xx/mu = (4/3)*u_x,x - (2/3)*u_y,y vstress_xy/mu = u_y,x + u_x,y vstress_xz/mu = u_z,x vstress_yx/mu = u_y,x + u_x,y vstress_yy/mu =(-2/3)*u_x,x + (4/3)*u_y,y vstress_yz/mu = u_z,y and (rho*u_x),t = vstress_xx,x + vstress_yx,y (rho*u_y),t = vstress_xy,x + vstress_yy,y (rho*u_z),t = vstress_xz,x + vstress_yz,y In case mu is constant we can handle the momentum equation with a constant-coefficient term, but the energy equation intrinsically requires the more generic framework. The energy equation reads energy,t = Div(vstress.u) - Div(heatFlux) or in coordinates energy,t = (u_j*vstress_ij),i - heatFlux_i,i. In two dimensions, energy,t = (u_x*vstress_xx),x + (u_x*vstress_yx),y + (u_y*vstress_xy),x + (u_y*vstress_yy),y + (u_z*vstress_xz),x + (u_z*vstress_yz),y - heatFlux_x,x - heatFlux_y,y. To handle the heat flux, recall that heatFlux = -k*Grad T. That is, heatFlux_x = -k*T,x and heatFlux_y = -k*T,y. So if we can take k := heatConductivity*rho as constant then we can handle heat conduction with the constant-coefficient framework with no cross-derivative terms (so merely one-sided boundary fluxes). It was the heat equation, after all, that inspired the one-sided boundary fluxes in the first place. + 10-moment diffusive term Recall the diffusive 10-moment system: -Energy,t = Div Q. where Q = a_1*3*Sym(Grad TI) + a_0*3*Sym(id*tr(3*Sym(Grad TI))), where we now have adopted shortened names: Q = HeatFlux TI = inverse of temperature tensor, id = identity tensor, tr = trace. In components, Q_ijk = Q1_ijk + Q0_ijk, where Q1_ijk = a_1*(TI_jk,i + TI_ij,k + TI_ik,j) and Q0_ijk = a_0*( delta_ij*(2*TI_nk,n + TI_nn,k) + delta_ik*(2*TI_nj,n + TI_nn,j) + delta_jk*(2*TI_ni,n + TI_nn,i)) -Energy_ij,t = Q_ijk,k. In two dimensions -Energy_ij,t = Q_xij,x + Q_yij,y and Q1_xxx/a_1 = 3*TI_xx,x Q1_xxy/a_1 = 2*TI_xy,x + TI_xx,y Q1_xxz/a_1 = 2*TI_xz,x Q1_xyy/a_1 = TI_yy,x + 2*TI_xy,y Q1_xyz/a_1 = TI_yz,x + TI_xz,y Q1_xzz/a_1 = TI_zz,x Q1_yxx/a_1 = 2*TI_xy,x + TI_xy,y Q1_yxy/a_1 = TI_yy,x + 2*TI_xy,y Q1_yxz/a_1 = TI_yz,x + TI_xz,y Q1_yyy/a_1 = 3*TI_yy,y Q1_yyz/a_1 = 2*TI_yz,y Q1_yzz/a_1 = TI_zz,y and Q0_xxx/a_0 = Dx (9*TI_xx+3*TI_yy+3*TI_zz) + Dy (6*TI_yx) Q0_xxy/a_0 = Dx (2*TI_xy) + Dy (TI_xx+TI_yy+TI_zz) Q0_xxz/a_0 = Dx (2*TI_xz) + Dy (2*TI_yz) Q0_xyy/a_0 = Dx (3*TI_xx+TI_yy+TI_zz) + Dy (2*TI_yx) Q0_xyz/a_0 = 0 Q0_xzz/a_0 = Dx (3*TI_xx+TI_yy+TI_zz) + Dy (2*TI_yx) Q0_yxx/a_0 = Dx (2*TI_xy) + Dy (TI_xx+3*TI_yy+TI_zz) Q0_yxy/a_0 = Dx (3*TI_xx+TI_yy+TI_zz) + Dy (2*TI_yx) Q0_yxz/a_0 = 0 Q0_yyy/a_0 = Dx (6*TI_xy) + Dy (3*TI_xx+9*TI_yy+3*TI_zz) Q0_yyz/a_0 = Dx (2*TI_xz) + Dy (2*TI_yz) Q0_yzz/a_0 = Dx (2*TI_xy) + Dy (TI_xx+3*TI_yy+TI_zz), where Dx := partial derivative in the x direction and Dy := partial derivative in the y direction. Note that there is a bit of redundancy here, since Q0_xxz = Q0_yyz, Q0_xyy = Q0_xzz, Q0_yxx = Q0_yzz, Q0_xyz = Q0_yxz = 0. The pattern here is that if there is a repeated index then it is the nonrepeated index that matters, and if there is no repeated index then the term is zero.