Conservative scheme compatible with some other conservation laws: conservation of the local angular momentum

1 Introduction

We are interested in the solution of the Euler equations for compressible flows

\frac{\partial u}{\partial t} + div f (u) = 0,

(1)

defined on a space-time domain $Ω \times T$ , where

u = (ρ, ρ v, E) and f = (ρ v, ρ v \otimes v + p Id, (E + p) v) = (f_{1}, \dots, f_{d}), d = 1, 2, 3.

The total energy writes

E = ρ e + \frac{1}{2} ρ v^{2},

and the pressure is given by the equation of state $p = p (ρ, e)$ . This function satisfies the usual convexity requirements. Equation (1) is complemented by initial and boundary conditions.

In addition to (1), the solution $u$ is known to satisfy other conservation relation. In the case of smooth flows, it is known that $u$ satisfies:

\frac{\partial S}{\partial t} + div (ρ S) = 0 (\leq 0 in % general .)

(2)

where $S = - ρ s$ and $s$ is the specific entropy defined by the Gibbs relation:

T d s = d (\frac{e}{ρ}) - \frac{p}{ρ^{2}} d ρ .

In order to fully define the entropy, we need a complete equation of state, $p = p (s, ρ)$ from which one can deduce $p = p (ρ, s)$ from standard thermodynamical relations.

There are indeed other relations, for example the preservation of the angular momentum. Defining $J = ρ x \land v$ , and from

\frac{\partial ρ v}{\partial t} + div (ρ v \otimes v) + \nabla p = 0,

we write, in 2D,

\begin{matrix} \frac{\partial ρ v}{\partial t} + \frac{\partial f_{1}}{\partial x} + \frac{\partial f_{2}}{\partial y} = 0, f_{1} = (\begin{matrix} ρ v_{x}^{2} + p ρ v_{x} v_{y} \end{matrix}) = (\begin{matrix} f_{1}^{1} f_{1}^{2} \end{matrix}), f_{2} = (\begin{matrix} ρ v_{x} v_{y} ρ v_{y}^{2} + p \end{matrix}) = (\begin{matrix} f_{2}^{1} f_{2}^{2} \end{matrix}), \end{matrix}

We note that $f_{1}^{2} = f_{2}^{1}$ .

\begin{matrix} \frac{\partial ρ x \land v}{\partial t} & = x \land \frac{\partial ρ v}{\partial t} = - x \land (\frac{\partial f_{1}}{\partial x_{1}} + \frac{\partial f_{2}}{\partial x_{2}}) = - x_{1} (\frac{\partial f_{1}^{2}}{\partial x_{1}} + \frac{\partial f_{2}^{2}}{\partial x_{2}}) + x_{2} (\frac{\partial f_{1}^{1}}{\partial x_{1}} + \frac{\partial f_{2}^{1}}{\partial x_{2}}) = - \frac{\partial}{\partial x_{1}} (f_{1}^{2} x_{1} - x_{2} f_{1}^{1}) - \frac{\partial}{\partial x_{2}} (x_{1} f_{2}^{2} - f_{2}^{1} x_{2}) + (f_{1}^{2} - f_{2}^{1}), \end{matrix}

Hence we have:

\frac{\partial J}{\partial t} + \frac{\partial}{\partial x_{1}} (f_{1}^{2} x_{1} - x_{2} f_{1}^{1}) + \frac{\partial}{\partial x_{2}} (x_{1} f_{2}^{2} - f_{2}^{1} x_{2}) = 0.

(3)

In 3 dimensions, we have a similar results: a direct computation of the time evolution of the angular momentum equation give:

and since

f_{2}^{3} = f_{3}^{2} = ρ u_{2} u_{3}, f_{3}^{1} = f_{1}^{3} = ρ u_{3} u_{1}, f_{1}^{2} = f_{2}^{1} = ρ u_{1} u_{2},

we obtain the additional conservation relation:

