HDG for advection-diffusion on the sphere

Project.toml

@@ -9,6 +9,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
+GridapHybrid = "94e6dd8f-308f-4a47-89a7-811b3127ae57"
 GridapP4est = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9"
 GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
 GridapPardiso = "34aa2546-dee6-11e9-014e-739fa02ec06f"
@@ -27,13 +28,12 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 GridapDistributed = "=0.2.4"
-Gridap = "=0.17.7"
 GridapP4est = "= 0.1.2"
 GridapPETSc = "=0.4.1"
-PartitionedArrays = "0.2.9"
 MPI = "=0.19"
-WriteVTK = "=1.14.0"
+PartitionedArrays = "0.2.9"
 SparseMatricesCSR = "=0.6.6"
+WriteVTK = "=1.14.0"
 julia = "1.5"
 
 [extras]

driver/sequential/SBRAdvectionHDG.jl

@@ -0,0 +1,50 @@
+module SBRAdvectionHDG
+
+using Gridap
+using GridapGeosciences
+using GridapHybrid
+
+function Gridap.Geometry.push_normal(invJt,n)
+  v = invJt⋅n
+  m = sqrt(inner(v,v))
+  if m < eps()
+    return zero(v)
+  else
+    return v/m
+  end
+end
+
+# solid body rotation on the sphere for the scalar flux form
+# advection equation using a hybridized discontinuous Galerkin
+# method for L2 elements
+
+order  = 1
+degree = 3
+n      = 4
+dt     = 0.5*2.0*π/((order+1)*4*n)
+# single rotation about the sphere
+nstep  = Int(2.0*π/dt)
+
+# solid body rotation velocity field
+function u₀(xyz)
+  θϕr = xyz2θϕr(xyz)
+  u   = cos(θϕr[2])
+  spherical_to_cartesian_matrix(θϕr)⋅VectorValue(u,0,0)
+end
+
+# Gaussian tracer initial condition
+function p₀(xyz)
+  rsq = (xyz[1] - 1.0)*(xyz[1] - 1.0) + xyz[2]*xyz[2] + xyz[3]*xyz[3]
+  exp(-4.0*rsq)
+end
+
+model = CubedSphereDiscreteModel(n,2; radius=1)
+
+advection_hdg(model, order, degree,
+              u₀, p₀, dt, nstep;
+              write_solution=true,
+              write_solution_freq=Int(nstep/8),
+              write_diagnostics=true,
+              write_diagnostics_freq=1,
+              dump_diagnostics_on_screen=true)
+end

src/AdvectionHDG.jl

