Merge branch 'master' into do/lts-baseline

makefile

@@ -125,7 +125,7 @@ EXAMPLE_TEST_DIRS := examples
 
 MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag autodiff hooke \
- multidomain dpg hdiv-linear-solver spde
+ multidomain dpg hdiv-linear-solver spde spinner
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \

miniapps/CMakeLists.txt

@@ -36,3 +36,4 @@ add_subdirectory(parelag)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
 add_subdirectory(hdiv-linear-solver)
+add_subdirectory(spinner)

miniapps/spinner/CMakeLists.txt

@@ -0,0 +1,29 @@
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+list(APPEND SPINNER_SOURCES
+    spinner-lts.cpp
+    cell-solver.hpp)
+list(APPEND SPINNER_HEADERS
+    cell-solver.hpp
+    new-kelly-estimator.hpp)
+
+convert_filenames_to_full_paths(SPINNER_SOURCES)
+convert_filenames_to_full_paths(SPINNER_HEADERS)
+
+set(SPINNER_COMMON_FILES
+    EXTRA_SOURCES ${SPINNER_SOURCES}
+    EXTRA_HEADERS ${SPINNER_HEADERS})
+
+add_mfem_miniapp(spinner-lts
+    MAIN spinner-lts.cpp
+    ${SPINNER_COMMON_FILES}
+    LIBRARIES mfem)

miniapps/spinner/cell-solver.hpp