\begin{matrix} \frac{\partial J_{1}}{\partial t} & + \frac{\partial}{\partial x_{1}} (x_{2} f_{1}^{3} - x_{3} f_{1}^{2}) + \frac{\partial}{\partial x_{2}} (x_{2} f_{2}^{3} - x_{3} f_{2}^{2}) + \frac{\partial}{\partial x_{3}} (x_{2} f_{3}^{3} - x_{3} f_{3}^{2}) = 0, \frac{\partial J_{2}}{\partial t} & + \frac{\partial}{} \partial x_{1} (x_{3} f_{1}^{1} - x_{1} f_{1}^{3}) + \frac{\partial}{\partial x_{2}} (x_{3} f_{2}^{1} - x_{1} f_{2}^{3}) + \frac{\partial}{\partial x_{3}} (x_{3} f_{3}^{1} - x_{1} f_{3}^{3}) = 0, \frac{\partial J_{3}}{\partial t} & + \frac{\partial}{} \partial x_{1} (x_{1} f_{1}^{2} - x_{2} f_{1}^{1}) + \frac{\partial}{\partial x_{2}} (x_{1} f_{2}^{2} - x_{2} f_{2}^{1}) + \frac{\partial}{\partial x_{3}} (x_{1} f_{3}^{2} - x_{2} f_{3}^{1}) = 0. \end{matrix}

(4)

The weak form of the relations (3) and (4) are also satisfied by $J$ . To see this, in the 2D case, one has to take the test functions $x_{1} φ$ and $x_{2} φ$ where $φ$ has compact support and is $C^{1}$ in space-time domain, and subtract to get the weak form of (3). The same method also works in 3D.

Remark 1.1.

The PDEs (3) and (4) are translational invariant. This means that if $J$ satisfies (3) or (4) in the frame of coordinates $x$ , then $J^{'} = J + ρ a \land v$ will satisfy (3) or (4) in the frame of coordinates $x^{'} = x + a$ .

Developing schemes that locally preserves the angular momentum has recently attracted some attention, see e.g [9, 10] and the reference therein. In [5, 7, 1, 4] has been developed a simple technique that enables to correct an initial scheme so that an additional conservation is satisfied. This has been applied to some non conservative form of the Euler equation, also to satisfy entropy conservation when this is relevant, and to various schemes: residual distribution schemes and discontinuous Galerkin ones. The goal of this paper is to show how one can extend the method to the problem of angular preservation, for second and higher order schemes, in a fully discrete manner, both in time and space.

The paper is organised as follows. In section 2 we briefly introduce a high order residual distribution scheme for the unsteady version of (1). In section 3 we construct a second order scheme for (1) that preserves locally the angular momentum; then we show this is possible for triangles and tetrahedrons. In section 4 we provide a possible solution to conservation of momentum for higher than second order schemes. In section 5 we study the validity and effectiveness of the proposed strategy, considering three different problems. In section 6 we describe how to extend this to discontinuous Galerkin schemes. Finally, section 7 provide some conclusions and future perspectives.

2 Numerical scheme

2.1 Some generalities

The goal is to construct a scheme for (1) that preserves locally the angular momentum. For this we use a high order residual distribution scheme for the unsteady version of (1) that we briefly describe now, see [6, 2] for more details.

