ScalarEq: Implement modal Bernstein DG evolution...

ScalarEq: Implement modal Bernstein DG evolution (!207)

MWE to check consistency of operators.

This is not intended to work with any HRSC schemes or 2d etc.

src/ScalarEq/ScalarEq.jl

@@ -3,6 +3,7 @@ module ScalarEq
 
 
 using dg1d
+using LinearAlgebra
 using LoopVectorization
 
 

src/ScalarEq/rhs.jl

@@ -267,7 +267,18 @@ end
 #######################################################################
 
 
-function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+function rhs!(env, mesh::Mesh1d{DGElement}, P::Project, hrsc::Nothing)
+  if mesh.element.kind === :nodal_lagrange
+    nodal_lagrange_rhs!(env, mesh, P, hrsc)
+  elseif mesh.element.kind === :modal_bernstein
+    modal_bernstein_rhs!(env, mesh, P, hrsc)
+  else
+    TODO(mesh.element.kind)
+  end
+end
+
+
+function nodal_lagrange_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
 
   @unpack cache, mesh = env
   @unpack equation = P
@@ -288,6 +299,34 @@ function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
 end
 
 
+function modal_bernstein_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+
+  @unpack cache, mesh = env
+  @unpack equation = P
+
+  @unpack u_modal, flx_u, src_u = get_static_variables(cache)
+  @unpack u                     = get_dynamic_variables(cache)
+  @unpack rhs_u                 = get_rhs_variables(cache)
+  @unpack nflx_u, bdry_u        = get_bdry_variables(cache)
+
+  u_modal .= u
+  for (v_u, v_u_modal) in zip(eachcell(mesh,u),eachcell(mesh,u_modal))
+    mul!(v_u, P.prms.V_bernstein, v_u_modal)
+  end
+
+  dg1d.interpolate_face_data!(mesh, u, bdry_u)
+  broadcast_volume!(flux_source, equation, mesh)
+  broadcast_faces!(llf_1d, equation, mesh)
+  # broadcast_bdry!(bdryllf_1d, equation, P.bdrycond, mesh)
+
+  compute_rhs_weak_form!(rhs_u, flx_u, src_u, nflx_u, mesh)
+
+  u .= u_modal
+
+  return
+end
+
+
 function rhs!(env, mesh::Mesh2d, P::Project, hrsc::Nothing)
 
   @unpack cache = env

src/ScalarEq/setup.jl

@@ -13,9 +13,11 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   bernstein = BernsteinReconstruction(mesh)
   analyze_error = prms["ScalarEq"]["analyze_error"]
   analyze_error_norm = prms["ScalarEq"]["analyze_error_norm"]
+  V_bernstein = dg1d.Bernstein.vandermonde_matrix(mesh.element.z)
   fixedprms = (; bernstein, av_derivative_scheme, av_drag, hrsc, muscl_omega,
                  slope_limiter_method, slope_limiter_tvb_M,
                  av_recompute_substeps,
+                 V_bernstein,
                  analyze_error, analyze_error_norm)
 
   # construct Project
@@ -123,6 +125,10 @@ function register_equation!(mesh::Mesh1d, P::Project)
     else
       TODO(P.prms.hrsc)
     end
+  elseif mesh.element isa DGElement
+    if mesh.element.kind === :modal_bernstein
+      register_variables!(mesh.cache, static_variablenames=(:u_modal,))
+    end
   end
 end
 function register_equation!(mesh::Mesh2d, P::Project)

src/dgelement.jl

@@ -6,7 +6,7 @@ struct DGElement
   # quadrature rule
   # quadr_w are such that data sampled on `z` can be directly passed to the integrate
   # functions below. In particular, any interpolation from `z` onto `quadr_z`
-  # (need kind === :nodal_bernstein) is already accounted for.
+  # (need kind === :modal_bernstein) is already accounted for.
   quadr::Symbol             # :lgl, :lgl_mass_lumped, :glgl
   quadr_z::Vector{Float64}  # quadrature points
   quadr_w::Vector{Float64}  # quadrature weights
@@ -61,7 +61,7 @@ struct DGElement
       l_rhs = zeros(Float64, N+1)
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
-    elseif kind === :nodal_bernstein
+    elseif kind === :modal_bernstein
       if quadr === :lgl || quadr === :lgl_lumped
         quadr_z, quadr_w = LGL.rule(Npts)
       elseif quadr === :glgl
@@ -82,7 +82,7 @@ struct DGElement
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
     else
-      error("Invalid quadrature rule specified")
+      error("Invalid DG formulation specified: $kind. Use on of nodal_lagrange, modal_bernstein.")
     end
 
     MDM    = invM * transpose(S)

src/parameters.jl

@@ -510,12 +510,14 @@ end
   - `nodal_lagrange` ... Nodal approximation using Lagrange polynomials collocated with
     the quadrature nodes provided by `basis`. This is the standard method discussed
     in Hesthaven, Warburton, 2007.
-  - `nodal_bernstein` ... (Quasi-)Nodal approximation using Bernstein polynomials.
-    The coefficients correspond to the nodal values sampled on an equidistant grid.
+  - `modal_bernstein` ... A modal approximation using Bernstein polynomials.
+    The coefficients correspond to the nodal values sampled on an equidistant grid,
+    but the Bernstein polynomials are not interpolating.
+    Hence, this Ansatz is also refered to as quasi-nodal.
     See Beisiegel, Behrens, Environ Earth Sci (2015) 74:7275–7284.
   """
   dg_kind = "nodal_lagrange"
-  @check scheme in [ "nodal_lagrange", "nodal_bernstein" ]
+  @check scheme in [ "nodal_lagrange", "modal_bernstein" ]
 
   """
   Periodicity in cartesian directions `x,y`

test/IntegrationTests/refs/modal_bernstein_advection_sine/modal_bernstein_advection_sine.toml

+[Evolution]
+cfl = 0.5
+tend = 2.2
+
+[ScalarEq]
+analyze_error = false
+equation = "advection"
+id = "sine"
+
+[Output]
+substep_every_iteration = 200
+interpolate_variables = ["u"]
+substep_variables = ["u"]
+interpolate_aligned_ts = [0.3, 0.6, 0.8]
+aligned_ts = [0.05, 0.1, 0.15, 0.2]
+interpolate_nodes = "$(collect(range(0.0,1.0,length=10)))"
+variables1d = ["u"]
+
+[Mesh]
+dg_kind = "modal_bernstein"
+range = [0.0, 1.0]
+k = 40
+basis = "lgl"
+n = 4

test/IntegrationTests/refs/testconfig.toml

@@ -7,6 +7,9 @@ default_abstol1d = 1e-8
 [advection_sine]
 variables1d = [ "u", "u_err" ]
 
+[modal_bernstein_advection_sine]
+variables1d = [ "u" ]
+
 [advection_sine_analysis]
 variables0d = [ "u_err_norm", "u_norm2" ]
 variables1d = [ "u" ]

test/UnitTests/src/test_dgelement.jl

@@ -209,7 +209,7 @@ end
     @test all(MD.≈S)
   end
 
-  el = dg1d.DGElement(N, :lgl, :nodal_bernstein)
+  el = dg1d.DGElement(N, :lgl, :modal_bernstein)
   T, invT = dg1d.Bernstein.bernstein_legendre_transformation(el.Npts)
   M_legendre = dg1d.Legendre.mass_matrix(el.Npts)
   TMT = T*M_legendre*transpose(T) |> dg1d.purge_zeros