@@ -0,0 +1,469 @@
+class AbstractCellSolver 
+{
+public:
+    virtual ~AbstractCellSolver() = default;
+
+    virtual void Init(mfem::GridFunction& Vgf, mfem::GridFunction& sgf) const = 0;
+
+    virtual double EvalReaction(const double V, const mfem::Vector &cells_s, int int_point, float t) const = 0;
+    virtual void InternalRHS(mfem::Vector&rhs, double V, mfem::Vector &cells_s) const = 0;
+
+    virtual void FullExplicitStep(mfem::Vector& cells_V, mfem::Vector& cells_s, float t, float Δt) const  = 0;
+    virtual void InternalEulerUpdate(const double V, mfem::Vector &cells_s, int int_point, float t, double Δt) const = 0;
+
+    virtual int InternalDim() const = 0;
+    virtual int HGateIndex() const = 0;
+    
+    virtual void RescaleVoltage(mfem::GridFunction& φₘgf) const = 0;
+};
+
+// Implementation of https://www.frontiersin.org/articles/10.3389/fphys.2019.00721/full
+class PathmanathanCordeiroGrayCellSolver final 
+: public AbstractCellSolver
+{
+private:
+    // Parameters
+    float C_m = 1.0f; // [µF/cm^-2]
+    //float C_m = 0.01f; // [µF/mm^-2]
+    // ------ I_Na -------
+    constexpr static float g_Na = 12.0f;    // [mS/µF]
+    constexpr static float E_m  = -52.244f; // [mV]
+    constexpr static float k_m  = 6.5472f;  // [mV]
+    constexpr static float τ_m  = 0.12f;    // [ms]
+    constexpr static float E_h  = -78.7f;   // [mV]
+    constexpr static float k_h  = 5.93f;    // [mV]
+    constexpr static float δ_h  = 0.799163; // dimensionless
+    constexpr static float τ_h0 = 6.80738;  // [ms]
+    // ------ I_K1 -------
+    constexpr static float g_K1 = 0.73893f;  // [mS/µF]
+    constexpr static float E_z  = -91.9655f; // [mV]
+    constexpr static float k_z  = 12.4997f;  // [mV]
+    // ------ I_to -------
+    float g_to = 0.1688f*1.9;   // [mS/µF]
+    constexpr static float E_r  = 14.3116f;  // [mV]
+    constexpr static float k_r  = 11.462f;   // [mV]
+    constexpr static float E_s  = -47.9286f; // [mV]
+    constexpr static float k_s  = 4.9314f;   // [mV]
+    constexpr static float τ_s  = 9.90669f;  // [ms]
+    // ------ I_CaL -------
+    constexpr static float g_CaL = 0.11503f; // [mS/µF]
+    constexpr static float E_d   = 0.7f;     // [mV]
+    constexpr static float k_d   = 4.3f;     // [mV]
+    constexpr static float E_f   = -15.7f;   // [mV]
+    constexpr static float k_f   = 4.6f;     // [mV]
+    constexpr static float τ_f   = 30.0f;    // [ms]
+    // ------ I_Kr -------
+    constexpr static float g_Kr = 0.056f; // [mS/µF]
+    constexpr static float E_xr = -26.6f; // [mV]
+    constexpr static float k_xr = 6.5f;   // [mV]
+    constexpr static float τ_xr = 334.0f; // [ms]
+    constexpr static float E_y  = -49.6f; // [mV]
+    constexpr static float k_y  = 23.5f;  // [mV]
+    // ------- I_Ks --------
+    constexpr static float g_Ks = 0.008f; // [mS/µF]
+    constexpr static float E_xs = 24.6f;  // [mV]
+    constexpr static float k_xs = 12.1f;  // [mV]
+    constexpr static float τ_xs = 628.0f; // [ms]
+    // ------- Other --------
+    constexpr static float E_Na = 65.0f;  // [mV]
+    constexpr static float E_K  = -85.0f; // [mV]
+    constexpr static float E_Ca = 50.0f;  // [mV]
+
+    // Helper
+    static inline float sigmoid(const float v, const float E_Y, const float k_Y,
+                                const float sign)
+    {
+        return 1.0f / (1.0f + exp(sign * (v - E_Y) / k_Y));
+    }
+
+    inline float rhs_h(const float v, const float H) const
+    {
+        const float τ = (2.0f * τ_h0 * exp(δ_h * (v - E_h) / k_h)) /
+                        (1.0f + exp((v - E_h) / k_h));
+        const float h_inf = sigmoid(v, E_h, k_h, 1.0f);
+
+        const float b = h_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a*H + b;
+    }
+
+    // derivative of (1/(1+exp((x-E)/k))-h)/((2*t*exp(delta*(x-E)/k))/(1+exp((x-E)/k))) in x
+    inline float rhs_hdV(const float v, const float H) const
+    {
+        double term1 = -δ_h/(2.0f*k_h*τ_h0) * (1.0f - H*(exp((v-E_h)/k_h)+1.0f)) * exp(δ_h * (v - E_h) / k_h);
+        double term2 = 1.0f/(2.0f*k_h*τ_h0) * (1.0f/(exp((v-E_h)/k_h)+1.0f) - H) * (exp((v-E_h-(δ_h*(v-E_h)))/k_h));
+        double term3 = exp((v-E_h)/k_h-(δ_h*(v-E_h))/k_h)/(2.0f*k_h*τ_h0*exp((v-E_h)/k_h));
+        return term1 + term2 + term3;
+    }
+
+    inline float rhs_hdh(const float v, const float H) const
+    {
+        const float τ = (2.0f * τ_h0 * exp(δ_h * (v - E_h) / k_h)) /
+                        (1.0f + exp((v - E_h) / k_h));
+        // const float h_inf = sigmoid(v, E_h, k_h, 1.0f);
+
+        // const float b = h_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a;
+    }
+
+    inline float update_h(const float v, const float H, const float Δt) const
+    {
+        const float τ = (2.0f * τ_h0 * exp(δ_h * (v - E_h) / k_h)) /
+                        (1.0f + exp((v - E_h) / k_h));
+        const float h_inf = sigmoid(v, E_h, k_h, 1.0f);
+
+        const float b = h_inf / τ;
+        const float a = -1.0f / τ;
+
+        return std_clamp<float>(exp(a * Δt) * (H + b / a) - b / a, 0.0f,
+                          1.0f);
+    }
+
+    inline float rhs_m(const float v, const float M) const
+    {
+        float τ     = τ_m;
+        float m_inf = sigmoid(v, E_m, k_m, -1.0f);
+
+        float b = m_inf / τ;
+        float a = -1.0f / τ;
+
+        return a*M + b;
+    }
+
+    inline float rhs_mdv(const float v, const float M) const
+    {
+        return exp(-(v-E_m)/k_m)/(k_m*τ_m*(exp(-(v-E_m)/k_m)+1)*(exp(-(v-E_m)/k_m)+1));
+    }
+
+    inline float rhs_mdm(const float v, const float M) const
+    {
+        float τ     = τ_m;
+        // float m_inf = sigmoid(v, E_m, k_m, -1.0f);
+
+        // float b = m_inf / τ;
+        float a = -1.0f / τ;
+
+        return a;
+    }
+
+    inline float update_m(const float v, const float M, const float Δt) const
+    {
+        float τ     = τ_m;
+        float m_inf = sigmoid(v, E_m, k_m, -1.0f);
+
+        float b = m_inf / τ;
+        float a = -1.0f / τ;
+
+        return std_clamp<float>(exp(a * Δt) * (M + b / a) - b / a, 0.0f,
+                          1.0f);
+    }
+
+    inline float rhs_s(const float v, const float S) const
+    {
+        const float τ     = τ_s;
+        const float s_inf = sigmoid(v, E_s, k_s, 1.0f);
+
+        float b = s_inf / τ;
+        float a = -1.0f / τ;
+
+        return a*S + b;
+    }
+
+    inline float rhs_sdv(const float v, const float S) const
+    {
+        return exp((v-E_s)/k_s)/(k_s*τ_s*(exp((v-E_s)/k_s)+1)*(exp((v-E_s)/k_s)+1));
+    }
+
+    inline float rhs_sds(const float v, const float S) const
+    {
+        const float τ     = τ_s;
+        // const float s_inf = sigmoid(v, E_s, k_s, 1.0f);
+
+        // float b = s_inf / τ;
+        float a = -1.0f / τ;
+
+        return a;
+    }
+
+    inline float update_s(const float v, const float S, const float Δt) const
+    {
+        const float τ     = τ_s;
+        const float s_inf = sigmoid(v, E_s, k_s, 1.0f);
+
+        float b = s_inf / τ;
+        float a = -1.0f / τ;
+
+        return std_clamp<float>(exp(a * Δt) * (S + b / a) - b / a, 0.0f,
+                          1.0f);
+    }
+
+    inline float rhs_f(const float v, const float F) const
+    {
+        const float τ     = τ_f;
+        const float f_inf = sigmoid(v, E_f, k_f, 1.0f);
+
+        const float b = f_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a*F + b;
+    }
+
+    inline float rhs_fdv(const float v, const float F) const
+    {
+        return exp((v-E_f)/k_f)/(k_f*τ_f*(exp((v-E_f)/k_f)+1)*(exp((v-E_f)/k_f)+1));
+    }
+
+    inline float rhs_fdf(const float v, const float F) const
+    {
+        const float τ     = τ_f;
+        // const float f_inf = sigmoid(v, E_f, k_f, 1.0f);
+
+        // const float b = f_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a;
+    }
+
+    inline float update_f(const float v, const float F, const float Δt) const
+    {
+        const float τ     = τ_f;
+        const float f_inf = sigmoid(v, E_f, k_f, 1.0f);
+
+        const float b = f_inf / τ;
+        const float a = -1.0f / τ;
+
+        return std_clamp<float>(exp(a * Δt) * (F + b / a) - b / a, 0.0f,
+                          1.0f);
+    }
+
+    inline float rhs_xr(const float v, const float Xr) const
+    {
+        const float τ      = τ_xr;
+        const float xr_inf = sigmoid(v, E_xr, k_xr, -1.0f);
+
+        const float b = xr_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a*Xr + b;
+    }
+
+    inline float rhs_xrdv(const float v, const float Xr) const
+    {
+        return exp(-(v-E_xr)/k_xr)/(k_xr*τ_xr*(exp(-(v-E_xr)/k_xr)+1)*(exp(-(v-E_xr)/k_xr)+1));
+    }
+
+    inline float rhs_xrdxr(const float v, const float Xr) const
+    {
+        const float τ      = τ_xr;
+        // const float xr_inf = sigmoid(v, E_xr, k_xr, -1.0f);
+
+        // const float b = xr_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a;
+    }
+
+    inline float update_xr(const float v, const float Xr, const float Δt) const
+    {
+        const float τ      = τ_xr;
+        const float xr_inf = sigmoid(v, E_xr, k_xr, -1.0f);
+
+        const float b = xr_inf / τ;
+        const float a = -1.0f / τ;
+
+        return std_clamp<float>(exp(a * Δt) * (Xr + b / a) - b / a, 0.0f,
+                          1.0f);
+    }
+
+    inline float rhs_xs(const float v, const float Xs) const
+    {
+        const float τ      = τ_xs;
+        const float xs_inf = sigmoid(v, E_xs, k_xs, -1.0f);
+
+        const float b = xs_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a*Xs + b;
+    }
+
+    inline float rhs_xsdv(const float v, const float Xs) const
+    {
+        return exp(-(v-E_xs)/k_xs)/(k_xs*τ_xs*(exp(-(v-E_xs)/k_xs)+1)*(exp(-(v-E_xs)/k_xs)+1));
+    }
+
+    inline float rhs_xsdxs(const float v, const float Xs) const
+    {
+        const float τ      = τ_xs;
+        // const float xs_inf = sigmoid(v, E_xs, k_xs, -1.0f);
+
+        // const float b = xs_inf / τ;
+        const float a = -1.0f / τ;
+
+        return a;
+    }
+
+    inline float update_xs(const float v, const float Xs, const float Δt) const
+    {
+        const float τ      = τ_xs;
+        const float xs_inf = sigmoid(v, E_xs, k_xs, -1.0f);
+
+        const float b = xs_inf / τ;
+        const float a = -1.0f / τ;
+
+        return std_clamp<float>(exp(a * Δt) * (Xs + b / a) - b / a, 0.0f,
+                          1.0f);
+    }
+
+public:
+    PathmanathanCordeiroGrayCellSolver(double C_m_)
+    : C_m(C_m_)
+    {
+    }
+    ~PathmanathanCordeiroGrayCellSolver() = default;
+
+    void RescaleVoltage(mfem::GridFunction& φₘgf) const override
+    {
+        for(int i=0; i<φₘgf.Size(); i++) {
+            double offset = 1-φₘgf(i);
+            φₘgf(i) = -85.0*offset + (1.0-offset)*-5.0;
+        }
+    }
+
+    int HGateIndex() const override { return 0; };
+
+    void Init(mfem::GridFunction& Vgf, mfem::GridFunction& sgf) const override
+    {
+        assert(sgf.FESpace()->GetOrdering() == mfem::Ordering::byVDIM);
+
+        for (int node = 0; node < Vgf.FESpace()->GetNDofs(); node++)
+        {
+            Vgf(node)    = E_K;
+        }
+
+        for (int node = 0; node < sgf.FESpace()->GetNDofs(); node++)
+        {
+            const auto V = E_K;
+
+            sgf(6 * node + 0) = sigmoid(V, E_h, k_h, 1.0f);
+            sgf(6 * node + 1) = sigmoid(V, E_m, k_m, -1.0f);
+            sgf(6 * node + 2) = sigmoid(V, E_f, k_f, 1.0f);
+            sgf(6 * node + 3) = sigmoid(V, E_s, k_s, 1.0f);
+            sgf(6 * node + 4) = sigmoid(V, E_xs, k_xs, -1.0f);
+            sgf(6 * node + 5) = sigmoid(V, E_xr, k_xr, -1.0f);
+        }
+    }
+
+    void FullExplicitStep(mfem::Vector& cells_V, mfem::Vector& cells_s, float t, float Δt) const override
+    {
+        for (int i = 0; i < cells_V.Size(); i++)
+        {
+            float actual_Δt = Δt;
+
+            const double V = cells_V(i);
+
+            const float h  = cells_s(6 * i + 0);
+            const float m  = cells_s(6 * i + 1);
+            const float f  = cells_s(6 * i + 2);
+            const float s  = cells_s(6 * i + 3);
+            const float xs = cells_s(6 * i + 4);
+            const float xr = cells_s(6 * i + 5);
+
+            // Update gates
+            const float h_  = update_h (V, h,  actual_Δt);
+            const float m_  = update_m (V, m,  actual_Δt);
+            const float f_  = update_f (V, f,  actual_Δt);
+            const float s_  = update_s (V, s,  actual_Δt);
+            const float xs_ = update_xs(V, xs, actual_Δt);
+            const float xr_ = update_xr(V, xr, actual_Δt);
+
+            // Instantaneous gates
+            const float r = sigmoid(V, E_r, k_r, -1.0);
+            const float d = sigmoid(V, E_d, k_d, -1.0);
+            const float z = sigmoid(V, E_z, k_z, 1.0);
+            const float y = sigmoid(V, E_y, k_y, 1.0);
+
+            // Currents
+            const float I_Na  = g_Na * m * m * m * h * h * (V - E_Na);
+            const float I_K1  = g_K1 * z * (V - E_K);
+            const float I_to  = g_to * r * s * (V - E_K);
+            const float I_CaL = g_CaL * d * f * (V - E_Ca);
+            const float I_Kr  = g_Kr * xr * y * (V - E_K);
+            const float I_Ks  = g_Ks * xs * (V - E_K);
+
+            const float I_total =
+                I_Na + I_K1 + I_to + I_CaL + I_Kr + I_Ks;
+
+            // Actual step
+            cells_V(i) = V - actual_Δt * I_total / C_m;
+
+            cells_s(6 * i + 0) = h_;
+            cells_s(6 * i + 1) = m_;
+            cells_s(6 * i + 2) = f_;
+            cells_s(6 * i + 3) = s_;
+            cells_s(6 * i + 4) = xs_;
+            cells_s(6 * i + 5) = xr_;
+        }
+    }
+
+    //   //! NOTE: int_point must fit to cells_s.
+    double EvalReaction(const double V, const mfem::Vector &cells_s, int int_point, float t) const override
+    {
+        const float h = cells_s(6 * int_point + 0);
+        const float m = cells_s(6 * int_point + 1);
+        const float f = cells_s(6 * int_point + 2);
+        const float s = cells_s(6 * int_point + 3);
+        const float xs = cells_s(6 * int_point + 4);
+        const float xr = cells_s(6 * int_point + 5);
+
+        // Instantaneous gates
+        const float r = sigmoid(V, E_r, k_r, -1.0);
+        const float d = sigmoid(V, E_d, k_d, -1.0);
+        const float z = sigmoid(V, E_z, k_z, 1.0);
+        const float y = sigmoid(V, E_y, k_y, 1.0);
+            
+        // Currents
+        const float I_Na = g_Na * m * m * m * h * h * (V - E_Na);
+        const float I_K1 = g_K1 * z * (V - E_K);
+        const float I_to = g_to * r * s * (V - E_K);
+        const float I_CaL = g_CaL * d * f * (V - E_Ca);
+        const float I_Kr = g_Kr * xr * y * (V - E_K);
+        const float I_Ks = g_Ks * xs * (V - E_K);
+
+        return I_Na + I_K1 + I_to + I_CaL + I_Kr + I_Ks;
+    }
+
+    void InternalRHS(mfem::Vector&rhs, double V, mfem::Vector &cells_s) const override
+    {
+        // h
+        rhs(0) = rhs_h(V, cells_s(0));
+        // m
+        rhs(1) = rhs_m(V, cells_s(1));
+        // f
+        rhs(2) = rhs_f(V, cells_s(2));
+        // s
+        rhs(3) = rhs_s(V, cells_s(3));
+        // xs
+        rhs(4) = rhs_xs(V, cells_s(4));
+        // xr
+        rhs(5) = rhs_xr(V, cells_s(5));
+    }
+
+    void InternalEulerUpdate(const double V, mfem::Vector &cells_s, int int_point, float t, double Δt) const override
+    {
+        cells_s(6 * int_point + 0) += Δt*rhs_h (V, cells_s(6 * int_point + 0));
+        cells_s(6 * int_point + 1) += Δt*rhs_m (V, cells_s(6 * int_point + 1));
+        cells_s(6 * int_point + 2) += Δt*rhs_f (V, cells_s(6 * int_point + 2));
+        cells_s(6 * int_point + 3) += Δt*rhs_s (V, cells_s(6 * int_point + 3));
+        cells_s(6 * int_point + 4) += Δt*rhs_xs(V, cells_s(6 * int_point + 4));
+        cells_s(6 * int_point + 5) += Δt*rhs_xr(V, cells_s(6 * int_point + 5));
+
+        for(int i=0;i<6;i++) {
+            cells_s(6 * int_point + i) = std_clamp<float>(cells_s(6 * int_point + i), 0.0f, 1.0f);
+        }
+    }
+
+    int InternalDim() const override {return 6;}
+};

