Merge pull request #116 from amrueda/subcell_positivity_nonconservative

Add support for non nonconservative equations to subcell-limiting

.github/workflows/FormatCheck.yml

@@ -30,7 +30,7 @@ jobs:
         #       format(".")
         run: |
           julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))'
-          julia  -e 'using JuliaFormatter; format(["benchmark", "ext", "src", "utils"])'
+          julia  -e 'using JuliaFormatter; format(["benchmark", "ext", "src", "test", "utils"])'
       - name: Format check
         run: |
           julia -e '

NEWS.md

@@ -14,6 +14,7 @@ for human readability.
 - Wetting and drying feature and examples for 1D and 2D shallow water equations
 - Implementation of the polytropic Euler equations in 2D
 - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh`
+- AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh`
 
 #### Changed
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.5.47-pre"
+version = "0.5.48-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl

@@ -0,0 +1,110 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations2D(1.4)
+
+"""
+    initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
+
+An MHD blast wave modified from:
+- Dominik Derigs, Gregor J. Gassner, Stefanie Walch & Andrew R. Winters (2018)
+  Entropy Stable Finite Volume Approximations for Ideal Magnetohydrodynamics
+  [doi: 10.1365/s13291-018-0178-9](https://doi.org/10.1365/s13291-018-0178-9)
+This setup needs a positivity limiter for the density.
+"""
+function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
+  # setup taken from Derigs et al. DMV article (2018)
+  # domain must be [-0.5, 0.5] x [-0.5, 0.5], γ = 1.4
+  r = sqrt(x[1]^2 + x[2]^2)
+
+  pmax = 10.0
+  pmin = 0.01
+  rhomax = 1.0
+  rhomin = 0.01
+  if r <= 0.09
+    p = pmax
+    rho = rhomax
+  elseif r >= 0.1
+    p = pmin
+    rho = rhomin
+  else
+    p = pmin + (0.1 - r) * (pmax - pmin) / 0.01
+    rho = rhomin + (0.1 - r) * (rhomax - rhomin) / 0.01
+  end
+  v1 = 0.0
+  v2 = 0.0
+  v3 = 0.0
+  B1 = 1.0/sqrt(4.0*pi)
+  B2 = 0.0
+  B3 = 0.0
+  psi = 0.0
+  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell_local_symmetric)
+volume_flux  = (flux_derigs_etal, flux_nonconservative_powell_local_symmetric)
+basis = LobattoLegendreBasis(3)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons=[1],
+                                positivity_variables_nonlinear=[pressure],
+                                positivity_correction_factor=0.1,
+                                bar_states=false)
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg=volume_flux,
+                                                volume_flux_fv=surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-0.5, -0.5)
+coordinates_max = ( 0.5,  0.5)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=10_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval) 
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+cfl = 0.4
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback())
+
+sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks);
+                  dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl → examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl

@@ -7,12 +7,20 @@ using Trixi
 
 # TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 1.0/3.0 * 10^(-3) # equivalent to Re = 3000
+mu() = 1.0/3.0 * 10^(-5) # equivalent to Re = 30,000
 
 equations = CompressibleEulerEquations2D(1.4)
 equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
                                                           Prandtl=prandtl_number())
 
+"""
+A compressible version of the double shear layer initial condition. Adapted from
+Brown and Minion (1995).
+
+- David L. Brown and Michael L. Minion (1995)
+  Performance of Under-resolved Two-Dimensional Incompressible Flow Simulations.
+  [DOI: 10.1006/jcph.1995.1205](https://doi.org/10.1006/jcph.1995.1205)
+"""
 function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D)
   # Shear layer parameters
   k = 80
@@ -22,8 +30,8 @@ function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquatio
   Ms = 0.1 # maximum Mach number
 
   rho = 1.0
-  v1  = x[2] <= 0.5 ? u0 * tanh(k*(x[2]*0.5 - 0.25)) : u0 * tanh(k*(0.75 -x[2]*0.5))
-  v2  = u0 * delta * sin(2*pi*(x[1]*0.5 + 0.25))
+  v1  = x[2] <= 0.5 ? u0 * tanh(k*(x[2] - 0.25)) : u0 * tanh(k*(0.75 -x[2]))
+  v2  = u0 * delta * sin(2*pi*(x[1]+ 0.25))
   p   = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms
 
   return prim2cons(SVector(rho, v1, v2, p), equations)
@@ -47,27 +55,43 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 2.0)
+tspan = (0.0, 1.0)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 50
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
-                                     extra_analysis_integrals=(energy_kinetic,
-                                                               energy_internal,
-                                                               enstrophy))
+analysis_interval = 2000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval=analysis_interval,)
 
+# This uses velocity-based AMR
+@inline function v1(u, equations::CompressibleEulerEquations2D)
+  rho, rho_v1, _, _ = u
+  return rho_v1 / rho
+end
+amr_indicator = IndicatorLöhner(semi, variable=v1)                                          
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 3,
+                                      med_level  = 5, med_threshold=0.2,
+                                      max_level  = 7, max_threshold=0.5)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=50,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.3)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
-                        alive_callback)
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
 
-time_int_tol = 1e-8
-sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
-            ode_default_options()..., callback=callbacks)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
 summary_callback() # print the timer summary

examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl

@@ -0,0 +1,93 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+diffusivity() = 5.0e-4
+equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations)
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0, -1.0)
+coordinates_max = ( 1.0,  1.0,  1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=80_000)
+
+# Define initial condition
+function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation3D)
+  # Store translated coordinate for easy use of exact solution
+  x_trans = x - equation.advection_velocity * t
+
+  nu = diffusivity()
+  c = 1.0
+  A = 0.5
+  L = 2
+  f = 1/L
+  omega = 2 * pi * f
+  scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t)
+  return SVector(scalar)
+end
+initial_condition = initial_condition_diffusive_convergence_test
+
+# define periodic boundary conditions everywhere
+boundary_conditions = boundary_condition_periodic
+boundary_conditions_parabolic = boundary_condition_periodic
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions,
+                                                                  boundary_conditions_parabolic))
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=3,
+                                      med_level=4, med_threshold=1.2,
+                                      max_level=5, max_threshold=1.45)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary