AMR for 1D Parabolic Eqs (Clean branch)

[DanielDoehring]

This is a cleaned-up version of #1602.

---

Related to #1147. To keep it simple this is for the moment only 1D.
Essentially, this is a lot copy-paste from the treatment for hyperbolic terms, with the addition of the cache_viscous datastructure for dynamic resizing.
I tried to mark the additions for parabolic terms compared to the hyperbolic ones wherever they occur.
MPI has not (yet) been tested.

Convergence test:
For this elixir

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal compressible Navier-Stokes equations

prandtl_number() = 0.72
mu() = 0.01

equations = CompressibleEulerEquations1D(1.4)
equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), Prandtl=prandtl_number(),
                                                          gradient_variables=GradientVariablesPrimitive())

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
               volume_integral=VolumeIntegralWeakForm())

coordinates_min = -1.0
coordinates_max =  1.0

# Create a uniformly refined mesh with periodic boundaries
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=3,
                periodicity=false,
                n_cells_max=30_000) # set maximum capacity of tree data structure

# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion1D`
#       since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion1D`)
#       and by the initial condition (which passes in `CompressibleEulerEquations1D`).
# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000)
function initial_condition_navier_stokes_convergence_test(x, t, equations)
  # Amplitude and shift
  A = 0.5
  c = 2.0

  # convenience values for trig. functions
  pi_x = pi * x[1]
  pi_t = pi * t

  rho = c + A * cos(pi_x) * cos(pi_t)
  v1  = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0)) ) * cos(pi_t)
  p   = rho^2

  return prim2cons(SVector(rho, v1, p), equations)
end

@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations)
  x = x[1]

  # TODO: parabolic
  # we currently need to hardcode these parameters until we fix the "combined equation" issue
  # see also https://github.com/trixi-framework/Trixi.jl/pull/1160
  inv_gamma_minus_one = inv(equations.gamma - 1)
  Pr = prandtl_number()
  mu_ = mu()

  # Same settings as in `initial_condition`
  # Amplitude and shift
  A = 0.5
  c = 2.0

  # convenience values for trig. functions
  pi_x = pi * x
  pi_t = pi * t

  # compute the manufactured solution and all necessary derivatives
  rho    =  c  + A * cos(pi_x) * cos(pi_t)
  rho_t  = -pi * A * cos(pi_x) * sin(pi_t)
  rho_x  = -pi * A * sin(pi_x) * cos(pi_t)
  rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t)

  v1    =       log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t)
  v1_t  = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t)
  v1_x  =       (A * log(x + 2.0) * exp(-A * (x - 1.0)) + (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t)
  v1_xx = (( 2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0)
                         - A * A * log(x + 2.0) * exp(-A * (x - 1.0))
                         - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t))

  p    = rho * rho
  p_t  = 2.0 * rho * rho_t
  p_x  = 2.0 * rho * rho_x
  p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x

  # Note this simplifies slightly because the ansatz assumes that v1 = v2
  E   = p * inv_gamma_minus_one + 0.5 * rho * v1^2
  E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t
  E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x

  # Some convenience constants
  T_const = equations.gamma * inv_gamma_minus_one / Pr
  inv_rho_cubed = 1.0 / (rho^3)

  # compute the source terms
  # density equation
  du1 = rho_t + rho_x * v1 + rho * v1_x

  # y-momentum equation
  du2 = ( rho_t * v1 + rho * v1_t 
         + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x
    # stress tensor from y-direction
         - v1_xx * mu_)

  # total energy equation
  du3 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x)
    # stress tensor and temperature gradient terms from x-direction
                                - v1_xx * v1   * mu_
                                - v1_x  * v1_x * mu_
         - T_const * inv_rho_cubed * (        p_xx * rho   * rho
                                      - 2.0 * p_x  * rho   * rho_x
                                      + 2.0 * p    * rho_x * rho_x
                                      -       p    * rho   * rho_xx ) * mu_ )

  return SVector(du1, du2, du3)
end

initial_condition = initial_condition_navier_stokes_convergence_test

# BC types
velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2])

heat_bc_left = Isothermal((x, t, equations) -> 
                          Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), 
                                            equations_parabolic))
heat_bc_right = Adiabatic((x, t, equations) -> 0.0)

boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_left)
boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_right)

# define inviscid boundary conditions
boundary_conditions = (; x_neg = boundary_condition_slip_wall,
                         x_pos = boundary_condition_slip_wall)

# define viscous boundary conditions
boundary_conditions_parabolic = (; x_neg = boundary_condition_left,
                                   x_pos = boundary_condition_right)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver;
                                             boundary_conditions=(boundary_conditions, boundary_conditions_parabolic),
                                             source_terms=source_terms_navier_stokes_convergence_test)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span `tspan`
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

amr_controller = ControllerThreeLevel(semi, 
                                      IndicatorMax(semi, variable=Trixi.density),
                                      base_level=Trixi.maximum_level(mesh.tree),
                                      med_level=Trixi.maximum_level(mesh.tree)+1, med_threshold=2.0,
                                      max_level=Trixi.maximum_level(mesh.tree)+2, max_threshold=2.2)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, amr_callback)