We are given a tessellation of $Ω$ by non overlapping simplex denoted generically by $K$ , and more precisely $Ω = \cup K$ . We assume the mesh to be conformal which means that $K \cap K^{'}$ is either empty, reduce to a full face or a full edge or a vertex. From this, we consider $V = {f : Ω \to R, f_{| K} \in P^{r} \forall K \in Ω} \cap C^{0} (Ω)$ . We are given a basis of $V$ , ${φ_{σ}}$ where the $σ$

are the degrees of freedom (DOFs). The restriction of

φ_{σ}

to any

K

is assumed here to be a Bézier polynomial, see appendix C for the notations we use. The reason is that for any

σ

, we have

| C_{σ} | = \int_{Ω} φ_{σ} d x > 0.

(5)

This is not necessarily true for Lagrange polynomials.

We are looking for an approximation of the solution of the form

u = \sum σ u_{σ} φ_{σ},

Then the solution $u^{n}$ at time $t_{n}$ can be written as follows

u^{n} = \sum σ u_{σ}^{n} φ_{σ} .

The coefficients $u_{σ}^{n}$ , are chosen by a numerical method and we use the following technique to calculate the coefficients in the approximation.

2.2 Residual distribution scheme for steady problems

First, we consider a steady version of system (1)

div f (u) = 0.

(6)

The main steps of the residual distribution approach can be summarized as follows

For any element $K \in Ω$ , compute a fluctuation term (total residual) (see figure 1(a))

$Φ^{K} (u) = \int_{\partial K} f (u) \cdot n d x (= \int_{K} div f (u) d x),$ (7)
For every DOF $σ$ within the element K, define the nodal residuals $Φ_{σ}^{K}$ as the contribution to the fluctuation term $Φ^{K}$ (see figure 1(b)) such that

$Φ^{K} (u) = \sum σ \in K Φ_{σ}^{K}, \forall K \in Ω,$ (8)
The resulting scheme is obtained for each DOF $σ$ by collecting all the nodal residual contributions $Φ_{σ}^{K}$ from all elements K surrounding a node $σ \in Ω$ (see figure 1(c))

$\sum K, σ \in K Φ_{σ}^{K} (u) = 0, \forall σ \in Ω .$ (9)

We can write a similar formulation for the boundary conditions. For any DOF $σ \in Ω$ we can split (6) into the internal and boundary contributions

\sum K, σ \in K Φ_{σ, x}^{K} (u) + \sum Γ, σ \in γ Φ_{σ, x}^{Γ} (u) = 0,

(10)

where $Φ_{σ, x}^{K}$ and $Φ_{σ, x}^{γ}$ are the residuals corresponding to the spatial discretization and $γ$ is an edge on the boundary $Γ$ of $Ω$ . If we suppose the Dirichlet boundary condition $u = g$ on $Γ$ , for any K and $Γ$ , $Φ_{σ, x}^{K}$ and $Φ_{σ, x}^{γ}$ satisfy the following conservation relations

	$\sum K, σ \in K Φ_{σ, x}^{K} (u) = \int_{\partial K} f (u) . n,$	(11a)
and for the boundary condition (here Dirichlet in weak form)
	$\sum Γ, σ \in Γ Φ_{σ, x}^{Γ} (u) = \int_{Γ} (F_{n} (u, g) - f (u) . n) .$	(11b)

In the appendix A, we provide several examples of such fluctuations. For now, and in order to simplify the discussion on the approximation of the unsteady problem (1), we introduce a ”variational” form of the fluctuation. In all the known cases, we can write $Φ_{σ}^{K}$ as

(12)

where

\sum σ \in K β_{σ}^{K} (u) = Id (= 1 % i n t h e s c a l a r c a s e) .

The distribution coefficients can be scalar in the scalar case, and matrices in the system case.

Then, we can rewrite the scheme in a Petrov-Galerkin fashion:

where $ξ_{σ}^{K} (u) = β_{σ}^{K} (u) - φ_{σ} Id$ . To avoid confusion, we will write $ξ_{σ}^{K} (u)$ as $ξ_{σ}^{K}$ , thus removing the functional dependency of this term with respect to the solution $u$ . We remove it, but we do not forget it!

We note that the functions $ξ_{σ}^{K}$ satisfies

\sum σ \in K ξ_{σ}^{K} = 0.

Then we define $Ψ_{σ, u} = φ_{σ} Id + ξ_{σ}$ with, for $x \in K$ ,

ξ_{σ} (x) = ξ_{σ}^{K} .

Defining $W = span (Ψ_{σ})$ , we can formally rewrite the scheme as: for any $v \in V$ , find $u$ such that

- \int_{Ω} \nabla v \cdot f (u) d x + \int_{\partial Ω} v F_{n} d γ + \sum K \int_{K} ξ (u, v) div f (u) d x = 0,

where

ξ (x) = \sum σ v_{σ} ξ_{σ} .

Here the functional dependence of $ξ$ with respect to $u$ exists but is implicit.

2.3 Residual distribution scheme for unsteady problems

Using the results presented above, it becomes possible to describe the unsteady version of the scheme. In order to get a consistent approximation, we simply ”multiply” (1) with the Petrov-Galerkin test function $v + ξ (v)$ , and integrate in time. Note that $v$ will not depend on time. Between $t_{n}$ and $s \in] t_{n}, t_{n + 1}]$ , we get