miniapps/spinner/inline-niederer-wedge.mesh

@@ -0,0 +1,9 @@
+MFEM INLINE mesh v1.0
+
+type = hex
+nx = 20
+ny = 7
+nz = 3
+sx = 20.0
+sy = 7.0
+sz = 3.0

miniapps/spinner/inline-test-quad.mesh

@@ -0,0 +1,7 @@
+MFEM INLINE mesh v1.0
+
+type = quad
+nx = 16
+ny = 16
+sx = 160.0
+sy = 160.0

miniapps/spinner/makefile

@@ -0,0 +1,71 @@
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/spinner/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SPINNER_SRC = spinner-lts.cpp
+SPINNER_OBJ = $(BLOCK_SOLVERS_SRC:.cpp=.o)
+
+SEQ_MINIAPPS = spinner-lts
+PAR_MINIAPPS = 
+
+ifeq ($(MFEM_USE_MPI),NO)
+	MINIAPPS = $(SEQ_MINIAPPS)
+else
+	MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
+all: $(MINIAPPS)
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+SPINNER_INCLUDES = $(if $(SRC),-I$(MFEM_DIR))
+
+%: %.o
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+spinner: $(SPINNER_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(SPINNER_OBJ) $(MFEM_LIBS)
+
+%.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS) $(SPINNER_OBJ)
+
+clean-exec:
+	@rm -rf mesh.* sol.* ParaView

miniapps/spinner/new-kelly-estimator.cxx

