Taal: first support for MHD

examples/2d/parameters_ec_mhd.jl

@@ -0,0 +1,55 @@
+# TODO: Taal refactor, rename to
+# - mhd_ec.jl
+# or something similar?
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = IdealGlmMhdEquations2D(1.4)
+
+initial_conditions = Trixi.initial_conditions_ec_test
+
+surface_flux = flux_derigs_etal
+volume_flux  = flux_derigs_etal
+solver = DGSEM(3, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2, -2)
+coordinates_max = ( 2,  2)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=10_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_conditions, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+save_solution = SaveSolutionCallback(interval=10,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+
+callbacks = CallbackSet(summary_callback, stepsize_callback, analysis_callback, save_solution, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

src/callbacks/analysis.jl

@@ -370,6 +370,12 @@ pretty_form_file(::typeof(energy_kinetic)) = "e_kinetic"
 pretty_form_repl(::typeof(energy_internal)) = "∑e_internal"
 pretty_form_file(::typeof(energy_internal)) = "e_internal"
 
+pretty_form_repl(::Val{:l2_divb}) = "L2 ∇⋅B"
+pretty_form_file(::Val{:l2_divb}) = "l2_divb"
+
+pretty_form_repl(::Val{:linf_divb}) = "L∞ ∇⋅B"
+pretty_form_file(::Val{:linf_divb}) = "linf_divb"
+
 
 # specialized implementations specific to some solvers
 include("analysis_dg2d.jl")

src/callbacks/analysis_dg2d.jl

@@ -89,8 +89,37 @@ function analyze(::typeof(entropy_timederivative), du::AbstractArray{<:Any,4}, u
   return dsdu_ut
 end
 
-# TODO: Taal implement
-# function analyze(::Val{:l2_divb}, du::AbstractArray{<:Any,4}, u, t,
-#   mesh::TreeMesh{2}, equations, dg::DG, cache)
-# function analyze(::Val{:linf_divb}, du::AbstractArray{<:Any,4}, u, t,
-#     mesh::TreeMesh{2}, equations, dg::DG, cache)
+
+
+function analyze(::Val{:l2_divb}, du::AbstractArray{<:Any,4}, u, t,
+                 mesh::TreeMesh{2}, equations, dg::DG, cache)
+  integrate(mesh, equations, dg, cache, u, cache, dg.basis.derivative_matrix) do u, i, j, element, equations, dg, cache, derivative_matrix
+      divb = zero(eltype(u))
+      for k in eachnode(dg)
+        divb += ( derivative_matrix[i, k] * u[6, k, j, element] +
+                  derivative_matrix[j, k] * u[7, i, k, element] )
+      end
+      divb *= cache.elements.inverse_jacobian[element]
+  end
+end
+
+function analyze(::Val{:linf_divb}, du::AbstractArray{<:Any,4}, u, t,
+                 mesh::TreeMesh{2}, equations, dg::DG, cache)
+  @unpack derivative_matrix, weights = dg.basis
+
+  # integrate over all elements to get the divergence-free condition errors
+  linf_divb = zero(eltype(u))
+  for element in eachelement(dg, cache)
+    for j in eachnode(dg), i in eachnode(dg)
+      divb = 0.0
+      for k in eachnode(dg)
+        divb += ( derivative_matrix[i, k] * u[6, k, j, element] +
+                  derivative_matrix[j, k] * u[7, i, k, element] )
+      end
+      divb *= cache.elements.inverse_jacobian[element]
+      linf_divb = max(linf_divb, abs(divb))
+    end
+  end
+
+  return linf_divb
+end

src/equations/2d/ideal_glm_mhd.jl

@@ -218,6 +218,39 @@ end
   end
 end
 
+@inline function calcflux_twopoint_nonconservative!(f1, f2, u::AbstractArray{<:Any,4}, element,
+                                                    equations::IdealGlmMhdEquations2D, dg, cache)
+  for j in eachnode(dg), i in eachnode(dg)
+    rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = get_node_vars(u, equations, dg, i, j, element)
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v3 = rho_v3 / rho
+
+    # Powell nonconservative term: Φ^Pow = (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+    phi_pow = 0.5 * SVector(0, B1, B2, B3, v1*B1 + v2*B2 + v3*B3, v1, v2, v3, 0)
+
+    # Galilean nonconservative term: Φ^Gal_{1,2} = (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+    # x-direction
+    phi_gal_x = 0.5 * SVector(0, 0, 0, 0, v1*psi, 0, 0, 0, v1)
+    # y-direction
+    phi_gal_y = 0.5 * SVector(0, 0, 0, 0, v2*psi, 0, 0, 0, v2)
+
+    # add both nonconservative terms into the volume
+    for l in eachnode(dg)
+      _, _, _, _, _, B1, _, _, psi = get_node_vars(u, equations, dg, l, j, element)
+      for v in eachvariable(equations)
+        f1[v, l, i, j] += phi_pow[v] * B1 + phi_gal_x[v] * psi
+      end
+      _, _, _, _, _, _, B2, _, psi = get_node_vars(u, equations, dg, i, l, element)
+      for v in eachvariable(equations)
+        f2[v, l, i, j] += phi_pow[v] * B2 + phi_gal_y[v] * psi
+      end
+    end
+  end
+
+  return nothing
+end
+
 
 """
     flux_derigs_etal(u_ll, u_rr, orientation, equation::IdealGlmMhdEquations2D)
@@ -491,8 +524,8 @@ end
   v1 = rho_v1 / rho
   v2 = rho_v2 / rho
   v3 = rho_v3 / rho
-  cf_x_direction = calc_fast_wavespeed(u_node, 1, equation)
-  cf_y_direction = calc_fast_wavespeed(u_node, 2, equation)
+  cf_x_direction = calc_fast_wavespeed(u, 1, equation)
+  cf_y_direction = calc_fast_wavespeed(u, 2, equation)
   cf_max = max(cf_x_direction, cf_y_direction)
   equation.c_h = max(equation.c_h, cf_max) # GLM cleaning speed = c_f