Implement AMR with p4est

examples/2d/elixir_advection_amr_p4est_unstructured_flag.jl

@@ -0,0 +1,95 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_condition = initial_condition_gauss
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(
+  :all => boundary_condition
+)
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+# Deformed rectangle that looks like a waving flag,
+# lower and upper faces are sinus curves, left and right are vertical lines.
+f1(s) = SVector(-5.0, 5 * s - 5.0)
+f2(s) = SVector( 5.0, 5 * s + 5.0)
+f3(s) = SVector(5 * s, -5.0 + 5 * sin(0.5 * pi * s))
+f4(s) = SVector(5 * s,  5.0 + 5 * sin(0.5 * pi * s))
+faces = (f1, f2, f3, f4)
+
+# This creates a mapping that transforms [-1, 1]^2 to the domain with the faces defined above.
+# It generally doesn't work for meshes loaded from mesh files because these can be meshes
+# of arbitrary domains, but the mesh below is specifically built on the domain [-1, 1]^2.
+Trixi.validate_faces(faces)
+mapping_flag = Trixi.transfinite_mapping(faces)
+
+# Unstructured mesh with 24 cells of the square domain [-1, 1]^n
+mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
+isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
+                              mesh_file)
+
+mesh = P4estMesh(mesh_file, polydeg=3,
+                 mapping=mapping_flag,
+                 initial_refinement_level=1)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+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_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+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=1,
+                                      med_level=2, med_threshold=0.1,
+                                      max_level=3, max_threshold=0.6)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=0.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, save_solution,
+                        amr_callback, stepsize_callback);
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1, # 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/2d/elixir_advection_amr_solution_independent.jl

@@ -17,8 +17,8 @@ function IndicatorSolutionIndependent(semi)
   return IndicatorSolutionIndependent{typeof(cache)}(cache)
 end
 function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any,4},
-                                             equations, dg, cache;
-                                             t, kwargs...)
+                                                   equations, dg, cache;
+                                                   t, kwargs...)
 
   mesh = indicator.cache.mesh
   alpha = indicator.cache.alpha
@@ -86,8 +86,8 @@ initial_condition = initial_condition_gauss
 
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
 
-coordinates_min = (-5, -5)
-coordinates_max = ( 5,  5)
+coordinates_min = (-5.0, -5.0)
+coordinates_max = ( 5.0,  5.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level=4,
                 n_cells_max=30_000)

examples/2d/elixir_advection_amr_solution_independent_p4est.jl

@@ -0,0 +1,153 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+# define new structs inside a module to allow re-evaluating the file
+module TrixiExtension
+
+using Trixi
+
+struct IndicatorSolutionIndependent{Cache<:NamedTuple} <: Trixi.AbstractIndicator
+  cache::Cache
+end
+
+function IndicatorSolutionIndependent(semi)
+  basis = semi.solver.basis
+  alpha = Vector{real(basis)}()
+  cache = (; semi.mesh, alpha)
+  return IndicatorSolutionIndependent{typeof(cache)}(cache)
+end
+
+function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any,4},
+                                                   equations, dg, cache;
+                                                   t, kwargs...)
+
+  mesh = indicator.cache.mesh
+  alpha = indicator.cache.alpha
+  resize!(alpha, nelements(dg, cache))
+
+  #Predict the theoretical center.
+  advectionvelocity = (1.0, 1.0)
+  center = t.*advectionvelocity
+
+  inner_distance = 1
+  outer_distance = 1.85
+
+  #Iterate over all elements
+  for element in 1:length(alpha)
+    # Calculate periodic distance between cell and center.
+    # This requires an uncurved mesh!
+    coordinates = SVector(0.5 * (cache.elements.node_coordinates[1, 1, 1, element] +
+                                 cache.elements.node_coordinates[1, end, 1, element]),
+                          0.5 * (cache.elements.node_coordinates[2, 1, 1, element] +
+                                 cache.elements.node_coordinates[2, 1, end, element]))
+
+    #The geometric shape of the amr should be preserved when the base_level is increased.
+    #This is done by looking at the original coordinates of each cell.
+    cell_coordinates = original_coordinates(coordinates, 5/8)
+    cell_distance = periodic_distance_2d(cell_coordinates, center, 10)
+    if cell_distance < (inner_distance+outer_distance)/2
+      cell_coordinates = original_coordinates(coordinates, 5/16)
+      cell_distance = periodic_distance_2d(cell_coordinates, center, 10)
+    end
+
+    #Set alpha according to cells position inside the circles.
+    target_level = (cell_distance < inner_distance) + (cell_distance < outer_distance)
+    alpha[element] = target_level/2
+  end
+  return alpha
+end
+
+# For periodic domains, distance between two points must take into account
+# periodic extensions of the domain
+function periodic_distance_2d(coordinates, center, domain_length)
+  dx = coordinates .- center
+  dx_shifted = abs.(dx .% domain_length)
+  dx_periodic = min.(dx_shifted, domain_length .- dx_shifted)
+  return sqrt(sum(dx_periodic.^2))
+end
+
+#This takes a cells coordinates and transforms them into the coordinates of a parent-cell it originally refined from.
+#It does it so that the parent-cell has given cell_length.
+function original_coordinates(coordinates, cell_length)
+  offset = coordinates .% cell_length
+  offset_sign = sign.(offset)
+  border = coordinates - offset
+  center = border + (offset_sign .* cell_length/2)
+  return center
+end
+
+end # module TrixiExtension
+
+import .TrixiExtension
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+# advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_condition = initial_condition_gauss
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-5.0, -5.0)
+coordinates_max = ( 5.0,  5.0)
+
+trees_per_dimension = (1, 1)
+
+mesh = P4estMesh(trees_per_dimension, polydeg=3,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 initial_refinement_level=4)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+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_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, TrixiExtension.IndicatorSolutionIndependent(semi),
+                                      base_level=4,
+                                      med_level=5, med_threshold=0.1,
+                                      max_level=6, max_threshold=0.6)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, 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

