<jupyter_start><jupyter_text>Finding Numerical Fluxes for DGCopyright (C) 2020 <NAME>MIT LicensePermission is hereby granted, free of charge, to any person obtaining a copyof this software and associated documentation files (the "Software"), to dealin the Software without restriction, including without limitation the rightsto use, copy, modify, merge, publish, distribute, sublicense, and/or sellcopies of the Software, and to permit persons to whom the Software isfurnished to do so, subject to the following conditions:The above copyright notice and this permission notice shall be included inall copies or substantial portions of the Software.THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS ORIMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THEAUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHERLIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS INTHE SOFTWARE.----NB: This note book uses the [Maxima-Jupyter kernel](https://github.com/robert-dodier/maxima-jupyter) to interface with [Maxima](http://maxima.sourceforge.net/), an open-source computer algebra system. I have found that the [Docker image](https://hub.docker.com/r/calyau/maxima-jupyter) is a fairly convenient way of using this kernel.Some Maxima version info for reproducibility:<jupyter_code>build_info();<jupyter_output><empty_output><jupyter_text>Load some packages:<jupyter_code>kill(all);
load("itensor");

assert(condition):=if not condition then error("Assertion violated") else true$

norm_2_squared(v):=v.v;

crossfunc(f):=makelist(
  sum(sum(
      levi_civita([i,j,k])*f(j,k),
   j,1,3),k,1,3),i,1,3)$

crossprod(a,b):=crossfunc(lambda([j,k], a[j]*b[k]));<jupyter_output>REDEFINITION-WITH-DEFMACRO: 
  redefining MAXIMA::IFNOT in DEFMACRO

REDEFINITION-WITH-DEFMACRO: 
  redefining MAXIMA::M+OR*OR^P in DEFMACRO

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::$IDUMMY-IMPL in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::$IDUMMY in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::ISPROD in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::DERAT in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::PLUSI in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::MINUSI in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::COVI in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::CONTI in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::DERI in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::NAME in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::$COVI-IMPL in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::$COVI in DEFUN

REDEFINITION-WITH-DEFUN: 
  redefining MAXIMA::$CONTI-IMPL in DEFUN

REDEFINITION-WITH-D[...]<jupyter_text>Simplification UtilitiesThese function simplify expressions coming from the (symbolic) simultaneous diagonalization, by letting Maxima realize that $n \cdot n =1$.<jupyter_code>/* ------------------------------------------------------------------------- */
/* Simplification for expressions stemming from hyperbolic systems */
/* ------------------------------------------------------------------------- */

hypsimp(x):=ratsimp(ratsubst(1,n.n,x))$

fullhypsimp(x):=hypsimp(
  ratsubst(
    last(n)^2,
    1-sum(n[i]^2,i,1,length(n)-1),
    x)
  )$
  
/* ------------------------------------------------------------------------- */
/* diagonalize a given hyperbolic operator A */
/* ------------------------------------------------------------------------- */

hypdiagonalize(A):=block([evA, V, invV,D],
  evA:hypsimp(apply(append, eigenvectors(A)[2])),
  V:transpose(apply(matrix, evA)),
  invV:hypsimp(invert(V)),
  assert(hypsimp(V.invV)=ident(length(A))),
  D:hypsimp(invV.A.V),
  [V, D, invV])$<jupyter_output><empty_output><jupyter_text>Wave Equation<jupyter_code>/* redefine this to change dimensionality: */
n:[nx,ny]; 
dims:length(n);

assume(c>0);

if dims = 1 then n:[1];<jupyter_output><empty_output><jupyter_text>Define the flux jacobian for the first-order form of the wave-equation:$$\begin{align*}\partial_t u &= \nabla \cdot \boldsymbol v\\\partial_t \boldsymbol v &= \nabla u\end{align*}$$projected onto a line with normal `n`:<jupyter_code>esymmatrix(n, v, i,j):=ematrix(n,n,v,i,j)+ematrix(n,n,v,j,i);
wave_A:sum(n[i]*esymmatrix(dims+1, -c, 1+i,1),i,1,dims);<jupyter_output><empty_output><jupyter_text>Find the eigenvalues of the flux Jacobian:<jupyter_code>[wave_V, wave_D, wave_invV]:hypdiagonalize(wave_A);<jupyter_output><empty_output><jupyter_text>Finding the Numerical FluxThis function sets up a system of Rankine-Hugoniot conditions across the flux fan and solves for the stagnation state:<jupyter_code>/* ------------------------------------------------------------------------- */
/* compute upwind flux for a given operator with eigenvalues evs, sorted
 * in ascending order.
 * Sign assumptions for all variables occuring in evs must be in place.
 */
/* ------------------------------------------------------------------------- */
hyp_upwind_flux(evs, D):=block([evvars, Dp, Dm, n, midstates, states, unknowns],
  evvars:listofvars(evs),

  add_evvars_suffix(suffix, x):=subst(makelist(v=concat(''v, suffix), v, evvars), x),

  evsm:add_evvars_suffix(m, evs),
  evsp:add_evvars_suffix(p, evs),

  Dm:add_evvars_suffix(m, D),
  Dp:add_evvars_suffix(p, D),

  midstates:makelist(makelist(concat(s,state,i), i, 1, length(D)),
      state, 1, length(evs)-1),

  states:append(
    [makelist(concat(sm, i), i, 1, length(D))],
    midstates,
    [makelist(concat(sp,i), i, 1, length(D))]),

  unknowns:apply(append, midstates),

  result:if member(0, evs) then
    block([biasedD, veceqns, eqns, soln],
      biasedD:makelist(
        if evs[i] = 0 then [Dp,Dm]
        else if evs[i] > 0 then [Dp,Dp]
        else [Dm,Dm],
        i, 1, length(evs)),

      veceqns:apply(append, makelist(
        -(if evs[i] > 0 then evsp[i] else evsm[i]) *(states[i+1]-states[i])
        +(biasedD[i][1].states[i+1]-biasedD[i][2].states[i]),
        i,1,length(evs))),

      eqns:makelist(veceqns[i,1], i, 1, length(veceqns)),

      soln:solve(eqns, unknowns),
      assert(length(soln)=1),

      for i: 1 thru length(evs) do
        if evs[i] = 0 then return(Dp.subst(soln[1], midstates[i]))
    )
  else
    block([straddle_idx, Dstates, veceqns, eqns, soln],
      straddle_idx:for i: 1 thru length(evs)-1 do
        if (evs[i] < 0) and (evs[i+1] > 0) then return(i),

      flux:makelist(concat(flux,i),i,1,length(D)),

      unknowns:append(unknowns, flux),

      Dstates:append(
        [Dm.first(states)],
        makelist(
          if i = straddle_idx then flux
          else if evs[i] > 0 then Dp.midstates[i]
          else Dm.midstates[i],
          i, 1, length(midstates)),
        [Dp.last(states)]),

      veceqns:apply(append, makelist(
        -(if evs[i] > 0 then evsp[i] else evsm[i]) *(states[i+1]-states[i])
        +(Dstates[i+1]-Dstates[i]),
        i,1,length(evs))),

      eqns:makelist(veceqns[i,1], i, 1, length(veceqns)),

      print(covect(eqns)),
      soln:solve(eqns, unknowns),
      assert(length(soln)=1),

      subst(soln[1], flux)
    ),
  subst(
    append(
      makelist(concat(sm, i)=sm[i,1], i, 1, length(D)),
      makelist(concat(sp, i)=sp[i,1], i, 1, length(D))
      ),
    result)
  )$<jupyter_output><empty_output><jupyter_text>Find an expression for the flux in characteristic variables.Note the `p` and `m` suffixes for the $+$ and $-$ sides of the interface.<jupyter_code>wave_eigenvalues:makelist(wave_D[i,i], i, 1, length(wave_D));

if member(0, wave_eigenvalues) then
  wave_sflux:hyp_upwind_flux([-c,0,c], wave_D)
else
  wave_sflux:hyp_upwind_flux([-c,c], wave_D);<jupyter_output>solve: dependent equations eliminated: (1 6 8)<jupyter_text>Convert back to conserved variables:<jupyter_code>wave_wflux:ratsimp(wave_V.ev(wave_sflux, [sm=wave_sminw,sp=wave_spinw]));<jupyter_output><empty_output><jupyter_text>Maxwell's EquationsFirst, set up some parameter assumptions:<jupyter_code>assume(c>0);
assume(mu>0);
assume(epsilon>0);
assume(epsinv>0);
assume(muinv>0);<jupyter_output><empty_output><jupyter_text>Some helper functions for matrix creation:<jupyter_code>/* A hyperbolic system matrix resulting from a curl */
curlmat(coord):=genmatrix(
  lambda ([i,j], levi_civita([coord,j,i])),
  3,3)$

vstack:append$

hstack(a,b):=transpose(append(transpose(a),transpose(b)))$

blockmat(a11,a12,a21,a22):=vstack(hstack(a11,a12),hstack(a21,a22))$
n:[nx,ny,nz];<jupyter_output><empty_output><jupyter_text>Next, write down the flux Jacobian on a line with normal `n`:<jupyter_code>max_submat(i):=blockmat(
  zeromatrix(3,3),
  -epsinv*curlmat(i), /* epsinv = 1/epsilon */
  muinv*curlmat(i), /* muinv = 1/mu */
  zeromatrix(3,3)
  )$

max_Asimp:sum(n[i]*max_submat(i),i,1,3);
max_A:subst([epsinv=1/epsilon,muinv=1/mu], max_Asimp);<jupyter_output><empty_output><jupyter_text>Next, diagonalize to obtain the transformation to/from characteristic variables:<jupyter_code>max_invsubst(x):=subst([epsinv=1/epsilon, muinv=1/mu], x)$

[max_V, max_D, max_invV]:max_invsubst(hypdiagonalize(max_Asimp));<jupyter_output><empty_output><jupyter_text>Now find the flux in characteristic variables:<jupyter_code>max_Dinc:subst([1/(sqrt(epsilon)*sqrt(mu))=c], max_D);
max_sflux:hyp_upwind_flux([-c,0,c], max_Dinc);

/* FIXME: max_V should not depend on epsilon and mu, but it does 
   For now, make cp and cm equal. */

max_sflux:subst(
  [cp=1/(sqrt(epsilon)*sqrt(mu)), cm=1/(sqrt(epsilon)*sqrt(mu))], 
  max_sflux);<jupyter_output>solve: dependent equations eliminated: (1 2 11 16 15 12)<jupyter_text>And in conserved variables:<jupyter_code>max_Em:makelist(Em[i],i,1,3)$
max_Ep:makelist(Ep[i],i,1,3)$
max_Hm:makelist(Hm[i],i,1,3)$
max_Hp:makelist(Hp[i],i,1,3)$
max_wm:vstack(max_Em,max_Hm)$
max_wp:vstack(max_Ep,max_Hp)$

max_sminw:hypsimp(max_invV.max_wm)$
max_spinw:hypsimp(max_invV.max_wp)$

max_wflux:fullhypsimp(max_V.ev(max_sflux, [sm=max_sminw,sp=max_spinw]));
max_stronglocalpart:max_A.max_wm;

max_strongwflux:max_stronglocalpart-max_wflux;<jupyter_output><empty_output><jupyter_text>Check against value from [the literature](https://doi.org/10.1016/0010-4655(91)90199-U):<jupyter_code>max_Z:sqrt(mu/epsilon)$
max_Y:sqrt(epsilon/mu)$

max_knownstrongwflux:ratsimp(vstack(
  -1/(2*epsilon)
  *(crossprod(n,(max_Hm-max_Hp)-1/max_Z*crossprod(n,max_Em-max_Ep))),
  1/(2*mu)
  *(crossprod(n,(max_Em-max_Ep)+1/max_Y*crossprod(n,max_Hm-max_Hp)))
  ));

assert(norm_2_squared(hypsimp(
  (max_strongwflux)
  -max_knownstrongwflux))=0);<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup