.editorconfig

+# https://editorconfig.org/
+# https://github.com/editorconfig/editorconfig-vim 
+# https://github.com/editorconfig/editorconfig-emacs 
+
+root = true
+
+[*]
+indent_style = space
+end_of_line = lf
+charset = utf-8
+trim_trailing_whitespace = true
+insert_final_newline = true
+
+[*.py]
+indent_size = 4
+
+[*.rst]
+indent_size = 4
+
+[*.cpp]
+indent_size = 2
+
+[*.hpp]
+indent_size = 2
+
+# There may be one in doc/
+[Makefile]
+indent_style = tab
+
+# https://github.com/microsoft/vscode/issues/1679
+[*.md]
+trim_trailing_whitespace = false

.github/dependabot.yml

+version: 2
+updates:
+    # Set update schedule for GitHub Actions
+    - package-ecosystem: "github-actions"
+      directory: "/"
+      schedule:
+          interval: "weekly"
+
+# vim: sw=4

.github/workflows/autopush.yml

+name: Gitlab mirror
+on:
+    push:
+        branches:
+        - main
+
+jobs:
+    autopush:
+        name: Automatic push to gitlab.tiker.net
+        if: startsWith(github.repository, 'inducer/')
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   run: |
+                curl -L -O https://tiker.net/ci-support-v0
+                . ./ci-support-v0
+                mirror_github_to_gitlab
+
+            env:
+                GITLAB_AUTOPUSH_KEY: ${{ secrets.GITLAB_AUTOPUSH_KEY }}
+
+# vim: sw=4

.github/workflows/ci.yml

+name: CI
+on:
+    push:
+        branches:
+        - main
+    pull_request:
+    schedule:
+        - cron:  '17 3 * * 0'
+
+concurrency:
+    group: ${{ github.head_ref || github.ref_name }}
+    cancel-in-progress: true
+
+jobs:
+    ruff:
+        name: Ruff
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   name: "Main Script"
+            run: |
+                pipx install ruff
+                ruff check
+
+    mypy:
+        name: Mypy
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   name: "Main Script"
+            run: |
+                curl -L -O https://tiker.net/ci-support-v0
+                . ./ci-support-v0
+                build_py_project_in_conda_env
+                python -m pip install mypy
+                ./run-mypy.sh
+
+    typos:
+        name: Typos
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   uses: crate-ci/typos@master
+
+    pylint:
+        name: Pylint
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   name: "Main Script"
+            run: |
+                echo "- matplotlib" >> .test-conda-env-py3.yml
+                echo "- scipy" >> .test-conda-env-py3.yml
+                echo "-------------------------------------------"
+                cat $CONDA_ENVIRONMENT
+                echo "-------------------------------------------"
+                USE_CONDA_BUILD=1
+
+                curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-pylint.sh
+
+                . ./prepare-and-run-pylint.sh "$(basename $GITHUB_REPOSITORY)" test/*.py \
+                    $(find examples -name '*.py')
+
+    pytest3:
+        name: Pytest on Py3
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   name: "Main Script"
+            run: |
+                curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh
+                . ./build-and-test-py-project-within-miniconda.sh
+
+    pyexamples3:
+        name: Examples on Py3
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   name: "Main Script"
+            run: |
+                curl -L -O https://tiker.net/ci-support-v0
+                . ci-support-v0
+                build_py_project_in_conda_env
+                run_examples
+
+                python wave/wave-op-mpi.py --lazy
+                python wave/wave-op-mpi.py --lazy --quad --nonaffine
+                python euler/acoustic_pulse.py --lazy
+                python euler/vortex.py --oi --lazy
+
+                # --oversubscribe is an option for Open MPI (which is what the CI uses)
+                # It allows the CI to succeed even if the CI runner does not
+                # have a sufficient number of cores.
+                mpiexec -np 2 --oversubscribe python wave/wave-op-mpi.py --lazy
+
+                mpiexec -np 2 --oversubscribe python wave/wave-op-mpi.py --numpy
+
+    docs:
+        name: Documentation
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   name: "Main Script"
+            run: |
+                echo "- matplotlib" >> .test-conda-env-py3.yml
+                curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/ci-support.sh
+
+                . ci-support.sh
+                build_py_project_in_conda_env
+                build_docs
+
+    downstream_tests:
+        strategy:
+            matrix:
+                downstream_project: [mirgecom]
+        name: Tests for downstream project ${{ matrix.downstream_project }}
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -   name: "Main Script"
+            env:
+                DOWNSTREAM_PROJECT: ${{ matrix.downstream_project }}
+            run: |
+                curl -L -O -k https://tiker.net/ci-support-v0
+                . ./ci-support-v0
+
+                # https://github.com/inducer/grudge/issues/211
+                export CISUPPORT_PARALLEL_PYTEST=no
+
+                test_downstream "$DOWNSTREAM_PROJECT"
+
+
+# vim: sw=4

.gitignore

@@ -20,3 +20,22 @@ examples/*.pdf
 
 *.vtu
 *.vts
+
+run-debug-*
+
+*.vtu
+*.pvtu
+*.pvd
+*.dot
+*.vtk
+*.silo
+*.dat
+
+.cache
+.pytest_cache
+
+*.json
+
+# pylint stuff
+.pylintrc.yml
+.run-pylint.py

.gitlab-ci.yml

+Python 3 POCL:
+  script:
+  - export PY_EXE=python3
+  - export PYOPENCL_TEST=portable:cpu
+  - export EXTRA_INSTALL="pybind11 numpy mako"
+  - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh
+  - ". ./build-and-test-py-project.sh"
+  tags:
+  - python3
+  - pocl
+  - mpi
+  except:
+  - tags
+  artifacts:
+    reports:
+      junit: test/pytest.xml
+
+Python 3 Intel:
+  script:
+  - export PY_EXE=python3
+  - export EXTRA_INSTALL="pybind11 numpy mako"
+  - source /opt/enable-intel-cl.sh
+  - export PYOPENCL_TEST="intel(r):pu"
+  - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh
+  - ". ./build-and-test-py-project.sh"
+  tags:
+  - python3
+  - pocl
+  - mpi
+  except:
+  - tags
+  artifacts:
+    reports:
+      junit: test/pytest.xml
+
+Python 3 POCL Examples:
+  script:
+  - export PY_EXE=python3
+  - export PYOPENCL_TEST=portable:cpu
+  - export EXTRA_INSTALL="pybind11 numpy mako mpi4py pyvisfile pymetis"
+  - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-py-project-and-run-examples.sh
+  - ". ./build-py-project-and-run-examples.sh"
+  tags:
+  - python3
+  - pocl
+  - large-node
+  except:
+  - tags
+
+Python 3 Intel Examples:
+  script:
+  - export PY_EXE=python3
+  - source /opt/enable-intel-cl.sh
+  - export PYOPENCL_TEST="intel(r):pu"
+  - export EXTRA_INSTALL="pybind11 numpy mako mpi4py pyvisfile pymetis"
+  - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-py-project-and-run-examples.sh
+  - ". ./build-py-project-and-run-examples.sh"
+  tags:
+  - python3
+  - pocl
+  - large-node
+  except:
+  - tags
+
+Python 3 Conda:
+  tags:
+  - linux
+  - large-node
+  script: |
+    export PYOPENCL_TEST=portable:cpu
+    CONDA_ENVIRONMENT=.test-conda-env-py3.yml
+    curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh
+
+    # Shut up ibverbs about fork(), e.g. https://gitlab.tiker.net/inducer/grudge/-/jobs/220796
+    export RDMAV_FORK_SAFE=1
+
+    . ./build-and-test-py-project-within-miniconda.sh
+
+Python 3 Conda Examples:
+  tags:
+  - linux
+  - large-node
+  script: |
+    CONDA_ENVIRONMENT=.test-conda-env-py3.yml
+    curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/ci-support.sh
+    . ci-support.sh
+    build_py_project_in_conda_env
+
+    # Shut up ibverbs about fork(), e.g. https://gitlab.tiker.net/inducer/grudge/-/jobs/220796
+    export RDMAV_FORK_SAFE=1
+
+    run_examples
+
+Documentation:
+  script: |
+    EXTRA_INSTALL="pybind11 numpy matplotlib"
+    curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-docs.sh
+    . ./build-docs.sh
+  tags:
+  - python3
+
+Ruff:
+  script:
+  - pipx install ruff
+  - ruff check
+  tags:
+  - docker-runner
+  except:
+  - tags
+
+Mypy:
+  script: |
+    EXTRA_INSTALL="Cython mpi4py"
+    curl -L -O https://tiker.net/ci-support-v0
+    . ./ci-support-v0
+    build_py_project_in_venv
+    python -m pip install mypy
+    ./run-mypy.sh
+  tags:
+  - python3
+  except:
+  - tags
+
+Pylint:
+  script: |
+    EXTRA_INSTALL="pybind11 make numpy scipy matplotlib mpi4py"
+    curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-pylint.sh
+    . ./prepare-and-run-pylint.sh "$CI_PROJECT_NAME" test/*.py \
+        $(find examples -name '*.py')
+  tags:
+  - python3
+  except:
+  - tags

.pylintrc-local.yml

+- arg: py-version
+  val: '3.10'
+
+- arg: ignore
+  val:
+    - mappers
+- arg: ignored-modules
+  val:
+  - sympy

.test-conda-env-py3.yml

+name: test-conda-env-py3
+channels:
+- conda-forge
+- nodefaults
+dependencies:
+- git
+- numpy
+- libhwloc=2
+# pocl 3.1 needed for SVM functionality
+- pocl>=3.1
+- islpy
+- pyopencl
+- python=3
+- gmsh
+
+# test scripts use ompi-specific arguments
+- openmpi
+- mpi4py

README.rst

 grudge
 ======
 
+.. image:: https://gitlab.tiker.net/inducer/grudge/badges/main/pipeline.svg
+    :alt: Gitlab Build Status
+    :target: https://gitlab.tiker.net/inducer/grudge/commits/main
+.. image:: https://github.com/inducer/grudge/workflows/CI/badge.svg?branch=main&event=push
+    :alt: Github Build Status
+    :target: https://github.com/inducer/grudge/actions?query=branch%3Amain+workflow%3ACI+event%3Apush
+
+..
+    .. image:: https://badge.fury.io/py/grudge.png
+        :alt: Python Package Index Release Page
+        :target: https://pypi.org/project/grudge/
+
 grudge helps you discretize discontinuous Galerkin operators, quickly
 and accurately.
 
 It relies on
 
-* `numpy <http://pypi.python.org/pypi/numpy>`_ for arrays
-* `modepy <http://pypi.python.org/pypi/modepy>`_ for modes and nodes on simplices
-* `meshmdoe <http://pypi.python.org/pypi/meshmdoe>`_ for modes and nodes on simplices
-* `loopy <http://pypi.python.org/pypi/loopy>`_ for fast array operations
-* `pytest <http://pypi.python.org/pypi/pytest>`_ for automated testing
+* `numpy <https://pypi.org/project/numpy>`_ for arrays
+* `modepy <https://pypi.org/project/modepy>`_ for modes and nodes on simplices
+* `meshmode <https://pypi.org/project/meshmode>`_ for modes and nodes on simplices
+* `loopy <https://pypi.org/project/loopy>`_ for fast array operations
+* `pytest <https://pypi.org/project/pytest>`_ for automated testing
 
 and, indirectly,
 
-* `PyOpenCL <http://pypi.python.org/pypi/pyopencl>`_ as computational infrastructure
+* `PyOpenCL <https://pypi.org/project/pyopencl>`_ as computational infrastructure
 
 PyOpenCL is likely the only package you'll have to install
 by hand, all the others will be installed automatically.
 
 .. image:: https://badge.fury.io/py/grudge.png
-    :target: http://pypi.python.org/pypi/grudge
+    :target: https://pypi..org/project/grudge
 
 Resources:
 
-* `documentation <http://documen.tician.de/grudge>`_
-* `wiki home page <http://wiki.tiker.net/Grudge>`_
-* `source code via git <http://gitlab.tiker.net/inducer/grudge>`_
+* `documentation <https://documen.tician.de/grudge>`_
+* `wiki home page <https://wiki.tiker.net/Grudge>`_
+* `source code via git <https://gitlab.tiker.net/inducer/grudge>`_

contrib/maxima/abarbanel-pml.mac

+kill(all);
+load("eigen");
+load("itensor");
+load("diag");
+
+/*
+See
+
+S. Abarbanel und D. Gottlieb, “On the construction and analysis of absorbing
+layers in CEM,” Applied Numerical Mathematics,  vol. 27, 1998, S. 331-340.
+(eq 3.7-3.11)
+
+E. Turkel und A. Yefet, “Absorbing PML
+boundary layers for wave-like equations,”
+Applied Numerical Mathematics,  vol. 27,
+1998, S. 533-557.
+(eq. 4.10) 
+
+*/
+
+/* -------------------------------------------------------------------------- */
+/* Variable declarations */
+/* -------------------------------------------------------------------------- */
+
+coords:[x,y,z];
+allvars:append([t],coords);
+
+max_E:makelist(concat(E,i),i,coords);
+max_H:makelist(concat(H,i),i,coords);
+max_w:append(max_E,max_H);
+
+max_P:makelist(concat(P,i),i,coords);
+max_Q:makelist(concat(Q,i),i,coords);
+aux_w:append(max_P, max_Q);
+
+depends(max_w,allvars);
+depends(aux_w,allvars);
+
+sig:makelist(concat(s,i),i,coords);
+
+depends(sig, coords);
+
+/* -------------------------------------------------------------------------- */
+/* Utilities */
+/* -------------------------------------------------------------------------- */
+
+crossfunc(f):=makelist(
+  sum(sum(
+      levi_civita([i,j,k])*f(j,k),
+   j,1,3),k,1,3),i,1,3)$
+
+curl(x):=crossfunc(lambda([i,j], diff(x[j], coords[i])));
+
+shift(l, amount):=makelist(l[mod(i-amount-1, length(l))+1], i, 1, length(l));
+
+norm(x):=sum(x[i]^2, i, 1, length(x));
+
+
+/* -------------------------------------------------------------------------- */
+/* Operator building */
+/* -------------------------------------------------------------------------- */
+
+/* 
+This here is *NOT* in conservation form
+  u_t + A u_x = 0,
+but in ODE form
+  u_t = (- A) u_x.
+        ---------- rhs
+*/
+
+max_rhs:append(1/epsilon*curl(max_H),-1/mu*curl(max_E));
+
+
+pml_rhs:makelist(0,i,1,6);
+aux_eqn:makelist(0,i,1,6);
+
+for mx:1 thru 3 do
+  block([my,mz,sgn],
+    my:mod(mx-1+1,3)+1,
+    mz:mod(mx-1+2,3)+1,
+    assert(levi_civita([mx,my,mz])=1),
+    pml_rhs[mx]:pml_rhs[mx]
+      +1/epsilon*diff(max_H[mz],coords[my])
+      -sig[my]/epsilon*(2*max_E[mx]+max_P[mx])
+      ,
+    pml_rhs[my]:pml_rhs[my]
+      -1/epsilon*diff(max_H[mz],coords[mx])
+      -sig[mx]/epsilon*(2*max_E[my]+max_P[my])
+      ,
+    pml_rhs[3+mz]:pml_rhs[3+mz]
+      +1/mu*diff(max_E[mx],coords[my])
+      -1/mu*diff(max_E[my],coords[mx])
+      +1/mu*diff(sig[mx]/epsilon,coords[mx])*max_Q[mx]
+      +1/mu*diff(sig[my]/epsilon,coords[my])*max_Q[my]
+      ,
+    aux_eqn[mx]:aux_eqn[mx]
+      -diff(max_P[mx],t)
+      +sig[my]/epsilon*max_E[mx]
+      ,
+    aux_eqn[my]:aux_eqn[my]
+      +sig[mx]/epsilon*max_E[my]
+      ,
+    aux_eqn[3+mx]:aux_eqn[3+mx]
+      -diff(max_Q[mx],t)-sig[mx]/epsilon*max_Q[mx]
+      -max_E[my] -max_E[mz]
+    );
+
+pml_eqn:makelist(diff(max_w[i],t)=pml_rhs[i], i, 1, 6);
+aux_eqn:makelist(aux_eqn[i]=0, i, 1, 6);
+print(expand(covect(pml_eqn)));
+print(expand(covect(aux_eqn)));
+slist:[
+  /*Qx=-Q[y],Qy=Q[x],*/Qz=0,
+  /*Px=P[y],Py=-P[x],*/Pz=0,
+  sz=0,
+  Hx=0,Hy=0,Ez=0,
+  epsilon=1,mu=1];
+
+print(expand(covect(subst(slist,pml_eqn))));
+print(expand(covect(subst(slist,aux_eqn))));
+
+load("em_units.mac");
+u_P:u_sigma*u_E/u_epsilon0*u_t;
+u_Q:u_E*u_t;
+
+assert(u_E/u_t - u_sigma/u_epsilon0*u_P=0);
+assert(u_H/u_t - (1/u_mu0)*(u_sigma/u_epsilon0)/u_x*u_Q=0);
+