@@ -0,0 +1,175 @@
+function advection_hdg_time_step!(
+     pn, um, un, model, dΩ, ∂K, d∂K, X, Y, dt,
+     assem=SparseMatrixAssembler(SparseMatrixCSC{Float64,Int},Vector{Float64},X,Y))
+
+  # Second order implicit advection
+  # References:
+  #   Muralikrishnan, Tran, Bui-Thanh, JCP, 2020 vol. 367
+  #   Kang, Giraldo, Bui-Thanh, JCP, 2020 vol. 401
+  γ     = 0.5*(2.0 - sqrt(2.0))
+  γdt   = γ*dt
+  γm1dt = (1.0 - γ)*dt
+
+  n  = get_cell_normal_vector(∂K)
+  nₒ = get_cell_owner_normal_vector(∂K)
+
+  # First stage
+  b₁((q,m)) = ∫(q*pn)dΩ - 
+              ∫(m*0.0)d∂K
+
+  a₁((p,l),(q,m)) = ∫(q*p - γdt*(∇(q)⋅un)*p)dΩ + ∫(((un⋅n) + abs(un⋅n))*γdt*q*p)d∂K -    # [q,p] block
+                    ∫(γdt*abs(un⋅n)*q*l)d∂K +                                            # [q,l] block
+                    ∫(((un⋅n) + abs(un⋅n))*p*m)d∂K -                                     # [m,p] block
+                    ∫(abs(un⋅n)*l*m)d∂K                                                  # [m,l] block
+
+  op₁   = HybridAffineFEOperator((u,v)->(a₁(u,v),b₁(v)), X, Y, [1], [2])
+  Xh    = solve(op₁)
+  ph,lh = Xh
+
+  # Second stage
+  b₂((q,m)) = ∫(q*pn + γm1dt*(∇(q)⋅un)*ph)dΩ - ∫(((un⋅n) + abs(un⋅n))*γm1dt*ph*q - abs(un⋅n)*γm1dt*lh*q)d∂K -
+              ∫(0.5*((un⋅n) + abs(un⋅n))*ph*m - 0.5*abs(un⋅n)*lh*m)d∂K
+
+  a₂((p,l),(q,m)) = ∫(q*p - γdt*(∇(q)⋅un)*p)dΩ + ∫(((un⋅n) + abs(un⋅n))*γdt*q*p)d∂K -    # [q,p] block
+                    ∫(γdt*abs(un⋅n)*q*l)d∂K +                                            # [q,l] block
+                    ∫(((un⋅n) + abs(un⋅n))*0.5*p*m)d∂K -                                 # [m,p] block
+                    ∫(0.5*abs(un⋅n)*l*m)d∂K                                              # [m,l] block
+
+  op₂   = HybridAffineFEOperator((u,v)->(a₂(u,v),b₂(v)), X, Y, [1], [2])
+  Xm    = solve(op₂)
+  pm,_  = Xm
+
+  get_free_dof_values(pn) .= get_free_dof_values(pm)
+end
+
+function project_initial_conditions(dΩ, P, Q, p₀, U, V, u₀, mass_matrix_solver)
+  # the tracer
+  b₁(q)    = ∫(q*p₀)dΩ
+  a₁(p,q)  = ∫(q*p)dΩ
+  rhs₁     = assemble_vector(b₁, Q)
+  L2MM     = assemble_matrix(a₁, P, Q)
+  L2MMchol = numerical_setup(symbolic_setup(mass_matrix_solver,L2MM),L2MM)
+  pn       = FEFunction(Q, copy(rhs₁))
+  pnv      = get_free_dof_values(pn)
+
+  solve!(pnv, L2MMchol, pnv)
+
+  # the velocity field
+  b₂(v)    = ∫(v⋅u₀)dΩ
+  a₂(u,v)  = ∫(v⋅u)dΩ
+  rhs₂     = assemble_vector(b₂, V)
+  MM       = assemble_matrix(a₂, U, V)
+  MMchol   = numerical_setup(symbolic_setup(mass_matrix_solver,MM),MM)
+  un       = FEFunction(V, copy(rhs₂))
+  unv      = get_free_dof_values(un)
+
+  solve!(unv, MMchol, unv)
+
+  pn, pnv, L2MM, un
+end
+
+function conservation(L2MM, pnv, p_tmp, mass_0, entropy_0)
+  mul!(p_tmp, L2MM, pnv)
+  mass_con    = sum(p_tmp)
+  entropy_con = p_tmp⋅pnv
+  mass_con    = (mass_con - mass_0)/mass_0
+  entropy_con = (entropy_con - entropy_0)/entropy_0
+  mass_con, entropy_con
+end
+
+function advection_hdg(
+        model, order, degree,
+        u₀, p₀, dt, N;
+        mass_matrix_solver::Gridap.Algebra.LinearSolver=Gridap.Algebra.BackslashSolver(),
+        write_diagnostics=true,
+        write_diagnostics_freq=1,
+        dump_diagnostics_on_screen=true,
+        write_solution=false,
+        write_solution_freq=N/10,
+        output_dir="adv_eq_ncells_$(num_cells(model))_order_$(order)")
+
+  # Forward integration of the advection equation
+  D = num_cell_dims(model)
+  Ω = Triangulation(ReferenceFE{D},model)
+  Γ = Triangulation(ReferenceFE{D-1},model)
+  ∂K = GridapHybrid.Skeleton(model)
+
+  reffeᵤ = ReferenceFE(lagrangian,VectorValue{3,Float64},order;space=:P)
+  reffeₚ = ReferenceFE(lagrangian,Float64,order;space=:P)
+  reffeₗ = ReferenceFE(lagrangian,Float64,order;space=:P)
+
+  # Define test FESpaces
+  V = TestFESpace(Ω, reffeᵤ; conformity=:L2)
+  Q = TestFESpace(Ω, reffeₚ; conformity=:L2)
+  M = TestFESpace(Γ, reffeₗ; conformity=:L2)
+  Y = MultiFieldFESpace([Q,M])
+
+  U = TrialFESpace(V)
+  P = TrialFESpace(Q)
+  L = TrialFESpace(M)
+  X = MultiFieldFESpace([P,L])
+
+  dΩ  = Measure(Ω,degree)
+  d∂K = Measure(∂K,degree)
+
+  @printf("time step: %14.9e\n", dt)
+  @printf("number of time steps: %u\n", N)
+
+  # Project the initial conditions onto the trial spaces
+  pn, pnv, L2MM, un = project_initial_conditions(dΩ, P, Q, p₀, U, V, u₀, mass_matrix_solver)
+
+  # Work array
+  p_tmp = copy(pnv)
+  p_ic = copy(pnv)
+
+  # Initial states
+  mul!(p_tmp, L2MM, pnv)
+  p01 = sum(p_tmp)  # total mass
+  p02 = p_tmp⋅pnv   # total entropy
+
+  function run_simulation(pvd=nothing)
+    diagnostics_file = joinpath(output_dir,"adv_diagnostics.csv")
+    if (write_diagnostics)
+      initialize_csv(diagnostics_file,"step", "mass", "entropy")
+    end
+    if (write_solution && write_solution_freq>0)
+      pvd[dt*Float64(0)] = new_vtk_step(Ω,joinpath(output_dir,"n=0"),["pn"=>pn,"un"=>un])
+    end
+
+    for istep in 1:N
+      advection_hdg_time_step!(pn, u₀, u₀, model, dΩ, ∂K, d∂K, X, Y, dt)
+
+      if (write_diagnostics && write_diagnostics_freq>0 && mod(istep, write_diagnostics_freq) == 0)
+        # compute mass and entropy conservation
+        pn1, pn2 = conservation(L2MM, pnv, p_tmp, p01, p02)
+        if dump_diagnostics_on_screen
+          @printf("%5d\t%14.9e\t%14.9e\n", istep, pn1, pn2)
+        end
+        append_to_csv(diagnostics_file; step=istep, mass=pn1, entropy=pn2)
+      end
+      if (write_solution && write_solution_freq>0 && mod(istep, write_solution_freq) == 0)
+        pvd[dt*Float64(istep)] = new_vtk_step(Ω,joinpath(output_dir,"n=$(istep)"),["pn"=>pn,"un"=>un])
+      end
+    end
+    # compute the L2 error with respect to the inital condition after one rotation
+    mul!(p_tmp, L2MM, p_ic)
+    l2_norm_sq = p_tmp⋅p_ic
+    p_ic .= p_ic .- pnv
+    mul!(p_tmp, L2MM, p_ic)
+    l2_err_sq = p_ic⋅p_tmp
+    l2_err = sqrt(l2_err_sq/l2_norm_sq)
+    @printf("L2 error: %14.9e\n", l2_err)
+    pn1, pn2 = conservation(L2MM, pnv, p_tmp, p01, p02)
+    l2_err, pn1, pn2
+  end
+  if (write_diagnostics || write_solution)
+    rm(output_dir,force=true,recursive=true)
+    mkdir(output_dir)
+  end
+  if (write_solution)
+    pvdfile=joinpath(output_dir,"adv_eq_ncells_$(num_cells(model))_order_$(order)")
+    paraview_collection(run_simulation,pvdfile)
+  else
+    run_simulation()
+  end
+end