examples/2d/elixir_advection_p4est_non_conforming.jl → examples/2d/elixir_advection_p4est_non_conforming_flag.jl

@@ -19,14 +19,13 @@ f2(s) = SVector( 1.0, s + 1.0)
 f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
 f4(s) = SVector(s,  1.0 + sin(0.5 * pi * s))
 
+# Create P4estMesh with 3 x 2 trees and 6 x 4 elements
 trees_per_dimension = (3, 2)
-
-# Create P4estMesh with 8 x 8 trees and 16 x 16 elements
 mesh = P4estMesh(trees_per_dimension, polydeg=3,
                  faces=(f1, f2, f3, f4),
                  initial_refinement_level=1)
 
-# Refine bottom left quadrant of each forest to level 4
+# Refine bottom left quadrant of each tree to level 4
 function refine_fn(p4est, which_tree, quadrant)
   if quadrant.x == 0 && quadrant.y == 0 && quadrant.level < 4
     # return true (refine)
@@ -37,7 +36,7 @@ function refine_fn(p4est, which_tree, quadrant)
   end
 end
 
-# Refine recursively until each bottom left quadrant of a forest has level 4
+# Refine recursively until each bottom left quadrant of a tree has level 4
 # The mesh will be rebalanced before the simulation starts
 refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))
 Trixi.p4est_refine(mesh.p4est, true, refine_fn_c, C_NULL)

examples/2d/elixir_advection_p4est_non_conforming_flag_unstructured.jl

@@ -0,0 +1,96 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_condition = initial_condition_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(
+  :all => boundary_condition
+)
+
+# 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)
+
+# Deformed rectangle that looks like a waving flag,
+# lower and upper faces are sinus curves, left and right are vertical lines.
+f1(s) = SVector(-1.0, s - 1.0)
+f2(s) = SVector( 1.0, s + 1.0)
+f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
+f4(s) = SVector(s,  1.0 + sin(0.5 * pi * s))
+faces = (f1, f2, f3, f4)
+
+Trixi.validate_faces(faces)
+mapping_flag = Trixi.transfinite_mapping(faces)
+
+# Unstructured mesh with 24 cells of the square domain [-1, 1]^n
+mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
+isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
+                              mesh_file)
+
+mesh = P4estMesh(mesh_file, polydeg=3,
+                 mapping=mapping_flag,
+                 initial_refinement_level=0)
+
+# Refine bottom left quadrant of each tree to level 3
+function refine_fn(p4est, which_tree, quadrant)
+  if quadrant.x == 0 && quadrant.y == 0 && quadrant.level < 3
+    # return true (refine)
+    return Cint(1)
+  else
+    # return false (don't refine)
+    return Cint(0)
+  end
+end
+
+# Refine recursively until each bottom left quadrant of a tree has level 4
+# The mesh will be rebalanced before the simulation starts
+refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))
+Trixi.p4est_refine(mesh.p4est, true, refine_fn_c, C_NULL)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 0.2));
+
+# 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_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed 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);
+
+# Print the timer summary
+summary_callback()