contrib/maxima/bgk-flow.mac

+kill(all);
+
+d:2;
+
+load("myhelpers.mac");
+assume(sqrt_rt>0);
+assume(tau>0);
+assume(r>0);
+
+coords:[x,y,z];
+n:makelist(concat(n, coords[k]), k, 1, d);
+
+sqrt_rt:r;
+/*
+w0:a1;
+w1:a2;
+w2:a3;
+w3:a4;
+w4:a5;
+w5:a6;
+*/
+
+/* ------------------------------------------------------------------------ */
+/* do not edit -- automatically generated by BGKFlowOperator */
+
+d:2;
+
+vars: [w0, w1, w2, w3, w4, w5];
+
+depends(vars, coords);
+
+bgk_neg_rhss: [
+sqrt_rt*(diff(w1, x) + diff(w2, y)),
+sqrt_rt*(diff(w0, x) + diff(w3, y) + sqrt(2)*diff(w4, x)),
+sqrt_rt*(diff(w0, y) + diff(w3, x) + sqrt(2)*diff(w5, y)),
+sqrt_rt*(diff(w1, y) + diff(w2, x)) + 1/tau*(w3 + (-1)*w1*w2/w0),
+sqrt(2)*sqrt_rt*diff(w1, x) + 1/tau*(w4 + (-1)*w1*w1/(sqrt(2)*w0)),
+sqrt(2)*sqrt_rt*diff(w2, y) + 1/tau*(w5 + (-1)*w2*w2/(sqrt(2)*w0))
+];
+
+/* ------------------------------------------------------------------------ */
+
+
+bgk_A:genmatrix(
+  lambda([i,j],
+    sum(n[k]*diff(bgk_neg_rhss[i], diff(vars[j], coords[k])), k, 1, d)),
+  length(vars), length(vars));
+
+/* diagonalize system ------------------------------------------------------- */
+[bgk_V, bgk_D, bgk_invV]:hypdiagonalize(bgk_A);
+
+/* auxiliary variables ------------------------------------------------------ */
+bgk_wm:makelist(concat(vars[i], m), i, 1, length(vars));
+bgk_wp:makelist(concat(vars[i], p), i, 1, length(vars));
+bgk_sminw:hypsimp(bgk_invV.bgk_wm);
+bgk_spinw:hypsimp(bgk_invV.bgk_wp);
+
+/* rankine-hugoniot flux ---------------------------------------------------- */
+bgk_sflux:hyp_upwind_flux(
+  [-sqrt(3)*sqrt_rt,-sqrt_rt,0,sqrt_rt,sqrt(3)*sqrt_rt], 
+  bgk_D);
+
+bgk_wflux:fullhypsimp(bgk_V.ev(bgk_sflux, [sm=bgk_sminw,sp=bgk_spinw]));
+
+bgk_stronglocalpart:subst(r=rm, bgk_A).bgk_wm;
+
+bgk_strongwflux:fullhypsimp(bgk_stronglocalpart-bgk_wflux);
+
+bgk_strongwflux:fullhypsimp(subst([rm=r,rp=r], bgk_strongwflux));
+
+/* display flux expression -------------------------------------------------- */
+display2d:false;
+
+dispsubst:append(
+  makelist(concat(w, i, p)=w[i][ext], i, 0, length(vars)-1),
+  makelist(concat(w, i, m)=w[i][int], i, 0, length(vars)-1),
+  makelist(n[i]=normal[i-1], i, 1, d)
+  );
+
+print(makelist(
+  subst(dispsubst, bgk_strongwflux[i,1]),
+  i, 1, length(bgk_strongwflux)));
+