src/GridapGeosciences.jl

@@ -1,5 +1,6 @@
 module GridapGeosciences
   using Gridap
+  using GridapHybrid
   using FillArrays
   using SparseArrays
   using LinearAlgebra
@@ -20,6 +21,7 @@ module GridapGeosciences
   include("ThermalShallowWaterExplicit_MatAdv.jl")
   include("ShallowWaterThetaMethodFullNewton.jl")
   include("Helpers.jl")
+  include("AdvectionHDG.jl")
   export rₑ, Ωₑ, g, f
   export CubedSphereDiscreteModel
   export perp,⟂
@@ -44,6 +46,7 @@ module GridapGeosciences
   export shallow_water_theta_method_full_newton_time_stepper
   export write_to_csv, get_scalar_field_from_csv, append_to_csv, initialize_csv
   export clone_fe_function, setup_mixed_spaces, setup_and_factorize_mass_matrices, new_vtk_step
+  export advection_hdg
 
   # MPI-parallel part
   using PartitionedArrays

src/mpi/CubedSphereDistributedDiscreteModels.jl

@@ -314,7 +314,8 @@ function Gridap.Geometry.get_glue(a::D2toD3AnalyticalMapCubedSphereTriangulation
   tface_to_mface=Gridap.Fields.IdentityVector(nc)
   tface_to_mface_map=Fill(Gridap.Fields.GenericField(identity),nc)
   mface_to_tface=tface_to_mface
-  Gridap.Geometry.FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
+  Dt=2
+  Gridap.Geometry.FaceToFaceGlue(Dt,tface_to_mface,tface_to_mface_map,mface_to_tface)
 end
 
 Gridap.Geometry.get_grid(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = trian.model.cube_grid_geo

test/sequential/SBRAdvectionHDGTests.jl

@@ -0,0 +1,59 @@
+module SBRAdvectionHDGTests
+
+using Test
+using Gridap
+using GridapGeosciences
+using GridapHybrid
+
+function Gridap.Geometry.push_normal(invJt,n)
+  v = invJt⋅n
+  m = sqrt(inner(v,v))
+  if m < eps()
+    return zero(v)
+  else
+    return v/m
+  end
+end
+
+# solid body rotation velocity field
+function u₀(xyz)
+  θϕr = xyz2θϕr(xyz)
+  u   = cos(θϕr[2])
+  spherical_to_cartesian_matrix(θϕr)⋅VectorValue(u,0,0)
+end
+
+# Gaussian tracer initial condition
+function p₀(xyz)
+  rsq = (xyz[1] - 1.0)*(xyz[1] - 1.0) + xyz[2]*xyz[2] + xyz[3]*xyz[3]
+  exp(-4.0*rsq)
+end
+
+# solid body rotation on the sphere for the scalar flux form
+# advection equation using a hybridized discontinuous Galerkin
+# method for L2 elements
+
+order  = 1
+degree = 3
+l2_err_i   = [2.990077828e-01 , 9.625369737e-02 , 1.918189011e-02 ]
+enst_con_i = [-2.533366481e-01, -7.330337526e-02, -1.289758067e-02]
+
+for i in 1:3
+  n      = 2*2^i
+  dt     = 0.5*2.0*π/((order+1)*4*n)
+  # single rotation about the sphere
+  nstep  = Int(2.0*π/dt)
+
+  model = CubedSphereDiscreteModel(n,2; radius=1)
+
+  l2_err, mass_con, enst_con = advection_hdg(model, order, degree,
+                                             u₀, p₀, dt, nstep;
+                                             write_solution=false,
+                                             write_solution_freq=Int(nstep/8),
+                                             write_diagnostics=true,
+                                             write_diagnostics_freq=1,
+                                             dump_diagnostics_on_screen=true)
+  @test abs(l2_err - l2_err_i[i])     < 10.0^-10
+  @test abs(mass_con)                 < 10.0^-12
+  @test abs(enst_con - enst_con_i[i]) < 10.0^-10
+end
+end

test/sequential/runtests.jl

@@ -1,6 +1,7 @@
 using GridapGeosciences
 using Test
 @testset "GridapGeosciences" begin
+  @time @testset "SBRAdvectionHDGTests" begin include("SBRAdvectionHDGTests.jl") end
   @time @testset "CubedSphereDiscreteModelsTests" begin include("CubedSphereDiscreteModelsTests.jl") end
   @time @testset "DarcyCubedSphereTests" begin include("DarcyCubedSphereTests.jl") end
   @time @testset "WeakDivPerpTests" begin include("WeakDivPerpTests.jl") end