\int_{t_{n}}^{s} \int_{Ω} (v + ξ (v)) \frac{\partial u}{\partial t} d x d t - \int_{t_{n}}^{s} \int_{Ω} \nabla v \cdot f (u) d x d t + \int_{t_{n}}^{s} \int_{\partial Ω} v F_{n} d γ d t + \sum K \int_{t_{n}}^{s} \int_{K} ξ (v) div f (u) d x d t = 0.

Since $v$ is independent of time, this can be equivalently rewritten as:

\int_{Ω} (v + ξ (v)) (u (x, s) - u (x, t_{n})) d x - \int_{t_{n}}^{s} \int_{Ω} \nabla v \cdot f (u) d x d t + \int_{t_{n}}^{s} \int_{\partial Ω} v F_{n} d γ d t + \sum K \int_{t_{n}}^{s} \int_{K} ξ (v) div f (u) d x d t = 0,

and then, for any $σ$ ,

This suggest to introduce the space-time fluctuations

Φ_{σ}^{K} (u, s) = \int_{K} (φ_{σ} + ξ_{σ}) (u (x, s) - u (x, t_{n})) d x      Φ_{σ, t}^{K} (u, s) + \int_{t_{n}}^{s} Φ_{σ, x}^{K} (u, s) d t

(13)

where

Φ_{σ, x}^{K} (u, s) = - \int_{K} \nabla φ_{σ} \cdot f (u) d x + \int_{K} ξ_{σ} div f (u) d x + \int_{\partial K} φ_{σ} f (u) \cdot n d γ .

2.4 Iterative timestepping method

In order to integrate the system, we proceed in two steps. First we introduce sub-time steps in the interval $[t_{n}, t_{n + 1}]$ : we subdivide the interval $[t_{n}, t_{n + 1}]$ into sub-intervals obtained from a partition

t_{n} = t_{(0)} < t_{(1)} < \dots < t_{(l)} < \dots < t_{(M)} = t_{n + 1} .

In all our examples, this partition will be regular: $t_{(l)} = t_{n} + \frac{l}{M} Δ t$ , but other choices could have been made, especially for very high order accuracy in time. We construct an approximation of $u$ at times $t_{(l)}$ , denoted by $u_{(l)} \approx u (t_{(l)})$ . Then we introduce $I_{M}$

the Lagrange interpolant (in time) of degree

M

defined from this subdivision. We also introduce the piece-wise constant interpolant

I_{0}

defined as

I_{0} (s) = u_{(0)} if s \in [t_{n}, t_{n + 1}] .

Another choice could have been made