contrib/maxima/bilinearmap.mac

+kill(all);
+lflatten(l):=apply(append, l);
+dims:3;
+origpoints: [[-1],[1]];
+for i:2 thru dims do
+  origpoints:lflatten(
+    makelist(
+      makelist(
+        append(opi, [j])
+        ,j, [-1, 1]
+        ),
+      opi, origpoints));
+points : makelist(ascii(97+i), i, 0, length(origpoints)-1);
+
+vars:makelist(concat(v,ascii(97+i)), i, 0, dims-1);
+mapargs:append(
+  [1], vars,
+  lflatten(
+    makelist(
+      makelist(vars[i]*vars[j], j, i+1, length(vars)), 
+      i, 1, length(vars)
+    )
+  )
+);
+
+/* Idea: (x,y)^T = MAT*(1,u,v,uv) */
+mapmat:genmatrix(mm, dims, length(mapargs));
+
+f:mapmat.mapargs;
+
+myjacobian(f, vars):=genmatrix(
+  lambda([i,j], diff(f[i], vars[j])),
+  length(f), length(vars));
+
+fjac:myjacobian(f, vars);
+
+print("bilinear map jacobian:");
+print(fjac);
+print("bilinear map jacobian det:");
+print(expand(determinant(fjac)));

contrib/maxima/cns.mac

+/* compare formulation in JSH/TW with 
+ * https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Compressible_flow_of_Newtonian_fluids
+ */
+
+kill(all);
+d:2;
+
+div(vec):=sum(diff(vec[i], coords[i]), i, 1, length(vec));
+grad(f):=makelist(diff(f, coords[i]), i, 1, d);
+vlaplace(vec):=makelist(div(grad(vec[i])), i, 1, length(vec));
+
+uvec:makelist([u,v,w][i], i, 1, d);
+coords:makelist([x,y,z][i], i, 1, d);
+
+depends(uvec, coords);
+
+dudx:genmatrix(
+  lambda([i,j], diff(uvec[i], coords[j])),
+  d, d);
+
+delta(i, j):=if i = j then 1 else 0;
+
+tau:mu*genmatrix(
+  lambda([i,j], 
+    dudx[i,j] + dudx[j,i] - 2/3 * delta(i,j) * mat_trace(dudx)),
+  d, d);
+
+rhou_rhs:makelist(div(tau[i]), i, 1, d);
+
+muv:0;
+rhou_rhs2:(1/3*mu+muv)*grad(div(uvec)) + mu*vlaplace(uvec);
+
+/*
+print(rhou_rhs);
+print(rhou_rhs2);
+*/
+print(ratsimp(rhou_rhs-rhou_rhs2));
+
+e_rhs:div(uvec.tau);
+
+e_rhs2:(1/3*mu+muv)*grad(div(uvec)) + mu*makelist(div(uvec*grad(vec[i])), i, 1, d);
+
+print(e_rhs);
+print(e_rhs2);
+print(ratsimp(e_rhs-e_rhs2));

contrib/maxima/eclean.mac