@@ -0,0 +1,279 @@
+template<typename Kernel>
+CustomKellyErrorEstimator<Kernel>::CustomKellyErrorEstimator(Kernel kernel_,
+                                         BilinearFormIntegrator& di_,
+                                         GridFunction& sol_,
+                                         FiniteElementSpace& flux_fespace_,
+                                         bool with_flux_,
+                                         const Array<int> &attributes_)
+   : kernel(kernel_)
+   , attributes(attributes_)
+   , flux_integrator(&di_)
+   , solution(&sol_)
+   , flux_space(&flux_fespace_)
+   , own_flux_fespace(false)
+   , with_flux(with_flux_)
+#ifdef MFEM_USE_MPI
+   , isParallel(dynamic_cast<ParFiniteElementSpace*>(sol_.FESpace()))
+#endif // MFEM_USE_MPI
+{
+   ResetCoefficientFunctions();
+}
+
+template<typename Kernel>
+CustomKellyErrorEstimator<Kernel>::CustomKellyErrorEstimator(Kernel kernel_,
+                                         BilinearFormIntegrator& di_,
+                                         GridFunction& sol_,
+                                         FiniteElementSpace* flux_fespace_,
+                                         bool with_flux_,
+                                         const Array<int> &attributes_)
+   : kernel(kernel_)
+   , attributes(attributes_)
+   , flux_integrator(&di_)
+   , solution(&sol_)
+   , flux_space(flux_fespace_)
+   , own_flux_fespace(true)
+   , with_flux(with_flux_)
+#ifdef MFEM_USE_MPI
+   , isParallel(dynamic_cast<ParFiniteElementSpace*>(sol_.FESpace()))
+#endif // MFEM_USE_MPI
+{
+   ResetCoefficientFunctions();
+}
+
+template<typename T>
+CustomKellyErrorEstimator<T>::~CustomKellyErrorEstimator()
+{
+   if (own_flux_fespace)
+   {
+      delete flux_space;
+   }
+}
+
+template<typename T>
+void CustomKellyErrorEstimator<T>::ResetCoefficientFunctions()
+{
+   compute_element_coefficient = [](Mesh* mesh, const int e)
+   {
+      return 1.0;
+   };
+
+   compute_face_coefficient = [](Mesh* mesh, const int f,
+                                 const bool shared_face)
+   {
+      auto FT = [&]()
+      {
+         #ifdef MFEM_USE_MPI
+            if (shared_face)
+            {
+               return dynamic_cast<ParMesh*>(mesh)->GetSharedFaceTransformations(f);
+            }
+         #endif // MFEM_USE_MPI
+         return mesh->GetFaceElementTransformations(f);
+      }();
+      const auto order = FT->GetFE()->GetOrder();
+
+      // Poor man's face diameter.
+      double diameter = 0.0;
+
+      Vector p1(mesh->SpaceDimension());
+      Vector p2(mesh->SpaceDimension());
+      // NOTE: We have no direct access to vertices for shared faces,
+      // so we fall back to compute the positions from the element.
+      // This can also be modified to compute the diameter for non-linear
+      // geometries by sampling along geometry-specific lines.
+      auto vtx_intrule = Geometries.GetVertices(FT->GetGeometryType());
+      const auto nip = vtx_intrule->GetNPoints();
+      for (int i = 0; i < nip; i++)
+      {
+         // Evaluate flux vector at integration point
+         auto fip1 = vtx_intrule->IntPoint(i);
+         FT->Transform(fip1, p1);
+
+         for (int j = i+1; j < nip; j++)
+         {
+            auto fip2 = vtx_intrule->IntPoint(j);
+            FT->Transform(fip2, p2);
+
+            diameter = std::max<double>(diameter, p2.DistanceTo(p1));
+         }
+      }
+      return diameter/(2.0*order);
+   };
+}
+
+template<typename Kernel>
+void CustomKellyErrorEstimator<Kernel>::ComputeEstimates()
+{
+   // Remarks:
+   // For some context you may have to consult the documentation of
+   // the FaceInfo class [1]. Also, the FaceElementTransformations
+   // documentation [2] may be helpful to grasp what is going on. Note
+   // that the FaceElementTransformations also works in the non-
+   // conforming case to transfer the Gauss points from the slave to
+   // the master element.
+   // [1]
+   // https://github.com/mfem/mfem/blob/02d0bfe9c18ce049c3c93a6a4208080fcfc96991/mesh/mesh.hpp#L94
+   // [2]
+   // https://github.com/mfem/mfem/blob/02d0bfe9c18ce049c3c93a6a4208080fcfc96991/fem/eltrans.hpp#L435
+
+   flux_space->Update(false);
+
+   auto xfes = solution->FESpace();
+   MFEM_ASSERT(xfes->GetVDim() == 1,
+               "Estimation for vector-valued problems not implemented yet.");
+   auto mesh = xfes->GetMesh();
+
+   this->error_estimates.SetSize(xfes->GetNE());
+   this->error_estimates = 0.0;
+
+   // 1. Compute fluxes in discontinuous space
+   GridFunction *flux =
+   #ifdef MFEM_USE_MPI
+      isParallel ? new ParGridFunction(dynamic_cast<ParFiniteElementSpace*>
+                                       (flux_space)) :
+   #endif // MFEM_USE_MPI
+   new GridFunction(flux_space);
+
+   *flux = 0.0;
+
+   // We pre-sort the array to speed up the search in the following loops.
+   if (attributes.Size())
+   {
+      attributes.Sort();
+   }
+
+   Array<int> xdofs, fdofs;
+   Vector el_x, el_f;
+   for (int e = 0; e < xfes->GetNE(); e++)
+   {
+      auto attr = xfes->GetAttribute(e);
+      if (attributes.Size() && attributes.FindSorted(attr) == -1)
+      {
+         continue;
+      }
+
+      xfes->GetElementVDofs(e, xdofs);
+      solution->GetSubVector(xdofs, el_x);
+
+      ElementTransformation* Transf = xfes->GetElementTransformation(e);
+      flux_integrator->ComputeElementFlux(*xfes->GetFE(e), *Transf, el_x,
+                                          *flux_space->GetFE(e), el_f, with_flux);
+
+      flux_space->GetElementVDofs(e, fdofs);
+      flux->AddElementVector(fdofs, el_f);
+   }
+
+   // 2. Add error contribution from local interior faces
+   for (int fi = 0; fi < mesh->GetNumFaces(); fi++)
+   {
+      if (mesh->FaceIsInterior(fi))
+      {
+         // Compute NC and face information
+         int FaceElement1, FaceElement2, NCFace;
+         mesh->GetFaceInfos(fi, &FaceElement1, &FaceElement2, &NCFace);
+         mesh->GetFaceElements(fi, &FaceElement1, &FaceElement2);
+
+         // We skip over master faces
+         bool isNCSlave    = FaceElement2 >= 0 && NCFace >= 0;
+         bool isConforming = FaceElement2 >= 0 && NCFace == -1;
+         if (isConforming || isNCSlave)
+         {
+            if (attributes.Size() &&
+                (attributes.FindSorted(mesh->GetAttribute(FaceElement1)) == -1
+                 || attributes.FindSorted(mesh->GetAttribute(FaceElement2)) == -1))
+            {
+               continue;
+            }
+
+            auto FT = mesh->GetFaceElementTransformations(fi);
+
+            const double kernel_value = kernel(FT, solution, flux);
+
+            // A local face is shared between two local elements, so we
+            // can get away with integrating the jump only once and add
+            // it to both elements. To minimize communication, the jump
+            // of shared faces is computed locally by each process.
+            auto h_k_face = compute_face_coefficient(mesh, fi, false);
+            error_estimates(FT->Elem1No) += h_k_face*kernel_value;
+            error_estimates(FT->Elem2No) += h_k_face*kernel_value;
+         }
+      }
+   }
+
+   current_sequence = solution->FESpace()->GetMesh()->GetSequence();
+
+#ifdef MFEM_USE_MPI
+   if (!isParallel)
+#endif // MFEM_USE_MPI
+   {
+      // Finalize element errors
+      for (int e = 0; e < xfes->GetNE(); e++)
+      {
+         auto factor = compute_element_coefficient(mesh, e);
+         // The sqrt belongs to the norm and hₑ to the indicator.
+         error_estimates(e) = sqrt(factor * error_estimates(e));
+      }
+
+      total_error = error_estimates.Norml2();
+      delete flux;
+      return;
+   }
+
+   #ifdef MFEM_USE_MPI
+
+   // 3. Add error contribution from shared interior faces
+   // Synchronize face data.
+
+   ParGridFunction *pflux = dynamic_cast<ParGridFunction*>(flux);
+   MFEM_VERIFY(pflux, "flux is not a ParGridFunction pointer");
+
+   ParMesh *pmesh = dynamic_cast<ParMesh*>(mesh);
+   MFEM_VERIFY(pmesh, "mesh is not a ParMesh pointer");
+
+   pflux->ExchangeFaceNbrData();
+
+   for (int sfi = 0; sfi < pmesh->GetNSharedFaces(); sfi++)
+   {
+      auto FT = pmesh->GetSharedFaceTransformations(sfi, true);
+
+      if (attributes.Size() &&
+          (attributes.FindSorted(FT->Elem1->Attribute) == -1
+           || attributes.FindSorted(FT->Elem2->Attribute) == -1))
+      {
+         continue;
+      }
+
+      // TODO Refactor this computation into a user-facing function.
+      // auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * xfes->GetFaceOrder(0)); // NOTE: This fails for DG
+      auto &int_rule = IntRules.Get(FT->FaceGeom,
+                                    2 * xfes->GetElementOrder(FT->Elem1No));
+      const auto nip = int_rule.GetNPoints();
+
+      const double kernel_value = kernel(FT, solution, flux);
+
+      auto h_k_face = compute_face_coefficient(mesh, sfi, true);
+      error_estimates(FT->Elem1No) += h_k_face*kernel_value;
+      // We skip "error_estimates(FT->Elem2No) += jump_integral"
+      // because the error is stored on the remote process and
+      // recomputed there.
+   }
+   delete flux;
+
+   // Finalize element errors
+   for (int e = 0; e < xfes->GetNE(); e++)
+   {
+      auto factor = compute_element_coefficient(mesh, e);
+      // The sqrt belongs to the norm and hₑ to the indicator.
+      error_estimates(e) = sqrt(factor * error_estimates(e));
+   }
+
+   // Finish by computing the global error.
+   auto pfes = dynamic_cast<ParFiniteElementSpace*>(xfes);
+   MFEM_VERIFY(pfes, "xfes is not a ParFiniteElementSpace pointer");
+
+   double process_local_error = pow(error_estimates.Norml2(),2.0);
+   MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
+                 MPI_SUM, pfes->GetComm());
+   total_error = sqrt(total_error);
+#endif // MFEM_USE_MPI
+}