I_{0} (s) = u_{(l)} if s \in [t_{(l)}, t_{(l + 1)} [.

for $l = 0, \dots, M - 1$ . The notation $U$

represents the vector

U = (u_{(0)}, u_{(1)}, \dots, u_{(M)})

i.e the vector of all the approximations for the sub-steps. Note that

u_{(0)} = u^{n}

and

u_{(M}) = u^{n + 1}

. We need residuals,

{Φ_{σ, t}^{K, p} (U), Φ_{σ, x}^{K, p} (U)}

computed for time

t_{p}

and that satisfies, for any

p

We set $Φ_{σ}^{K} (U) = Φ_{σ, t}^{K, p} (U) + Φ_{σ, x}^{K, p} (U)$ .

Then we introduce two approximations of (1):

A first order approximation in time: for any $σ \in K$ and $l = 1, \dots, M$ ,

$[L_{1} (U; u^{n})]_{σ, (l)} := | C_{σ} | u_{σ, (l)} - | C_{σ} | u_{σ, (0)} + \sum K, σ \in K \int_{t_{(0)}}^{t_{(l)}} I_{0} (Φ_{σ, x}^{K} (U), s) d s,$ (14)

where

$I_{0} (Φ_{σ, x}^{K} (U), s) = I_{0} (Φ_{σ, x}^{K} (u_{(0)}), Φ_{σ, x}^{K} (u_{(1)}), \dots, Φ_{σ, x}^{K} (u_{(M)}), s),$

i.e. here

$[L_{1} (U; u^{n})]_{σ, (l)}$ $= | C_{σ} | u_{σ, (l)} - | C_{σ} | u_{σ, (0)} + (t_{(l)} - t_{(0)}) \sum K, σ \in K I_{0} (Φ_{σ, x}^{K} (u_{(0)}))$

This is a first order explicit approximation in time.
A high order approximation: for any $σ \in K$ and $l = 1, \dots, M$ ,

$[L_{2} (U; u^{n})]_{σ, (l)} := \sum K, σ \in K \int_{t_{(0)}}^{t_{(l)}} I_{M} (Φ_{σ}^{K} (U), s) d s,$ (15)

After having performed exact integration in time to obtain the approximation for every sub-steps, the time integration can be written in the form

$\int_{t_{(0)}}^{t_{(l)}} I_{M} (Φ_{σ, x}^{K} (U)) d s = M \sum k = 0 θ_{k}^{l} Φ_{σ, x}^{K} (U) .$ (16)

Ideally, we would like to solve for each $σ$ and each $l$ ,

[L_{2} (U; u^{n})]_{σ, (l)} = 0,

(17)

but this is very difficult in general and the resulting scheme derived by $L_{2}$ operator is implicit. Instead we use a defect correction method and proceed within the time interval $[t_{n}, t_{n + 1}]$ as follows:

Set $U^{(0)} = (u^{n}, u^{n}, \dots, u^{n})$ ,
For any $p \geq 0$ , define $U^{(p + 1)}$ by

(18)

Since $L_{1}$ is explicit, $U^{(p + 1)}$ can be obtained explicitly. In our case, this amounts to a multi-step method where each step writes as

\begin{matrix} | C_{σ} | (U_{σ}^{(p + 1)} - U_{σ}^{(p)}) & = - \sum K, σ \in K Φ_{σ}^{K} (U^{(p)}), \end{matrix}

(19)

The scheme (18) is completely explicit. One has the following result, see [6]:

Proposition 2.1.

If two operators $L_{1}$ and $L_{2}$ depending on the discretization scale $Δ = Δ t$ , are such that:

There exists a unique $U_{Δ}^{⋆}$ such that $L_{2} (U_{Δ}^{⋆}) = 0.$
$L_{1}$ is coercive, i.e., there exists $α_{1} > 0$ independent of $Δ$ , such that for any U and V,

$α_{1} ∥ U - V ∥ \leq ∥ L_{1} (U) - L_{1} (V) ∥,$
$L_{1} - L_{2}$ is uniformly Lipschitz continuous with Lipschitz constant $α_{2} Δ$ , i.e., there exists $α_{2} > 0$ independent of $Δ$ , such that for any U and V,

$∥ (L_{1} (U) - L_{2} (U)) - (L_{1} (V) - L_{2} (V)) ∥ \leq α_{2} Δ ∥ U - V ∥ .$

Then if $ν = \frac{α_{2} Δ}{α_{1}} < 1$ the defect correction method is convergent, and after p iterations the error is smaller than $ν^{p} ∥ U^{(0)} - U_{Δ}^{⋆} ∥ .$

Therefore if there is a unique solution $U_{Δ}^{⋆}$ of (17), then after $k$ steps, $u_{(M)}^{(k)} = U_{Δ}^{⋆} + O (Δ t^{k})$ , so that only $M + 1$ steps are needed. Of course this is true thanks to the conditions (5), (14), (15) and the condition that $L_{1}$ is a first order approximation of $L_{2}$ :

L_{1} - L_{2} = O (Δ t) .

(20)

3 Angular momentum preservation: second order case

It is clear that the residuals depends where the order in time. We can compute the variation of the angular momentum from (19). For any $σ$ , we have

| C_{σ} | (x_{σ} \land (m_{σ}^{(p + 1)} - m_{σ}^{(p)})) + \sum K, σ \in K x_{σ} \land Φ_{m, σ}^{K} (U^{(p)}) + \sum Γ, σ \in Γ x_{σ} \land Φ_{m, σ}^{Γ} (U^{(p)}) = 0,

so that

where here $Φ_{m, σ}^{K}$ is the momentum component of the residual at the DOF $σ$ in the element $K$ and $Φ_{m, σ}^{Γ}$ is the momentum component of the residual at the DOF $σ$ in the boundary $Γ$ . An easy condition for having local conservation is:

\sum σ \in K x_{σ} \land Φ_{m, σ}^{K} (U^{(p)}) = \int_{K} (J^{(p)} - J^{(0)}) d x + Δ t \int_{\partial K} I_{2} (G (U^{(p)})) \cdot n d γ := Φ_{J}^{K}

(21)

\sum σ \in Γ x_{σ} \land Φ_{m, σ}^{Γ} (U^{(p)}) = Δ t \int_{\partial K} (^G - I_{2} (G (U^{(p)})) \cdot n) d γ := Φ_{J}^{Γ}

(22)

with $G$ the flux for the angular momentum and

G (U^{(p)}) = \frac{1}{2} (G (U^{(p)}) + G (U^{(0)}))

depending on how the system (1) has been discretized.

Of course, in general, the relation (21) and (22) cannot be satisfied. In order to be satisfied, as well as keeping

(23)

\sum σ \in Γ Φ_{m, σ}^{Γ} = Δ t \int_{Γ} F_{m} (U^{(p)}) \cdot n d γ

(24)

with a similar definition of the momentum flux, we will present a perturbation of the momentum residuals. At this level, it is important to provide the quadrature formula. For (21) and (23), we consider

\int_{K} f (x) d x \approx \frac{| K |}{3} \sum σ \in K f_{σ},

so it is exact for (23) and only approximate for (21), but second order.

To achieve this, following [7], we introduce a perturbation of the momentum residual, $r_{σ}^{K}$ , such that the new momentum residual is $Φ_{m, σ}^{K} + r_{σ}^{K}$ . We must have:

\begin{matrix} \sum σ \in K & r_{σ}^{K} = 0 \sum σ \in K & x_{σ} \land r_{σ}^{K} = Φ_{J}^{K} - \sum σ \in K x_{σ} \land Φ_{m, σ}^{K} (U^{(p)}) := Ψ^{K} \end{matrix}

(25)

Since (22) is not true, as above, we introduce a vectorial correction of the momentum residual $r_{σ}^{Γ}$ . We need:

\begin{matrix} \sum σ \in Γ & r_{σ}^{Γ} = 0 \sum σ \in Γ & x_{σ} \land r_{σ}^{Γ} = Φ_{J}^{Γ} - \sum σ \in Γ x_{σ} \land Φ_{m, σ}^{Γ} (U^{(p)}) := Ψ^{Γ} \end{matrix}

(26)

$r_{σ}^{Γ}$ can be derived in the same way.

In dimension $d$ and for an element with $p$ DOFs, we have $2 \times d$ equations for $2 \times p$ unknowns. Since $p$ is always larger than $d + 1$ , the system has a priori solutions. However, finding a close form formula is element dependent. In the following, we show this is possible for triangles and tetrahedrons. Once these solutions are described, we have a scheme that locally preserves the angular momentum $J_{σ} = x_{σ} \land m_{σ}$ and globally conserves

up to boundary terms.

Remark 3.1 (Translation invariance).

The angular momentum is defined after a frame has been defined, and if one makes a translation of vector $a$ , the scheme is translational invariant as in the continuous case.

3.1 Solution for triangular elements

We first have $r_{1} = - r_{2} - r_{3}$ , so

(x_{2} - x_{1}) \land r_{2} + (x_{3} - x_{1}) \land r_{3} = Ψ \in R

We define

r_{2} = r (x_{3} - x_{1}), r_{3} = r (x_{1} - x_{2}), r \in R

and get

r (det (x_{2} - x_{1}, x_{3} - x_{1}) + det (x_{3} - x_{1}, x_{1} - x_{2})) = Ψ,

i.e. since $det (x_{2} - x_{1}, x_{3} - x_{1}) = 2 | T |$ , we have:

r = \frac{Ψ}{4 | T |}, r_{2} = r (x_{3} - x_{1}), r_{3} = r (x_{1} - x_{2}), r_{1} = r (x_{2} - x_{3}) .

(27)

3.2 Solution for tetrahedrons

Again, $r_{1} = - \sum_{j = 2}^{4} r_{j}$ , so that

4 \sum j = 2 (x_{j} - x_{1}) \land r_{j} = Ψ .

In order to simplify the notations, we introduce $x_{i j} = x_{i} - x_{j}$ . We are looking for

\begin{matrix} r_{2} & = α_{2} x_{13} + β_{2} x_{14} r_{3} & = α_{3} x_{12} + β_{3} x_{14} r_{4} & = α_{4} x_{12} + β_{4} x_{13} \end{matrix}

This gives:

\begin{matrix} Ψ & = α_{2} x_{21} \land x_{13} + β_{2} x_{21} \land x_{14} + α_{3} x_{31} \land x_{12} + β_{3} x_{31} \land x_{14} + α_{4} x_{41} \land x_{12} + β_{4} x_{41} \land x_{13} \end{matrix}

Then, since $a \cdot (b \land c) = det (a, b, c)$ ,

\begin{matrix} Ψ \cdot x_{21} & = β_{3} det (x_{21}, x_{31}, x_{14}) + β_{4} det (x_{21}, x_{41}, x_{13}) Ψ \cdot x_{31} & = β_{2} det (x_{31}, x_{21}, x_{14}) + α_{4} det (x_{31}, x_{41}, x_{12}) Ψ \cdot x_{41} & = α_{2} det (x_{41}, x_{21}, x_{13}) + α_{3} det (x_{41}, x_{31}, x_{12}) \end{matrix}

Now,

\begin{matrix} 3 | T | & = det (x_{21}, x_{31}, x_{14}) = - det (x_{21}, x_{41}, x_{13}) = - det (x_{31}, x_{21}, x_{14}) = det (x_{31}, x_{41}, x_{12}) = det (x_{41}, x_{21}, x_{13}) = - det (x_{41}, x_{31}, x_{12}) \end{matrix}

so that we have:

\begin{matrix} Ψ \cdot x_{21} & = 3 | T | (β_{3} - β_{4}) Ψ \cdot x_{31} & = 3 | T | (- β_{2} + α_{4}) Ψ \cdot x_{41} & = 3 | T | (α_{2} - α_{3}) . \end{matrix}

This suggests to assume $β_{3} = - β_{4}$ , $α_{4} = - β_{2}$ , $α_{2} = - α_{3}$ , i.e

\begin{matrix} β_{3} = - β_{4} = \frac{Ψ \cdot x_{21}}{6 | T |} α_{4} = - β_{2} = \frac{Ψ \cdot x_{31}}{6 | T |} α_{2} = - α_{3} = \frac{Ψ \cdot x_{41}}{6 | T |} \end{matrix}

Remembering that $u \land (v \land w) = (u \cdot w) v - (u \cdot v) w$ , we see that

\begin{matrix} r_{2} & = \frac{1}{6 | T |} (Ψ \cdot x_{41} x_{13} - Ψ \cdot x_{31} x_{14}) = \frac{1}{6 | T |} Ψ \land (x_{14} \land x_{31}) r_{3} & = \frac{1}{6 | T |} (Ψ \cdot x_{21} x_{14} - Ψ \cdot x_{41} x_{12}) = \frac{1}{6 | T |} Ψ \land (x_{21} \land x_{14}) r_{4} & = \frac{1}{6 | T |} (Ψ \cdot x_{31} x_{12} - Ψ \cdot x_{21} x_{13}) = \frac{1}{6 | T |} Ψ \land (x_{31} \land x_{12}) . \end{matrix}

Last,

\begin{matrix} x_{14} \land x_{31} + x_{21} \land x_{14} + x_{31} \land x_{12} & = x_{14} \land x_{31} + x_{14} \land x_{12} + x_{31} \land x_{12} = x_{14} \land x_{31} + (x_{14} + x_{31}) \land x_{12} = x_{14} \land x_{31} + x_{34} \land x_{12} = (x_{13} + x_{34}) \land x_{31} + x_{34} \land x_{12} = x_{34} \land x_{31} + x_{34} \land x_{12} = x_{34} \land x_{32} \end{matrix}

and then we have:

(28)

4 Angular momentum preservation: the high order case

The idea is similar, the key point is to characterize the quadrature formula that describes $\int_{K} J d x$ . The additional difficulty is that $f_{σ} = f (σ) + O (h^{2})$ if we make the geometrical identification of the DOFs with the Greville points. Here to simplify the notations, $Φ_{σ}^{K}$ denotes the residual at $σ \in K$ evaluated for the momentum, it does not contain the contribution for the density or the energy. In this section we always use this short hand notation, except at the end.

We start again from (19) that we rewrite as

U_{σ}^{(p + 1)} = U_{σ}^{(p)} + δ U_{σ}^{(p)} .

Then

U^{(p + 1)} = \sum σ U_{σ}^{(p + 1)} φ_{σ} .

The angular momentum, integrated, is

We introduce the notations

y_{σ} = \frac{1}{| C_{σ} |} \int_{Ω} x B_{σ} d x,

z_{σ}^{K} = \frac{1}{| K |} \int_{K} x B_{σ} d x,

so that we can rewrite the total kinetic momentum as:

\int_{Ω} x \land m^{(p + 1)} d x = \sum σ | C_{σ} | y_{σ} \land m_{σ}^{(p + 1)},

and then we can write the update:

\begin{matrix} \sum σ | C_{σ} | y_{σ} \land m_{σ}^{(p + 1)} & = \sum σ | C_{σ} | y_{σ} \land m_{σ}^{(p)} - \sum σ y_{σ} \land (\sum K, σ \in K Φ_{m, σ}^{K}) = \sum σ | C_{σ} | y_{σ} \land m_{σ}^{(p)} - \sum σ y_{σ} \land (\sum K, σ \in K Φ_{m, σ}^{K} + \sum Γ, σ \in Γ Φ_{m, σ}^{Γ}) = \sum σ | C_{σ} | y_{σ} \land m_{σ}^{(p)} - \sum \end{matrix}

Conservative scheme compatible with some other conservation laws: conservation of the local angular momentum

Authors

Invariant finite-difference schemes with conservation laws preservation for one-dimensional MHD equations

Memory erasure with finite-sized spin reservoir

A Harten's Multiresolution Framework for Subdivision Schemes

The rotating rigid body model based on a non-twisting frame

Lagrangian Neural Network with Differentiable Symmetries and Relational Inductive Bias

Locally conservative and flux consistent iterative methods

A novel formulation for the evolution of relativistic rotating stars

1 Introduction

Remark 1.1.