+kill(all);
+
+load("ecleanbase.mac");
+
+print(extract_diagonal(eclean_D));
+
+/* ------------------------------------------------------------------------- */
+/* flux */
+/* ------------------------------------------------------------------------- */
+
+eclean_Dinc:ratsubst(c,1/(sqrt(%epsilon)*sqrt(%mu)), eclean_D);
+eclean_sflux:hyp_upwind_flux([-%chi*c,-c,0,c,%chi*c], eclean_Dinc);
+
+/*
+print("e-clean system flux in terms of characteristic variables:");
+print(eclean_sflux);
+*/
+
+/* FIXME: eclean_V should not depend on epsilon and mu, but it does 
+   For now, make cp and cm equal. */
+
+eclean_sflux:subst(
+  [cp=1/(sqrt(%epsilon)*sqrt(%mu)), 
+   cm=1/(sqrt(%epsilon)*sqrt(%mu)),
+   %chip=%chi,
+   %chim=%chi], 
+  eclean_sflux);
+
+eclean_wflux:fullhypsimp(eclean_V.ev(eclean_sflux, 
+  [sm=eclean_sminw,sp=eclean_spinw]));
+
+/*
+print("e-clean system weak flux in terms of physical variables:");
+print(eclean_wflux);
+*/
+
+eclean_strongwflux:eclean_A.eclean_wm - eclean_wflux;
+
+print("e-clean system strong flux in terms of physical variables:");
+print(eclean_strongwflux);
+
+eclean_maxstrongdiff:fullhypsimp(
+  eclean_strongwflux
+  -append(max_strongwflux,matrix([0]))
+  );
+
+/*
+print("e-clean additional strong-form flux wrt Maxwell system:");
+print(eclean_maxstrongdiff);
+*/
+
+eclean_simple_maxstrongdiff:1/2*append(
+  %chi*c*n*(%phi[m]-%phi[p] - n.(max_Em-max_Ep)),
+  [0,0,0],
+  [c*%chi*(-(%phi[m]-%phi[p]) + n.(max_Em-max_Ep))]);
+
+assert(norm_2_squared(hypsimp(
+  eclean_maxstrongdiff
+  - ev(eclean_simple_maxstrongdiff, [c=1/sqrt(%epsilon*%mu)])))=0);
+  
+/* ------------------------------------------------------------------------- */
+/* Radiation BCs */
+/* ------------------------------------------------------------------------- */
+
+eclean_radbdryspinw:makelist(
+  if eclean_D[i,i] >= 0 then eclean_sminw[i,1] else 0,
+i, 1, length(eclean_D))$
+
+/*
+print("Radiation boundary condition for E-divclean system:");
+print(fullhypsimp(eclean_V.eclean_radbdryspinw));
+*/
+
+/* ------------------------------------------------------------------------- */
+/* Better PEC */
+/* ------------------------------------------------------------------------- */
+/* PEC means: n x E = 0. (no tangential field)
+
+For normal Maxwell PEC, you prescribe E+ = - E-, even though the 
+normal component is wrong (it should be equal to the E- normal component).
+That doesn't matter because the Maxwell flux only looks at the tangential
+component. But here, it suddenly ends up mattering, so we need to prescribe
+the normal component correctly.
+*/
+
+eclean_simplepecwp:covect(append(max_pecbdrywp, [%phi[m]]));
+eclean_betterpecwp:covect(append(
+    normalcomp(max_Em)-tangentialcomp(max_Em), max_Hm, [-%phi[m]]));
+eclean_betterpecsp:fullhypsimp(
+  eclean_invV.ev(eclean_betterpecwp, 
+    [Em=eclean_Emins, Hm=eclean_Hmins, %phi[m]=eclean_phimins]));
+
+/*
+print("Better PEC condition in characteristic variables:");
+print(eclean_betterpecsp);
+*/
+
+print("Better PEC condition in physical variables:");
+print(eclean_betterpecwp);
+
+
+eclean_flux_at_betterpec:fullhypsimp(ev(append(max_strongwflux,matrix([0]))/*eclean_strongwflux*/,
+  [Ep=subrange(eclean_betterpecwp,1,3),
+   Hp=subrange(eclean_betterpecwp,4,6),
+   %phi[p]=eclean_betterpecwp[7,1]]));
+
+eclean_flux_at_simplepec:fullhypsimp(ev(append(max_strongwflux,matrix([0]))/*eclean_strongwflux*/,
+  [Ep=subrange(eclean_simplepecwp,1,3),
+   Hp=subrange(eclean_simplepecwp,4,6),
+   %phi[p]=eclean_simplepecwp[7,1]]));
+
+assert(norm_2_squared(hypsimp(eclean_flux_at_betterpec-eclean_flux_at_simplepec))=0);
+
+/* ------------------------------------------------------------------------- */
+
+eclean_munzpecwithphi(phival):=fullhypsimp(eclean_invV.
+        vstack(vstack(-eclean_Emins, eclean_Hmins), [phival]));
+
+/*
+print("Munz et al PEC boundary with phi+=0:");
+print(eclean_munzpecwithphi(0));
+print("Munz et al PEC boundary with phi+=phi-:");
+print(eclean_munzpecwithphi(eclean_phimins));
+print("Munz et al PEC boundary with phi+=-phi-:");
+print(eclean_munzpecwithphi(-eclean_phimins));
+*/
+
+/* ------------------------------------------------------------------------- */
+/* chi-radiation + PEC BC */
+/* ------------------------------------------------------------------------- */
+eclean_chiradbdryspinw:[
+  -eclean_sminw[3,1],
+  -eclean_sminw[4,1],
+  -eclean_sminw[1,1],
+  -eclean_sminw[2,1],
+  0*eclean_sminw[5,1],
+  eclean_sminw[6,1],
+  eclean_sminw[7,1]
+  ]$
+
+eclean_chiradwp:fullhypsimp(eclean_V.eclean_chiradbdryspinw);
+
+print("PEC+chirad BC for div E cleaning system");
+print(eclean_chiradwp);
+
+eclean_simple_expr_chiradwp:append(
+  -max_Em+3/2*n*(n.max_Em) + 1/2*%phi[m]*n,
+  max_Hm,
+  1/2*[%phi[m]+max_Em.n]
+  );
+
+eclean_diff_chiradwp:hypsimp(eclean_chiradwp
+    - subst([c=1/sqrt(%epsilon*%mu)], eclean_simple_expr_chiradwp));
+    
+assert(norm_2_squared(eclean_diff_chiradwp)=0);
+  
+/*
+print("Limit of PEC+chirad BC as %phi, %chi -> 0");
+print(fullhypsimp(limit(limit(
+  subst([%phi[m]=phi],eclean_chiradwp),phi,0),%chi,0)));
+  */
+
+eclean_strongw_chirad_flux:fullhypsimp(ev(eclean_strongwflux,
+  [Ep=subrange(eclean_chiradwp,1,3),
+   Hp=subrange(eclean_chiradwp,4,6),
+   %phi[p]=eclean_chiradwp[7,1]]));
+print("Flux at PEC+chirad:");
+print(eclean_strongw_chirad_flux);
+
+eclean_strongw_pec_flux:fullhypsimp(ev(eclean_strongwflux,
+  [Ep=subrange(eclean_betterpecwp,1,3),
+   Hp=subrange(eclean_betterpecwp,4,6),
+   %phi[p]=eclean_betterpecwp[7,1]]));
+print("Flux at pure PEC:");
+print(eclean_strongw_pec_flux);
+
+eclean_strongw_flux_diff:fullhypsimp(
+  eclean_strongw_pec_flux-eclean_strongw_chirad_flux);
+
+print("PEC Flux - Chirad Flux");
+print(eclean_strongw_flux_diff);
+
+
+/*
+f1:subst([
+  Normal(0)=nx,Normal(1)=ny,Normal(2)=nz,
+  Int[0]=Em[1],Int[1]=Em[2],Int[2]=Em[3],
+  Ext[0]=Ep[1],Ext[1]=Ep[2],Ext[2]=Ep[3],
+  Int[3]=Hm[1],Int[4]=Hm[2],Int[5]=Hm[3],
+  Ext[3]=Hp[1],Ext[4]=Hp[2],Ext[5]=Hp[3],
+  Int[6]=%phi[m],Ext[6]=%phi[p]] 
+  ,f0)
+*/