miniapps/spinner/new-kelly-estimator.hpp

@@ -0,0 +1,289 @@
+//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+//!!!!!!!!!!!!!!!!!! TODO Planned to be merged into #3693 !!!!!!!!!!!!!!!!!!!!!
+//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+template<typename Kernel>
+class CustomKellyErrorEstimator final : public ErrorEstimator
+{
+public:
+   /// Function type to compute the local coefficient hₑ of an element.
+   using ElementCoefficientFunction =
+      std::function<double(Mesh*, const int)>;
+   /** @brief Function type to compute the local coefficient hₖ of a face. The
+       third argument is true for shared faces and false for local faces. */
+   using FaceCoefficientFunction =
+      std::function<double(Mesh*, const int, const bool)>;
+
+private:
+   Kernel kernel;
+
+   int current_sequence = -1;
+
+   Vector error_estimates;
+
+   double total_error = 0.0;
+
+   Array<int> attributes;
+
+   /** @brief A method to compute hₑ on per-element basis.
+
+       This method weights the error approximation on the element level.
+
+       Defaults to hₑ=1.0.
+   */
+   ElementCoefficientFunction compute_element_coefficient;
+
+   /** @brief A method to compute hₖ on per-face basis.
+
+       This method weights the error approximation on the face level. The
+       background here is that classical Kelly error estimator implementations
+       approximate the geometrical characteristic hₖ with the face diameter,
+       which should be also be a possibility in this implementation.
+
+       Defaults to hₖ=diameter/2p.
+   */
+   FaceCoefficientFunction compute_face_coefficient;
+
+   BilinearFormIntegrator* flux_integrator; ///< Not owned.
+   GridFunction* solution;               ///< Not owned.
+
+   FiniteElementSpace*
+   flux_space; /**< @brief Ownership based on own_flux_fes. */
+   bool own_flux_fespace; ///< Ownership flag for flux_space.
+
+   bool with_flux; ///< Use flux or gradient.
+   bool isParallel;
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = solution->FESpace()->GetMesh()->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence,
+                  "improper mesh update sequence");
+      return (mesh_sequence > current_sequence);
+   }
+
+   /** @brief Compute the element error estimates.
+
+       Algorithm outline:
+       1. Compute flux field for each element
+       2. Add error contribution from local interior faces
+       3. Add error contribution from shared interior faces
+       4. Finalize by computing hₖ and scale errors.
+   */
+   void ComputeEstimates();
+
+public:
+   /** @brief Construct a new CustomKellyErrorEstimator object for a scalar field.
+       @param di_         The bilinearform to compute the interface flux.
+       @param sol_        The solution field whose error is to be estimated.
+       @param flux_fes_   The finite element space for the interface flux.
+       @param attributes_ The attributes of the subdomain(s) for which the
+                          error should be estimated. An empty array results in
+                          estimating the error over the complete domain.
+   */
+   CustomKellyErrorEstimator(Kernel kernel_, BilinearFormIntegrator& di_, GridFunction& sol_,
+                       FiniteElementSpace& flux_fes_,
+                       bool with_flux,
+                       const Array<int> &attributes_ = Array<int>());
+
+   /** @brief Construct a new CustomKellyErrorEstimator object for a scalar field.
+       @param di_         The bilinearform to compute the interface flux.
+       @param sol_        The solution field whose error is to be estimated.
+       @param flux_fes_   The finite element space for the interface flux.
+       @param attributes_ The attributes of the subdomain(s) for which the
+                          error should be estimated. An empty array results in
+                          estimating the error over the complete domain.
+   */
+   CustomKellyErrorEstimator(Kernel kernel_, BilinearFormIntegrator& di_, GridFunction& sol_,
+                       FiniteElementSpace* flux_fes_,
+                       bool with_flux,
+                       const Array<int> &attributes_ = Array<int>());
+
+   ~CustomKellyErrorEstimator();
+
+   /// Get a Vector with all element errors.
+   const Vector& GetLocalErrors() override
+   {
+      if (MeshIsModified())
+      {
+         ComputeEstimates();
+      }
+      return error_estimates;
+   }
+
+   /// Reset the error estimator.
+   void Reset() override { current_sequence = -1; };
+
+   virtual double GetTotalError() const override { return total_error; }
+
+   /** @brief Change the method to compute hₑ on a per-element basis.
+       @param compute_element_coefficient_
+                        A function taking a mesh and an element index to
+                        compute the local hₑ for the element.
+   */
+   void SetElementCoefficientFunction(ElementCoefficientFunction
+                                      compute_element_coefficient_)
+   {
+      compute_element_coefficient = compute_element_coefficient_;
+   }
+
+   /** @brief Change the method to compute hₖ on a per-element basis.
+       @param compute_face_coefficient_
+                        A function taking a mesh and a face index to
+                        compute the local hₖ for the face.
+   */
+   void SetFaceCoefficientFunction(
+      FaceCoefficientFunction
+      compute_face_coefficient_)
+   {
+      compute_face_coefficient = compute_face_coefficient_;
+   }
+
+   /// Change the coefficients back to default as described above.
+   void ResetCoefficientFunctions();
+};
+
+#include "new-kelly-estimator.cxx" // Template impl
+
+struct DGFluxKernel {
+   // Eval buffers
+   IntegrationPoint ip;
+   FiniteElementSpace* flux_space;
+   Vector val;
+   Vector normal;
+   Vector ref_normal;
+
+   DGFluxKernel(FiniteElementSpace* flux_space_)
+   : flux_space(flux_space_)
+   , val(Vector(flux_space_->GetVDim()))
+   , normal(Vector(flux_space_->GetMesh()->SpaceDimension()))
+   , ref_normal(Vector(flux_space_->GetMesh()->Dimension()))
+   {
+      
+   }
+
+   double operator()(FaceElementTransformations *FT, GridFunction *solution, GridFunction *flux) {
+      auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * std::min(flux_space->GetElementOrder(FT->Elem1No), flux_space->GetElementOrder(FT->Elem2No)));
+      const auto nip = int_rule.GetNPoints();
+
+      double jump_integral = 0.0;
+
+      // Numerical integration of ∫ [[flux ⋅ n]] dS
+      for (int i = 0; i < nip; i++)
+      {
+         // Set up integration point
+         auto &fip = int_rule.IntPoint(i);
+         FT->Face->SetIntPoint(&fip);
+
+         // Compute normal - note that the normals match at each face integration point up to the sign!
+         if (flux_space->GetMesh()->Dimension() == flux_space->GetMesh()->SpaceDimension())
+         {
+            // This computes a weighted normal, so we divide by the weight to get back the normal.
+            CalcOrtho(FT->Face->Jacobian(), normal);
+            normal /= FT->Face->Weight();
+         }
+         else
+         {
+            // This computes a weighted normal, so we divide by the weight to get back the normal.
+            FT->Loc1.Transf.SetIntPoint(&fip);
+            FT->Loc1.Transform(fip, ip);
+            CalcOrtho(FT->Loc1.Transf.Jacobian(), ref_normal);
+            ref_normal /= FT->Face->Weight();
+            auto &e1 = FT->GetElement1Transformation();
+            e1.SetIntPoint(&ip);
+            e1.AdjugateJacobian().MultTranspose(ref_normal, normal);
+            // We have to cancel the additional weighting from the
+            // reference to spatial transformation in the line above.
+            normal /= e1.Weight();
+         }
+
+         // Evaluate flux jump at IP on element 1
+         FT->Loc1.Transf.SetIntPoint(&fip);
+         FT->Loc1.Transform(fip, ip);
+         flux->GetVectorValue(FT->Elem1No, ip, val);
+         double jump = val * normal;
+
+         // Evaluate flux jump at IP on element 2
+         FT->Loc2.Transf.SetIntPoint(&fip);
+         FT->Loc2.Transform(fip, ip);
+         flux->GetVectorValue(FT->Elem2No, ip, val);
+         jump -= val * normal;
+
+         // Finalize integral
+         jump_integral += jump*jump*fip.weight * FT->Face->Weight();
+      }
+
+      return jump_integral;
+   }
+};
+
+
+struct DGWeightedFluxKernel {
+   // Eval buffers
+   IntegrationPoint ip;
+   FiniteElementSpace* flux_space;
+   Vector val;
+   Vector normal;
+   Vector ref_normal;
+
+   DGWeightedFluxKernel(FiniteElementSpace* flux_space_)
+   : flux_space(flux_space_)
+   , val(Vector(flux_space_->GetVDim()))
+   , normal(Vector(flux_space_->GetMesh()->SpaceDimension()))
+   , ref_normal(Vector(flux_space_->GetMesh()->Dimension()))
+   {
+   }
+
+   double operator()(FaceElementTransformations *FT, GridFunction *solution, GridFunction *flux) {
+      auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * std::min(flux_space->GetElementOrder(FT->Elem1No), flux_space->GetElementOrder(FT->Elem2No)));
+      const auto nip = int_rule.GetNPoints();
+
+      double jump_integral = 0.0;
+
+      // Numerical integration of ∫ [[flux ⋅ n]] dS
+      for (int i = 0; i < nip; i++)
+      {
+         // Set up integration point
+         auto &fip = int_rule.IntPoint(i);
+         FT->Face->SetIntPoint(&fip);
+
+         // Compute normal - note that the normals match at each face integration point up to the sign!
+         if (flux_space->GetMesh()->Dimension() == flux_space->GetMesh()->SpaceDimension())
+         {
+            // This computes a weighted normal, so we divide by the weight to get back the normal.
+            CalcOrtho(FT->Face->Jacobian(), normal);
+         }
+         else
+         {
+            // This computes a weighted normal, so we divide by the weight to get back the normal.
+            FT->Loc1.Transf.SetIntPoint(&fip);
+            FT->Loc1.Transform(fip, ip);
+            CalcOrtho(FT->Loc1.Transf.Jacobian(), ref_normal);
+            auto &e1 = FT->GetElement1Transformation();
+            e1.SetIntPoint(&ip);
+            e1.AdjugateJacobian().MultTranspose(ref_normal, normal);
+            // We have to cancel the additional weighting from the
+            // reference to spatial transformation in the line above.
+            normal /= e1.Weight();
+         }
+
+         // Evaluate flux jump at IP on element 1
+         FT->Loc1.Transf.SetIntPoint(&fip);
+         FT->Loc1.Transform(fip, ip);
+         flux->GetVectorValue(FT->Elem1No, ip, val);
+         double jump = val * normal;
+
+         // Evaluate flux jump at IP on element 2
+         FT->Loc2.Transf.SetIntPoint(&fip);
+         FT->Loc2.Transform(fip, ip);
+         flux->GetVectorValue(FT->Elem2No, ip, val);
+         jump -= val * normal;
+
+         // Finalize integral
+         jump_integral += jump*jump*fip.weight * FT->Face->Weight();
+      }
+
+      return jump_integral;
+   }
+};

