Parabolic mortars for P4estMesh

examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl

@@ -0,0 +1,117 @@
+using OrdinaryDiffEq, Plots
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 1.0e-3
+advection_velocity = (1.0, 0.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+# Define initial condition (copied from "examples/tree_1d_dgsem/elixir_advection_diffusion.jl")
+function initial_condition_bubble(x, t, equation)
+    return SVector((x[1]-pi) * (x[2]-pi) * (x[1]+pi) * (x[2]+pi))
+end
+
+initial_condition = initial_condition_bubble
+
+boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
+                           :y_neg => BoundaryConditionDirichlet(initial_condition),
+                           :y_pos => BoundaryConditionDirichlet(initial_condition),
+                           :x_pos => boundary_condition_do_nothing)
+
+boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), 
+                                     :x_pos => BoundaryConditionDirichlet(initial_condition), 
+                                     :y_neg => BoundaryConditionDirichlet(initial_condition), 
+                                     :y_pos => BoundaryConditionDirichlet(initial_condition))
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=2, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-pi, -pi) # minimum coordinates (min(x), min(y))
+coordinates_max = ( pi,  pi) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=2, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=false)
+
+# 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.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, .5)
+ode = semidiscretize(semi, tspan)
+
+# u = sol.u[end]
+
+# du = similar(u)
+# Trixi.rhs_parabolic!(du, u, semi, 0.0)
+
+# x, y = [semi.cache.elements.node_coordinates[i, :, :, :] for i in 1:2]
+# for i in eachindex(x)
+#     u[i] = initial_condition_bubble((x[i], y[i]), 0.0, equations)[1]
+# end
+# u = Trixi.wrap_array(u, semi)
+# fill!(cache_parabolic.elements.surface_flux_values, NaN);
+# dg = solver
+# parabolic_scheme = semi.solver_parabolic
+# t = 0.0
+# (; cache, cache_parabolic, boundary_conditions_parabolic) = semi
+# @unpack viscous_container = cache_parabolic
+# @unpack u_transformed, gradients, flux_viscous = viscous_container
+
+# Trixi.transform_variables!(u_transformed, u, mesh, equations_parabolic,
+#                            dg, parabolic_scheme, cache, cache_parabolic)
+
+# Trixi.calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic,
+#                      boundary_conditions_parabolic, dg, cache, cache_parabolic)
+
+# grad_x, grad_y = gradients
+# @show any(isnan.(grad_x))
+# @show any(isnan.(grad_y))
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=1,
+                                      med_level=2, med_threshold=0.5,
+                                      max_level=3, max_threshold=0.75)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5)
+
+# 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
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+time_int_tol = 1.0e-11
+sol = solve(ode, dt = 1e-7, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+
+# Print the timer summary
+summary_callback()
+plot(sol)
+# pd = PlotData2D(sol)
+# plot!(getmesh(pd))
+

examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_amr.jl

@@ -0,0 +1,101 @@
+using OrdinaryDiffEq, Plots
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 1.0e-2
+advection_velocity = (1.0, 1.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+function x_trans_periodic(x, domain_length=SVector(2 * pi), center=SVector(0.0))
+    x_normalized = x .- center
+    x_shifted = x_normalized .% domain_length
+    x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* domain_length
+    return center + x_shifted + x_offset
+end
+
+# Define initial condition (copied from "examples/tree_1d_dgsem/elixir_advection_diffusion.jl")
+function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D)
+    # Store translated coordinate for easy use of exact solution
+    # Assumes that advection_velocity[2] = 0 (effectively that we are solving a 1D equation)
+    x_trans = x_trans_periodic(x[2] - equation.advection_velocity[2] * t)
+    # y_trans = x_trans_periodic(x[1] - equation.advection_velocity[1] * t)
+
+    nu = diffusivity()
+    c = 0.0
+    A = 1.0
+    omega = 1.0
+    scalar = c + A * sin(omega * (sum(x_trans))) * exp(-nu * omega^2 * t)
+    return SVector(scalar)
+end
+
+# Define initial condition (copied from "examples/tree_1d_dgsem/elixir_advection_diffusion.jl")
+function initial_condition_new_test(x, t, equation::LinearScalarAdvectionEquation2D)
+    return SVector(2 * x[1] + x[2] > 0)
+end
+# initial_condition = initial_condition_diffusive_convergence_test
+initial_condition = initial_condition_new_test
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=2, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-pi, -pi) # minimum coordinates (min(x), min(y))
+coordinates_max = ( pi,  pi) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=2, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=true)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=2,
+                                      med_level=3, med_threshold=0.5,
+                                      max_level=4, max_threshold=0.75)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5)
+
+# 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
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+time_int_tol = 1.0e-11
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+# sol = solve(ode, ROCK4(eigen_est=eigen_est); abstol=time_int_tol, reltol=time_int_tol,
+#             ode_default_options()..., callback=callbacks)
+# Print the timer summary
+summary_callback()
+plot(sol)
+pd = PlotData2D(sol)
+plot!(getmesh(pd))
+

examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl

@@ -0,0 +1,102 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 0.001
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
+                                                          Prandtl=prandtl_number())
+
+# 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)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh
+trees_per_dimension = (6, 6)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=(false, false))
+
+function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D)
+  Ma = 0.5
+  rho = 1.0
+  u, v = 0.0, 0.0
+  p = 1.0 / (Ma^2 * equations.gamma)
+  return prim2cons(SVector(rho, u, v, p), equations)
+end
+initial_condition = initial_condition_cavity
+
+# BC types
+velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0))
+velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0))
+heat_bc = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
+boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)
+
+# define periodic boundary conditions everywhere
+boundary_conditions = Dict( :x_neg => boundary_condition_slip_wall,
+                            :y_neg => boundary_condition_slip_wall,
+                            :y_pos => boundary_condition_slip_wall,
+                            :x_pos => boundary_condition_slip_wall)
+
+boundary_conditions_parabolic = Dict( :x_neg => boundary_condition_cavity,
+                                      :y_neg => boundary_condition_cavity,
+                                      :y_pos => boundary_condition_lid,
+                                      :x_pos => boundary_condition_cavity)
+
+# 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))
+
+# semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+#                                              initial_condition, solver;)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 25.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=100)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi, variable=Trixi.density)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level=0,
+                                      med_level=1, med_threshold=0.02,
+                                      max_level=3, max_threshold=0.05)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+callbacks = CallbackSet(summary_callback, alive_callback,analysis_callback, amr_callback)
+# callbacks = CallbackSet(summary_callback, alive_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)
+summary_callback() # print the timer summary
+
+
+pd = PlotData2D(sol)
+plot(pd["rho"])
+plot!(getmesh(pd))