contrib/maxima/ecleanbase.mac

+load("maxwellbase.mac");
+
+assume(%chi>0);
+
+/*
+eclean_Asimp:blockmat(
+  max_Asimp,
+  vstack(epsinv*muinv*covect(n),constmatrix(3,1,0)),
+  hstack(%chi^2*n,constmatrix(1,3,0)),
+  zeromatrix(1,1)
+);
+*/
+
+eclean_Asimp:blockmat(
+  max_Asimp,
+  vstack(%chi*sqrt(epsinv*muinv)*covect(n),constmatrix(3,1,0)),
+  hstack(%chi*sqrt(epsinv*muinv)*n,constmatrix(1,3,0)),
+  zeromatrix(1,1)
+);
+
+[eclean_V, eclean_D, eclean_invV]:max_invsubst(hypdiagonalize(eclean_Asimp));
+eclean_A:max_invsubst(eclean_Asimp);
+
+eclean_wm:vstack(max_wm,[%phi[m]]);
+eclean_wp:vstack(max_wp,[%phi[p]]);
+
+eclean_sm:makelist(sm[i],i,1,length(eclean_D));
+eclean_sp:makelist(sp[i],i,1,length(eclean_D));
+
+eclean_sminw:hypsimp(eclean_invV.eclean_wm);
+eclean_spinw:hypsimp(eclean_invV.eclean_wp);
+
+eclean_wmins:hypsimp(eclean_V.eclean_sm);
+eclean_wpins:hypsimp(eclean_V.eclean_sp);
+
+eclean_Emins:makelist(eclean_wmins[i,1],i,1,3);
+eclean_Hmins:makelist(eclean_wmins[i,1],i,4,6);
+eclean_phimins:eclean_wmins[7,1];

contrib/maxima/ehclean.mac

+kill(all);
+load("ehcleanbase.mac");
+
+/* ------------------------------------------------------------------------- */
+/* flux */
+/* ------------------------------------------------------------------------- */
+
+ehclean_Dinc:ratsubst(c,1/(sqrt(%epsilon)*sqrt(%mu)), ehclean_D);
+print(ehclean_Dinc);
+ehclean_sflux:hyp_upwind_flux([-%chi*c,-c,c,%chi*c], ehclean_Dinc);
+
+print("eh-clean system flux in terms of characteristic variables:");
+print(covect(ehclean_sflux));
+
+/* FIXME: ehclean_V should not depend on epsilon and mu, but it does 
+   For now, make cp and cm equal. */
+/*
+ehclean_sflux:subst(
+  [cp=1/(sqrt(%epsilon)*sqrt(%mu)), 
+   cm=1/(sqrt(%epsilon)*sqrt(%mu)),
+   %chip=%chi,
+   %chim=%chi], 
+  ehclean_sflux);
+
+print("e-clean system flux in terms of physical variables:");
+ehclean_wflux:fullhypsimp(ehclean_V.ev(ehclean_sflux, 
+  [sm=ehclean_sminw,sp=ehclean_spinw]));
+print(ehclean_wflux);
+*/

contrib/maxima/ehcleanbase.mac

+load("maxwellbase.mac");
+
+assume(%chi>0);
+
+ehclean_Asimp:blockmat(
+  max_Asimp,
+  hstack(
+    vstack(epsinv*muinv*covect(n),zeromatrix(3,1)),
+    vstack(zeromatrix(3,1),epsinv*muinv*covect(n))
+  ),
+
+  vstack(
+    hstack(%chi^2*n,zeromatrix(1,3)),
+    hstack(zeromatrix(1,3),%chi^2*n)
+  ),
+  zeromatrix(2,2)
+);
+ehclean_A:max_invsubst(ehclean_Asimp);
+[ehclean_V, ehclean_D, ehclean_invV]:max_invsubst(hypdiagonalize(ehclean_Asimp));

contrib/maxima/em-modes.mac