miniapps/spinner/utils.hpp

@@ -0,0 +1,390 @@
+// Buffers to store all the necessary
+// information to advance an element in time
+struct EvalCache {
+  // Buffers for the element action
+  mfem::Vector phimelvals;
+  mfem::Vector phimelvals2;
+  mfem::Vector phimelvalsfacerhs;
+  mfem::Vector phimelvals_rhs;
+  mfem::Vector phimelvals_rhs2;
+  mfem::Vector selvals;
+  mfem::Vector selvals2;
+
+  mfem::DenseMatrix a_el;
+  mfem::DenseMatrix a2_el;
+
+  mfem::DenseMatrix minv_el;
+  mfem::DenseMatrix f_el;
+  mfem::Array<int> e_faces;
+  mfem::Array<int> f_ori;
+
+  // Buffers for the face action
+  mfem::Array<int> facedofblocks;
+  mfem::BlockVector phimrhsfacevals; // Block vector for face action (element pairs)
+  mfem::Vector phimrhsfacevals_e1view;
+  mfem::Vector phimrhsfacevals_e2view;
+  mfem::Vector phimrhsfacevals_view_prev; // Last time step value
+  mfem::Vector phimrhsfacevals_view_pred; // Predictor value
+  mfem::BlockVector phimrhsfacevals_buf;
+  mfem::BlockVector phimrhsfacevals_buf2;
+
+  // Dof buffers
+  mfem::Array<int> phimdofs;
+  mfem::Array<int> sdofs;
+  mfem::Array<int> vdofs1, vdofs2;
+
+  // ndofs = number of field dofs per element
+  // idim  = dimension of the internal problem
+  // nfaces = number of faces per elements
+  EvalCache(int ndofs, int idim, int nfaces) {
+    phimelvals        = mfem::Vector(ndofs);
+    phimelvals2       = mfem::Vector(ndofs);
+    phimelvalsfacerhs = mfem::Vector(ndofs);
+    phimelvals_rhs    = mfem::Vector(ndofs);
+    phimelvals_rhs2   = mfem::Vector(ndofs);
+
+    e_faces = mfem::Array<int>(nfaces);
+    f_ori   = mfem::Array<int>(nfaces);
+    
+    a_el  = mfem::DenseMatrix(ndofs);
+    a2_el = mfem::DenseMatrix(ndofs);
+
+    facedofblocks.Append(0);
+    facedofblocks.Append(ndofs);
+    facedofblocks.Append(2 * ndofs);
+
+    phimrhsfacevals.Update(facedofblocks);
+    phimrhsfacevals_buf.Update(facedofblocks);
+    phimrhsfacevals_buf2.Update(facedofblocks);
+
+    phimrhsfacevals_e1view    = mfem::Vector(ndofs);
+    phimrhsfacevals_e2view    = mfem::Vector(ndofs);
+    phimrhsfacevals_view_prev = mfem::Vector(ndofs);
+    phimrhsfacevals_view_pred = mfem::Vector(ndofs);
+
+    phimdofs   = mfem::Array<int>(ndofs);
+    sdofs    = mfem::Array<int>(idim * ndofs);
+    selvals  = mfem::Vector(idim * ndofs);
+    selvals2 = mfem::Vector(idim * ndofs);
+  }
+};
+
+// AMR update for fe spaces, gridfunctions and the forms
+inline void Update(std::vector<mfem::FiniteElementSpace *> &fespaces, std::vector<mfem::GridFunction *> &xs,
+            mfem::BilinearForm &a, mfem::BilinearForm &m) {
+  MFEM_PERF_FUNCTION;
+
+  // Update the space: recalculate the number of DOFs and construct a matrix
+  // that will adjust any GridFunctions to the new mesh state.
+  MFEM_PERF_BEGIN("FESpaces");
+  for (auto fespace : fespaces) {
+    fespace->Update();
+  }
+  MFEM_PERF_END("FESpaces");
+
+  // Interpolate the solution on the new mesh by applying the transformation
+  // matrix computed in the finite element space. Multiple GridFunctions could
+  // be updated here.
+  MFEM_PERF_BEGIN("GridFunctions");
+  for (auto x : xs) {
+    x->Update();
+  }
+  MFEM_PERF_END("GridFunctions");
+
+  // Inform the linear and bilinear forms that the space has changed.
+  MFEM_PERF_BEGIN("Forms");
+  a.Update();
+  m.Update();
+  MFEM_PERF_END("Forms");
+
+  // Free any transformation matrices to save memory.
+  MFEM_PERF_BEGIN("FESpaces");
+  for (auto fespace : fespaces) {
+    fespace->UpdatesFinished();
+  }
+  MFEM_PERF_END("FESpaces");
+}
+
+// Could not find this function in MFEM
+inline void GetFaceIndices(mfem::Array<int>&faces, mfem::Array<int>&orientations, mfem::Mesh* mesh, int ei) {
+  if(mesh->Dimension() == 1) {
+    mesh->GetElementVertices(ei, faces);
+  } else if(mesh->Dimension() == 2) {
+    mesh->GetElementEdges(ei, faces, orientations);
+  } else {
+    mesh->GetElementFaces(ei, faces, orientations);
+  }
+}
+
+
+// Linearize vector dofs to simplify the iteration order
+inline void GetElementVDofsLinear(const mfem::FiniteElementSpace &fespace, const int element, mfem::Array<int> &dofs)
+{
+  const int vdim = fespace.GetVDim();
+  fespace.GetElementDofs(element, dofs);
+  const int size = dofs.Size();
+  dofs.SetSize(size*vdim);
+  // Stride dofs first
+  for(int node = size-1; node >= 0; node--)
+  {
+    dofs[node*vdim] = dofs[node];
+  }
+  for(int node = 0; node < size; node++)
+  {
+    for (int vd = 1; vd < vdim; vd++)
+    {
+      dofs[node*vdim+vd] = -1;
+    }
+  }
+  // Fill gaps
+  for (int node = 0; node < size; node++)
+  {
+    for (int vd = 1; vd < vdim; vd++)
+    {
+      dofs[node*vdim+vd] = fespace.DofToVDof(dofs[node*vdim], vd);
+    }
+    dofs[node*vdim] = fespace.DofToVDof(dofs[node*vdim], 0);
+  }
+}
+
+// Helper for the initial condition
+template<class GF>
+void FibrillationInit(GF& phimgf, GF& sgf, int h_internal_offset) {
+  auto phimfespace = phimgf.FESpace();
+  auto sfespace  = sgf.FESpace();
+  auto mesh = phimfespace->GetMesh();
+  const int sdim = mesh->SpaceDimension();
+
+  mfem::Vector mesh_min(sdim), mesh_max(sdim);
+//   auto pmesh = dynamic_cast<mfem::ParMesh*>(mesh);
+//   if(pmesh)
+//   {
+//     pmesh->GetBoundingBox(mesh_min, mesh_max, 0);
+//   }
+//   else 
+  {
+    mesh->GetBoundingBox(mesh_min, mesh_max, 0);
+  }
+
+  int el = -1;
+  mfem::ElementTransformation *T = NULL;
+  const mfem::FiniteElement *fe = NULL;
+
+  mfem::Vector pos(sdim);
+  phimfespace->BuildDofToArrays();
+  for (int dof = 0; dof < phimfespace->GetNDofs(); dof++)
+  {
+      int j = phimfespace->GetElementForDof(dof);
+      if (el != j)
+      {
+          el = j;
+          T = phimfespace->GetElementTransformation(el);
+          fe = phimfespace->GetFE(el);
+      }
+      int ld = phimfespace->GetLocalDofForDof(dof);
+      const mfem::IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
+      T->SetIntPoint(&ip);
+      T->Transform(ip, pos);
+
+      const double width = mesh_max(0)-mesh_min(0);
+      const double offset = (pos(0)-mesh_min(0))/width;
+      phimgf(dof) = (1.0-offset);
+  }
+
+  if(sdim == 1) return;
+
+  sfespace->BuildDofToArrays();
+  for (int dof = 0; dof < sfespace->GetNDofs(); dof++)
+  {
+      int j = sfespace->GetElementForDof(dof);
+      if (el != j)
+      {
+          el = j;
+          T = sfespace->GetElementTransformation(el);
+          fe = sfespace->GetFE(el);
+      }
+      int vdof = sfespace->DofToVDof(dof, h_internal_offset);
+      int ld = sfespace->GetLocalDofForDof(dof);
+      const mfem::IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
+      T->SetIntPoint(&ip);
+      T->Transform(ip, pos);
+
+      const double width = mesh_max(1)-mesh_min(1);
+      const double offset = (pos(1)-mesh_min(1))/width;
+      sgf(vdof) = 0.1*offset + (1.0-offset)*0.6;
+  }
+}
+
+// Compute the element Jacobian for an SIPG diffusion element in ODE form
+inline void ComputeElementDiffusion(mfem::DenseMatrix& element_jac, mfem::DenseMatrix& element_jac_massinv, mfem::Mesh* mesh, int ei, const mfem::NCMesh::NCList& ncfacelist, mfem::BilinearForm& minv, mfem::BilinearForm& a, const std::vector<mfem::DenseMatrix>& face_matrices, EvalCache& ec)
+{
+  GetFaceIndices(ec.e_faces, ec.f_ori, mesh, ei);
+
+  element_jac = 0.0;
+  a.ComputeElementMatrix(ei, element_jac);
+
+  for (int fi : ec.e_faces) {
+    int fe1, fe2;
+    mesh->GetFaceElements(fi,&fe1,&fe2);
+    DEBUG_PRINT("    Visiting face " << fi << " with elements " << fe1 << " " << fe2);
+    if (fe2 >= 0) { // Not on boundary
+      MFEM_PERF_FINE_BEGIN("EvaluateFaceCall");
+      //AddFace
+      const auto& face_mat = face_matrices[fi];
+      const auto ndofs = ec.facedofblocks[1];
+      if(fe1 == ei) {
+        for(int i=0;i<ndofs;i++) {
+            for(int j=0;j<ndofs;j++) {
+                element_jac(i,j) += face_mat(i,j);
+            }
+        }
+      } else {
+        for(int i=0;i<ndofs;i++) {
+            for(int j=0;j<ndofs;j++) {
+                element_jac(i,j) += face_mat(i+ndofs,j+ndofs);
+            }
+        }
+      }
+      MFEM_PERF_FINE_END("EvaluateFaceCall");
+    } else {
+      int Inf1, Inf2, NCFace;
+      mesh->GetFaceInfos(fi, &Inf1, &Inf2, &NCFace);
+      if(NCFace < 0) {
+        // No-flux
+        // const auto& face_mat = face_matrices[fi];
+        // const auto ndofs = ec.facedofblocks[1];
+        // for(int i=0;i<ndofs;i++) {
+        //     for(int j=0;j<ndofs;j++) {
+        //         element_jac(i,j) += face_mat(i,j);
+        //     }
+        // }
+        continue; // Not a master face (i.e. a boundary face)
+      }
+      // MFEM_PERF_FINE_BEGIN("NCInfo");
+      DEBUG_PRINT("      NCFace=" << NCFace);
+      auto& masterinfo = ncfacelist.masters[NCFace];
+      DEBUG_PRINT("      master index=" << masterinfo.index << " local=" << int(masterinfo.local));
+      DEBUG_PRINT("      slave range=" << masterinfo.slaves_begin << ":" << masterinfo.slaves_end);
+      // MFEM_PERF_FINE_END("NCInfo");
+      for(int slave = masterinfo.slaves_begin; slave < masterinfo.slaves_end; slave++) {
+        auto& slaveinfo = ncfacelist.slaves[slave];
+        DEBUG_PRINT("        slave index=" << slaveinfo.index << " local=" << int(slaveinfo.local));
+        if(slaveinfo.index < 0) { // Degenerate face-edge constraint
+          continue;
+        }
+        mesh->GetFaceElements(slaveinfo.index, &fe1, &fe2);
+        const auto& face_mat = face_matrices[slaveinfo.index];
+        const auto ndofs1 = ec.facedofblocks[1];
+        const auto ndofs2 = ec.facedofblocks[2]-ec.facedofblocks[1];
+        MFEM_PERF_FINE_BEGIN("EvaluateFaceCall2");
+        if(ei == fe1) {
+          for(int i=0;i<ndofs1;i++) {
+            for(int j=0;j<ndofs1; j++) {
+              element_jac(i,j) += face_mat(i,j);
+            }
+          }
+        } else if(ei == fe2) {
+          for(int i=0;i<ndofs2;i++) {
+            for(int j=0;j<ndofs2; j++) {
+              element_jac(i,j) += face_mat(i+ndofs1,j+ndofs1);
+            }
+          }
+        } else {
+          std::cout << "Face-element table corrupted? (Case 1) ei=" << ei << " fi=" << fi << std::endl;
+          std::exit(-1);
+        }
+        MFEM_PERF_FINE_END("EvaluateFaceCall2");
+      }
+    }
+  }
+
+  minv.ComputeElementMatrix(ei, ec.minv_el); // M_e^-1
+  mfem::Mult(ec.minv_el, element_jac, element_jac_massinv); // J_e = M_e^-1 * K_e
+}
+
+// Wrapper to make clear what we intend do to
+inline void ComputeElementJacobian(mfem::DenseMatrix& element_jac, mfem::DenseMatrix& element_jac_massinv, mfem::Mesh* mesh, int ei, const mfem::NCMesh::NCList& ncfacelist, mfem::BilinearForm& minv, mfem::BilinearForm& a, std::vector<mfem::DenseMatrix>& face_matrices, std::shared_ptr<const AbstractCellSolver> stepper, EvalCache& ec)
+{
+  ComputeElementDiffusion(element_jac, element_jac_massinv, mesh, ei, ncfacelist, minv, a, face_matrices, ec);
+}
+
+// Dimension independent getter for faces
+inline const mfem::NCMesh::NCList& GetNCFaceList(const mfem::Mesh* mesh) {
+  static mfem::NCMesh::NCList emptylist;
+  if(mesh->ncmesh == nullptr || mesh->Dimension() == 1) {
+    return emptylist;
+  }
+  return mesh->Dimension() == 2 ? mesh->ncmesh->GetEdgeList() : mesh->ncmesh->GetFaceList();
+}
+
+// TODO try to integrate this with some kind of EABilinearForm
+inline std::vector<mfem::DenseMatrix> ComputeFaceMatrices(mfem::Mesh* mesh, mfem::FiniteElementSpace* fespace, mfem::BilinearFormIntegrator* dginteg, mfem::BilinearFormIntegrator* dgbdrinteg)
+{
+    std::vector<mfem::DenseMatrix> face_matrices(mesh->GetNumFaces());
+    mfem::DenseMatrix f_el;
+    for (int fi = 0; fi < mesh->GetNumFaces(); fi++) {
+        face_matrices[fi] = 0.0;
+        int fe1, fe2;
+        mesh->GetFaceElements(fi,&fe1,&fe2);
+        if (fe2 < 0) { // Maybe on boundary
+            // Logic taken from https://github.com/mfem/mfem/blob/541f10f9b44fd52bacde793a08dc19e53a8972ec/fem/bilinearform.cpp#L615-L628
+            // const auto fbi = face_to_be[fi];
+            // if(fbi < 0) { // Filter ncface
+            //     continue;
+            // }
+            // auto tr = mesh->GetBdrFaceTransformations(fbi);
+            // if (tr == nullptr) { 
+            //     continue;
+            // }
+            // dgbdrinteg->AssembleFaceMatrix(*fespace->GetFE(tr->Elem1No),
+            //                             *fespace->GetFE(tr->Elem1No), *tr, f_el);
+            // face_matrices[fi] = f_el;
+            continue;
+        }
+        auto tr = mesh->GetInteriorFaceTransformations(fi);
+        if (tr == nullptr) {
+            continue;
+        }
+        dginteg->AssembleFaceMatrix(*fespace->GetFE(tr->Elem1No),
+                                    *fespace->GetFE(tr->Elem2No), *tr, f_el);
+        face_matrices[fi] = f_el;
+    }
+    return face_matrices;
+}
+
+inline void EigenvaluesGershgorin(const mfem::DenseMatrix& A, mfem::Vector& evs)
+{
+  for(int row=0;row<A.Height();row++) {
+    evs(row) = A(row,row);
+    for(int col=0;col<row;col++) {
+      evs(row) += std::abs(A(row,col));
+    }
+    for(int col=row+1;col<A.Width();col++) {
+      evs(row) += std::abs(A(row,col));
+    }
+  }
+}
+
+struct Stimulus {
+  double dist_max = 0.0;
+  double t_max = 0.0;
+  double stim_max = 0.0;
+  double capacitance = 1.0;
+};
+
+void PrintBanner(std::ostream &out) {
+  out << "+-----------------------------------------------------------------+"<< std::endl;
+  out << "                                                                  "<< std::endl;
+  out << "        __     __       _    _      ____       @@@@       @@@@    "<< std::endl;
+  out << "       /..|   /..|     /./  /./   /. ___/    @@@@@@@@   @@@@@@@@  "<< std::endl;
+  out << "      /. .|  /. .|    /./__/./   /. /__     @@@@@@@@@@@@@@@@@@@@@ "<< std::endl;
+  out << "     /./|.| /./|.|   /.____./   /_... /     @@@@@@@@@@@@@@@@@@@@@ "<< std::endl;
+  out << "    /./ |.|/./ |.|  /./  /./   ___/. /       @@@@@@@@@@@@@@@@@@@  "<< std::endl;
+  out << "   /_/  |___/  |_| /_/  /_/   /_____/          @@@@@@@@@@@@@@@    "<< std::endl;
+  out << "                                                 @@@@@@@@@@@      "<< std::endl;
+  out << "       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                @@@@@         "<< std::endl;
+  out << "       ~~~~ Spinner LTS Miniapp ~~~~                  @           "<< std::endl;
+  out << "       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                              "<< std::endl;
+  out << "                                                                  "<< std::endl;
+  out << "+-----------------------------------------------------------------+"<< std::endl;
+}