###############################################################################
# run the simulation

time_int_tol = 1e-10
sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5,
            ode_default_options()..., callback=callbacks)
summary_callback() # print the timer summary

The convergence study returns

####################################################################################################
l2
rho                 rho_v1              rho_e               
error     EOC       error     EOC       error     EOC       
5.95e-05  -         7.00e-05  -         3.06e-04  -         
5.59e-06  3.41      5.54e-06  3.66      2.47e-05  3.63      
4.15e-07  3.75      3.89e-07  3.83      1.80e-06  3.78      
3.04e-08  3.77      2.61e-08  3.90      1.33e-07  3.76      
2.15e-09  3.82      1.69e-09  3.95      9.63e-09  3.78      

mean      3.69      mean      3.84      mean      3.74      
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_e               
error     EOC       error     EOC       error     EOC       
2.39e-04  -         5.45e-04  -         2.30e-03  -         
3.00e-05  2.99      5.18e-05  3.40      2.59e-04  3.15      
2.44e-06  3.62      3.00e-06  4.11      1.73e-05  3.91      
1.61e-07  3.92      1.93e-07  3.96      1.09e-06  3.98      
9.47e-09  4.09      1.22e-08  3.98      6.78e-08  4.01      

mean      3.66      mean      3.86      mean      3.76      
----------------------------------------------------------------------------------------------------

[DanielDoehring]

As an outlook: cache_viscous is for 2D currently of the form

mutable struct CacheViscous2D{uEltype <: Real}
    u_transformed::Array{uEltype, 4}
    # IDEA: Use SVector for fixed sized vectors?
    gradients::Vector{Array{uEltype, 4}}
    flux_viscous::Vector{Array{uEltype, 4}}

    # internal `resize!`able storage
    _u_transformed::Vector{uEltype}
    _gradients::Vector{Vector{uEltype}}
    _flux_viscous::Vector{Vector{uEltype}}

    function CacheViscous2D{uEltype}(n_vars::Integer, n_nodes::Integer, n_elements::Integer) where {uEltype <: Real}
        new(Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements),
            [Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements) for _ in 1:2],
            [Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements) for _ in 1:2],
            Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements),
            [Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements) for _ in 1:2],
            [Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements) for _ in 1:2])
    end   
end

with

# Only one-dimensional `Array`s are `resize!`able in Julia.
# Hence, we use `Vector`s as internal storage and `resize!`
# them whenever needed. Then, we reuse the same memory by
# `unsafe_wrap`ping multi-dimensional `Array`s around the
# internal storage.
function Base.resize!(cache_viscous::Union{CacheViscous2D, CacheViscous3D}, capacity)
    resize!(cache_viscous._u_transformed, capacity)
    for dim in 1:length(cache_viscous._gradients)
      resize!(cache_viscous._gradients[dim], capacity)
      resize!(cache_viscous._flux_viscous[dim], capacity)
    end

    return nothing
end

[codecov]

Codecov Report

Patch coverage: 92.31% and project coverage change: -8.14% ⚠️

Comparison is base (daf18a5) 96.15% compared to head (a74ec8d) 88.00%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1605      +/-   ##
==========================================
- Coverage   96.15%   88.00%   -8.14%     
==========================================
  Files         414      416       +2     
  Lines       33956    34070     +114     
==========================================
- Hits        32647    29983    -2664     
- Misses       1309     4087    +2778     

Flag	Coverage Δ
unittests	88.00% <92.31%> (-8.14%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files Changed	Coverage Δ
src/solvers/dgsem_tree/dg.jl	94.12% <ø> (ø)
src/callbacks_step/amr.jl	82.48% <80.95%> (-13.50%)	⬇️
...dgsem/elixir_navierstokes_convergence_walls_amr.jl	100.00% <100.00%> (ø)
src/callbacks_step/amr_dg1d.jl	97.55% <100.00%> (+0.60%)	⬆️
src/solvers/dgsem_tree/container_viscous_1d.jl	100.00% <100.00%> (ø)
src/solvers/dgsem_tree/dg_1d_parabolic.jl	94.34% <100.00%> (-0.10%)	⬇️

... and 117 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[sloede]

Thanks for starting to tackle this! Note that convergence is very hard to show with AMR, since you essentially would need an analytical AMR controller that does the exact same refinement for different levels. So getting an EOC of 3.5 where 4 is the expectation seems reasonable to me.

Note that we do not support parallel execution with MPI in 1D since it does not pay off (1D is sooooo cheap)

[DanielDoehring]

Note that convergence is very hard to show with AMR, since you essentially would need an analytical AMR controller that does the exact same refinement for different levels.

Yeah, that is also why I chose the IndicatorMax for the study which refines always in the (approximately) same region and not something like IndicatorLöhner or similar.

[ranocha]

Since @sloede is taking care of this, I removed the review requests for Jesse and me.

[sloede]

LGTM! Thanks a lot for this addition!