+kill(all);
+load("eigen");
+load("itensor");
+load("diag");
+
+/* -------------------------------------------------------------------------- */
+/* Utilities */
+/* -------------------------------------------------------------------------- */
+
+curl(x):=crossfunc(lambda([i,j], diff(x[j], coords[i])));
+div(f):=sum(diff(f[j], coords[j]), j, 1, length(coords));
+
+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]));
+
+/* -------------------------------------------------------------------------- */
+/* Variable declarations */
+/* -------------------------------------------------------------------------- */
+
+coords:[x,y,z];
+allvars:append([t],coords);
+
+mu: 1;
+epsilon: 1;
+c: 1/sqrt(mu*epsilon);
+
+epsilon_0: 1;
+mu_0: 1;
+c_0: 1/sqrt(mu_0*epsilon_0);
+
+max_B(max_H):= mu*max_H;
+max_D(max_E):= epsilon*max_E;
+
+/* SI conventions, assumed time dep: exp(%i*omega*t) */
+faraday(max_E, max_H):= curl(max_E) + %i * omega * max_B(max_H);
+ampere(max_E, max_H):= curl(max_B(max_H)) - %i * omega / c_0**2 * max_E;
+
+div_e(max_E, max_H):= div(max_E);
+div_h(max_E, max_H):= div(max_H);
+
+maxwell_pde(max_E, max_H):=append(
+    faraday(max_E, max_H),
+    ampere(max_E, max_H),
+    [div_e(max_E, max_H), div_h(max_E, max_H)]);
+
+/*
+spatial_deps:[
+  exp(%i*m*x)*exp(%i*n*y),
+  exp(%i*m*x)*exp(-%i*n*y),
+  exp(-%i*m*x)*exp(%i*n*y),
+  exp(-%i*m*x)*exp(-%i*n*y)
+  ];
+*/
+
+spatial_deps:[
+  exp(+%i*m*x)*exp(+%i*n*y)*exp(+%i*l*z),
+  exp(+%i*m*x)*exp(+%i*n*y)*exp(-%i*l*z),
+  exp(+%i*m*x)*exp(-%i*n*y)*exp(+%i*l*z),
+  exp(+%i*m*x)*exp(-%i*n*y)*exp(-%i*l*z)
+
+  /*
+  exp(-%i*m*x)*exp(+%i*n*y)*exp(+%i*l*z),
+  exp(-%i*m*x)*exp(+%i*n*y)*exp(-%i*l*z),
+  exp(-%i*m*x)*exp(-%i*n*y)*exp(+%i*l*z),
+  exp(-%i*m*x)*exp(-%i*n*y)*exp(-%i*l*z)
+  */
+  ];
+
+declare(m, integer, n, integer, l, integer);
+
+get_maxwell_solution(spatial_dep):=block([
+  max_B, max_D, coeffs, max_E, max_H, faraday, ampere, div_e, div_h, maxwell_pde, soln
+  ],
+  max_B: mu*max_H,
+  max_D: epsilon*max_E,
+
+  coeffs: [
+    Exr, Eyr, Ezr, Hxr, Hyr, Hzr,
+    Exi, Eyi, Ezi, Hxi, Hyi, Hzi
+    ],
+
+  max_E: [
+      (Exr+%i*Exi)*spatial_dep,
+      (Eyr+%i*Eyi)*spatial_dep,
+      (Ezr+%i*Ezi)*spatial_dep
+      ],
+  max_H: [
+      (Hxr+%i*Hxi)*spatial_dep,
+      (Hyr+%i*Hyi)*spatial_dep,
+      (Hzr+%i*Hzi)*spatial_dep
+      ],
+
+  soln:solve(
+    maxwell_pde(max_E, max_H),
+    append(coeffs, [omega])),
+
+  family1: soln[1],
+  omega_eq: family1[length(family1)],
+  assert(lhs(omega_eq) = omega),
+
+  [subst(family1, [max_E, max_H]), rhs(omega_eq)]
+  );
+
+maxwell_solutions:makelist(
+  get_maxwell_solution(spatial_deps[i]),
+  i, 1, length(spatial_deps));
+
+omegas:makelist(
+  maxwell_solutions[i][2],
+  i, 1, length(maxwell_solutions));
+display(omegas);
+
+max_E:ratsimp(sum(
+  maxwell_solutions[i][1][1],
+  i, 1, length(maxwell_solutions)));
+max_H:ratsimp(sum(
+  maxwell_solutions[i][1][2],
+  i, 1, length(maxwell_solutions)));
+
+print("Check Maxwell:");
+print(ratsimp(subst([omega=omegas[1]],maxwell_pde(max_E,max_H))));
+
+pec_bcs:append(
+    realpart(crossprod([-1,0,0], subst(x=0, max_E))),
+    realpart(crossprod([1,0,0], subst(x=%pi, max_E))),
+    realpart(crossprod([0,-1,0], subst(y=0, max_E))),
+    realpart(crossprod([0,1,0], subst(y=%pi, max_E))),
+    [
+      realpart([-1,0,0].subst(x=0, max_H)),
+      realpart([1,0,0].subst(x=%pi, max_H)),
+      realpart([0,-1,0].subst(y=0, max_H)),
+      realpart([0,1,0].subst(y=%pi, max_H))
+    ]);
+
+freevars: sublist(
+    listofvars([max_E, max_H]),
+    lambda([rvar], substring(string(rvar),1,2) = "%"));
+
+ev_pec_bcs:append(
+  subst([x=0, y=0, z=0], pec_bcs),
+  subst([x=0, y=0, z=%pi], pec_bcs),
+  subst([x=0, y=%pi, z=%pi], pec_bcs)
+  );
+
+/*
+Fails:
+
+pec_soln:linsolve(pec_bcs, freevars);
+*/

contrib/maxima/em_units.mac

+kill(all);
+load("myhelpers.mac");
+assume(m>0);
+assume(N>0);
+assume(s>0);
+
+N:kg*m/s^2;
+
+C : A*s;
+J : N*m;
+V : J/C;
+F : C/V;
+T : V*s/m^2;
+N : kg*m/s^2;
+
+u_c0 : m/s;
+u_mu0 : N / A^2;
+u_epsilon0 : 1/(u_c0^2*u_mu0);
+
+
+c0 : 299792458 * u_c0;
+mu0 : 4*%pi * 10^(-7) * u_mu0;
+epsilon0 : 1/(c0^2*mu0);
+
+u_E : V/m;
+u_D : C/m^2;
+u_B : T;
+u_H : A/m;
+u_J : A/m^2;
+u_t : s;
+u_x : m;
+
+u_sigma : u_J/u_E;
+
+u_ndE : u_E;
+u_ndH : sqrt(u_mu0/u_epsilon0) * u_H;
+

contrib/maxima/euler.mac

+kill(all);
+
+load("myhelpers.mac");
+
+assume(%gamma>1);
+
+d:2;
+usyms:[u,v,w];
+coords:[x,y,z];
+uvec:makelist(usyms[i], i, 1, d);
+rhouvec:rho*uvec;
+E_expr:p/(%gamma-1)+rho/2*uvec.uvec;
+
+p_solved:rhs(solve(E=E_expr,p)[1]);
+p_used: p; /* p or p_solved */
+
+/* fluxes ------------------------------------------------------------------- */
+vars:append([rho, E], rhouvec);
+depends(append([rho, E, p], uvec), [x,y,z,t]);
+
+rho_flux:rhouvec;
+E_flux:uvec*(E+p_used);
+rhou_flux:makelist(makelist(
+  rhouvec[i]*uvec[j] + if i = j then p_used else 0,
+  j, 1, d), i, 1, d);
+
+all_fluxes:makelist(
+  vstack([rho_flux[i]], [E_flux[i]], rhou_flux[i]),
+  i, 1, d);
+
+euler_eqns:makelist(
+  -'diff(vars[nr],t)=sum('diff(all_fluxes[i][nr], coords[i]), i, 1, d),
+  nr, 1, length(vars));
+
+/* linearization ------------------------------------------------------------ */
+u0vec:makelist(concat(usyms[i], 0), i, 1, d);
+duvec:makelist(concat('d, usyms[i]), i, 1, d);
+drhouvec:rho0*duvec+drho*u0vec;
+
+assume(%gamma>1);
+assume(p0>0);
+assume(rho0>0);
+
+zero_subst:append(
+    [rho=rho0, p=p0], 
+    makelist(usyms[i] =concat(usyms[i], 0), i, 1, d));
+lin_subst:append(
+    [rho=rho0+drho, p=p0+dp], 
+    makelist(usyms[i] =concat(usyms[i], 0) + duvec[i], i, 1, d));
+
+E0_expr:subst(zero_subst,E_expr);
+dE:subst(lin_subst,E_expr)-E0_expr;
+
+lin_subst:append(lin_subst, [E=E0+dE]);
+
+kill_second_order_terms(x):=block(
+  for lv1 in all_lin_vars do
+    for lv2 in all_lin_vars do
+      block(
+        x:ratsubst(0,lv1*lv2, x),
+        for co in coords do
+          x:ratsubst(0,diff(lv1, co)*lv2, x)),
+  x);
+
+lin_vars:append([drho, dE], drhouvec);
+all_lin_vars:append(duvec, [drho, dp]);
+
+depends(all_lin_vars, [x,y,z,t]);
+
+lin_euler_eqns:kill_second_order_terms(makelist(
+  /* ref solution is time-invariant*/
+  -diff(lin_vars[nr],t)
+  =sum(diff(subst(lin_subst, all_fluxes[i][nr]), coords[i]), i, 1, d),
+  nr, 1, length(vars)));
+
+lin_euler_eqns:ratsimp(lin_euler_eqns);
+
+/* convert to primitive variables */
+ddrho_dt: rhs(solve(lin_euler_eqns[1],diff(drho,t))[1]);
+
+lin_euler_eqns_p:makelist(lin_euler_eqns[i], i, 1, length(lin_euler_eqns));
+
+for i:2 thru d+2 do
+  lin_euler_eqns_p[i]:subst(diff(drho,t)=ddrho_dt, lin_euler_eqns_p[i]),
+
+for i:1 thru d do
+    lin_euler_eqns_p[2+i]:-solve(lin_euler_eqns_p[2+i], diff(duvec[i], t))[1];
+dduvec_dt: makelist(
+  rhs(solve(lin_euler_eqns_p[i+2],diff(duvec[i],t))[1]),
+  i, 1, d);
+
+lin_euler_eqns_p[2]:subst(E0=E0_expr, lin_euler_eqns_p[2]);
+for i:1 thru d do
+  lin_euler_eqns_p[2]:subst(diff(duvec[i],t)=dduvec_dt[i], lin_euler_eqns_p[2]);
+lin_euler_eqns_p[2]:-solve(lin_euler_eqns_p[2], diff(dp, t))[1];
+
+lin_euler_eqns_p:kill_second_order_terms(lin_euler_eqns_p);
+
+/* matrix building ---------------------------------------------------------- */
+prim_lin_vars:append([drho,dp],duvec);
+n:makelist(concat(n, coords[k]), k, 1, d);
+euler_mat:genmatrix(
+  lambda([i,j],
+    sum(n[k]
+      *diff(rhs(lin_euler_eqns_p[i]), diff(prim_lin_vars[j], coords[k])),
+      k, 1, d)),
+  length(lin_vars), length(lin_vars));
+
+/* diagonalize system ------------------------------------------------------- */
+[euler_V, euler_D, euler_invV]:hypdiagonalize(euler_mat);
+
+rel_D:ratsimp(euler_D - n.u0vec*ident(d+2));
+
+/* auxiliary variables, external and internal states ------------------------ */
+c0:sqrt(%gamma*p0/rho0);
+
+euler_wm:makelist(concat(prim_lin_vars[i], m), i, 1, length(prim_lin_vars));
+euler_wp:makelist(concat(prim_lin_vars[i], p), i, 1, length(prim_lin_vars));
+
+euler_sminw:hypsimp(euler_invV.euler_wm);
+euler_spinw:hypsimp(euler_invV.euler_wp);
+
+dumvec:makelist(euler_wm[i+2], i, 1, d);
+dupvec:makelist(euler_wp[i+2], i, 1, d);
+
+/* subsonic outflow bc------------------------------------------------------- */
+euler_outflowbdryspinw:makelist(
+  /* outer state equals free-stream flow, about which we have linearized. 
+  Hence the linearized state variable, which is the difference to free-stream
+  flow, is set to zero. */
+  /* 
+  = 0: convection: from interior
+  > 0: supersonic outflow: from interior
+  < 0: supersonic inflow: from exterior
+  */
+  if rel_D[i,i] >= 0 then euler_sminw[i,1] else 0,
+    i, 1, d+2);
+
+euler_woutflowbdry:fullhypsimp(euler_V.euler_outflowbdryspinw);
+
+euler_known_woutflowbdry:vstack([
+  drhom + n.dumvec*rho0/(2*c0) - dpm/(2*c0^2),
+  c0*rho0*(n.dumvec)/2 + dpm/2],
+  dumvec - n*(n.dumvec)/2 + dpm*n/(2*c0*rho0)
+  );
+
+euler_diff_woutflowbdry:hypsimp(euler_woutflowbdry-euler_known_woutflowbdry);
+assert(norm_2_squared(euler_diff_woutflowbdry)=0);
+
+/* subsonic inflow bc ------------------------------------------------------- */
+euler_inflowbdryspinw:makelist(
+  /* outer state equals free-stream flow, about which we have linearized. 
+  Hence the linearized state variable, which is the difference to free-stream
+  flow, is set to zero. */
+  /* 
+  > 0: supersonic outflow: from interior
+  = 0: convection: from exterior
+  < 0: supersonic inflow: from exterior
+  */
+  if rel_D[i,i] <= 0 then 0 else euler_sminw[i,1],
+    i, 1, d+2);
+
+euler_winflowbdry:fullhypsimp(euler_V.euler_inflowbdryspinw);
+
+euler_known_winflowbdry:vstack([
+  n.dumvec*rho0/(2*c0) + dpm/(2*c0^2),
+  c0*rho0*(n.dumvec)/2 + dpm/2],
+  n*(n.dumvec)/2 + dpm*n/(2*c0*rho0)
+  );
+
+euler_diff_winflowbdry:hypsimp(euler_winflowbdry-euler_known_winflowbdry);
+assert(norm_2_squared(euler_diff_winflowbdry)=0);