diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a14c19f..8832b6d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -5,8 +5,8 @@ on: - create jobs: - tutorials_set_1: - name: Tutorials1 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + tutorials: + name: Tutorials ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - ${{ matrix.test_set.name }} runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -17,6 +17,21 @@ jobs: - ubuntu-latest arch: - x64 + test_set: + - name: Basics + files: "poisson.jl elasticity.jl lagrange_multipliers.jl dg_discretization.jl darcy.jl stokes.jl advection_diffusion.jl" + - name: Nonlinear + files: "p_laplacian.jl hyperelasticity.jl inc_navier_stokes.jl" + - name: Applications + files: "fsi_tutorial.jl emscatter.jl isotropic_damage.jl TopOptEMFocus.jl" + - name: Transient + files: "isotropic_damage.jl transient_linear.jl transient_nonlinear.jl " + - name: Hybrid + files: "poisson_hdg.jl poisson_hho.jl" + - name: Plugins + files: "stokes_blocks.jl poisson_amr.jl poisson_unfitted.jl" + - name: Other + files: "validation.jl validation_DrWatson.jl interpolation_fe.jl poisson_dev_fe.jl geometry_dev.jl" steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 @@ -25,72 +40,11 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl poisson.jl validation.jl elasticity.jl lagrange_multipliers.jl - tutorials_set_2: - name: Tutorials2 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1.10' - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl p_laplacian.jl hyperelasticity.jl dg_discretization.jl fsi_tutorial.jl poisson_amr.jl - tutorials_set_3: - name: Tutorials3 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1.10' - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl darcy.jl inc_navier_stokes.jl stokes.jl poisson_dev_fe.jl emscatter.jl - tutorials_set_4: - name: Tutorials4 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1.10' - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl isotropic_damage.jl validation_DrWatson.jl interpolation_fe.jl transient_linear.jl transient_nonlinear.jl TopOptEMFocus.jl + - run: | + cd test/ + julia --project=.. --color=yes --check-bounds=yes runtests.jl ${{ matrix.test_set.files }} tutorials_mpi: - name: TutorialsMPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + name: Tutorials MPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -109,7 +63,9 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests_mpi.jl + - run: | + cd test/ + julia --project=.. --color=yes --check-bounds=yes runtests_mpi.jl docs: name: Documentation runs-on: ubuntu-latest diff --git a/Project.toml b/Project.toml index 506480e..e2b01a6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Santiago Badia ", "Francesc Verdugo "TopOptEMFocus.jl", "Poisson on unfitted meshes"=>"poisson_unfitted.jl", "Poisson with AMR"=>"poisson_amr.jl", + "Poisson with HDG"=>"poisson_hdg.jl", + "Poisson with HHO on polytopal meshes"=>"poisson_hho.jl", + "Block assembly and solvers: Incompressible Stokes example"=>"stokes_blocks.jl", "Lagrange multipliers" => "lagrange_multipliers.jl", "Low-level API - Poisson equation"=>"poisson_dev_fe.jl", "Low-level API - Geometry" => "geometry_dev.jl", diff --git a/docs/make.jl b/docs/make.jl index cf642a0..8bd6553 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -60,9 +60,9 @@ for (i,(title,filename)) in enumerate(Tutorials.files) end makedocs( - sitename = "Gridap tutorials", - format = Documenter.HTML(), - pages = pages + sitename = "Gridap tutorials", + format = Documenter.HTML(), + pages = pages ) # Documenter can also automatically deploy documentation to gh-pages. @@ -70,5 +70,5 @@ makedocs( # for more information. deploydocs( - repo = "github.com/gridap/Tutorials.git" + repo = "github.com/gridap/Tutorials.git" ) diff --git a/src/advection_diffusion.jl b/src/advection_diffusion.jl index af58d8e..6900ef3 100644 --- a/src/advection_diffusion.jl +++ b/src/advection_diffusion.jl @@ -82,20 +82,25 @@ # The meshed model is like following: -# ![Pentagon_2D_Mesh](Pentagon_2D_Mesh.png) +# ![Pentagon_2D_Mesh](../assets/advection_diffusion/Pentagon_2D_Mesh.png) # ## Numerical Scheme -# The weak form of a PDE is a reformulation that allows the problem to be solved in a broader function space. Instead of requiring the solution to satisfy the PDE at every point (as in the strong form), the weak form requires the solution to satisfy an integral equation. This makes it particularly suitable for numerical methods such as the Finite Element Method. +# The weak form of a PDE is a reformulation that allows the problem to be solved in a broader +# function space. Instead of requiring the solution to satisfy the PDE at every point (as in the strong form), +# the weak form requires the solution to satisfy an integral equation. This makes it particularly suitable for +# numerical methods such as the Finite Element Method. -# Since we already hcave the original PDE (strong form), we can multiply each side by a test function $v \in H^1(\Omega)$(functions with square-integrable first derivatives). $v$ satisfies the Dirichlet boundary condition of $0$. Make integral on both sides. +# Since we already hcave the original PDE (strong form), we can multiply each side by a test +# function $v \in H^1(\Omega)$(functions with square-integrable first derivatives). $v$ satisfies +# the Dirichlet boundary condition of $0$. Make integral on both sides. # The weak form associated with this formulation is, find $T \in H^1(\Omega)$, such that: -# $$ +# ```math # \int_\Omega v (\mathbf{u} \cdot \nabla T) \, \mathrm{d}\Omega # + \int_\Omega D \nabla T \cdot \nabla v \, \mathrm{d}\Omega # = 0 -# $$ +# ``` # ## FE spaces @@ -106,7 +111,7 @@ using GridapGmsh # Import the model from file using GmshDiscreteModel: -model_pentagon = GmshDiscreteModel("pentagon_mesh.msh") +model_pentagon = GmshDiscreteModel("../models/pentagon_mesh.msh") # Set up the test FE space $V_0$, which conforms the zero value boundary conditions. @@ -165,8 +170,8 @@ writevtk(Ω,"results_non_zero",cellfields=["uh_non_zero"=>uh_non_zero]) # We can use the ParaView to preview the results clearly. Here is the temperature distribution without any flow velocity: -# ![Result_zero](Result_zero.png) +# ![Result_zero](../assets/advection_diffusion/Result_zero.png) # Here is the temperature distribution with a flow velocity of 2 going up: -# ![Result_zero](Result_non_zero.png) +# ![Result_zero](../assets/advection_diffusion/Result_non_zero.png) diff --git a/src/poisson_hdg.jl b/src/poisson_hdg.jl new file mode 100644 index 0000000..aa2b9b5 --- /dev/null +++ b/src/poisson_hdg.jl @@ -0,0 +1,164 @@ +# In this tutorial, we will implement a Hybridizable Discontinuous Galerkin (HDG) method +# for solving the Poisson equation. The HDG method is an efficient variant of DG methods +# that introduces an auxiliary variable m on mesh interfaces to reduce the global system size. +# +# ## HDG Discretization +# +# We consider the Poisson equation with Dirichlet boundary conditions: +# +# ```math +# \begin{aligned} +# -\Delta u &= f \quad \text{in} \quad \Omega\\ +# u &= g \quad \text{in} \quad \partial\Omega +# \end{aligned} +# ``` +# +# The HDG method first rewrites the problem as a first-order system: +# +# ```math +# \begin{aligned} +# \boldsymbol{q} + \nabla u &= \boldsymbol{0} \quad \text{in} \quad \Omega\\ +# \nabla \cdot \boldsymbol{q} &= f \quad \text{in} \quad \Omega\\ +# u &= g \quad \text{on} \quad \partial\Omega +# \end{aligned} +# ``` +# +# The HDG discretization introduces three variables: +# - ``\boldsymbol{q}_h``: the approximation to the flux ``\boldsymbol{q}`` +# - ``u_h``: the approximation to the solution ``u`` +# - ``m_h``: the approximation to the trace of ``u`` on element faces +# +# Numerical fluxes are defindes as +# +# ```math +# \widehat{\boldsymbol{q}}_h = \boldsymbol{q}_h + \tau(u_h - m_h)\boldsymbol{n} +# ``` +# +# where $\tau$ is a stabilization parameter. +# +# First, let's load the required Gridap packages: + +using Gridap +using Gridap.Geometry +using Gridap.FESpaces +using Gridap.MultiField +using Gridap.CellData + +# ## Manufactured Solution +# +# We use the method of manufactured solutions to verify our implementation. +# We choose a solution u and derive the corresponding source term f: + +u(x) = sin(2*π*x[1])*sin(2*π*x[2]) +q(x) = -∇(u)(x) # Define the flux q = -∇u +f(x) = (∇ ⋅ q)(x) # Source term f = -Δu = -∇⋅(∇u)$ + +# ## Geometry +# +# We generate a D-dimensional simplicial mesh from a Cartesian grid: + +D = 2 # Problem dimension +nc = Tuple(fill(8, D)) # 4 cells in each direction +domain = Tuple(repeat([0, 1], D)) # Unit cube domain +model = simplexify(CartesianDiscreteModel(domain,nc)) + +# From this mesh, we will require two triangulations where to define our HDG spaces: +# +# 1. A cell triangulation $\Omega$, for the volume variables +# 2. A face triangulation $\Gamma$, for the skeleton variables +# +# These are given by + +Ω = Triangulation(ReferenceFE{D}, model) # Volume triangulation +Γ = Triangulation(ReferenceFE{D-1}, model) # Skeleton triangulation + +# ## PatchTopology and PatchTriangulation +# +# A key aspect of hybrid methods is the use of static condensation, which is the +# elimination of cell unknowns to reduce the size of the global system. +# To achieve this, we need to be able to assemble and solve local problems on each cell, that +# involve +# - contributions from the cell itself +# - contributions from the cell faces +# To this end, Gridap provides a general framework for patch-assembly and solves. The idea +# is to define a patch decomposition of the mesh (in this case, a patch is a cell and its sourrounding +# faces). We can then gather contributions for each patch, solve the local problems, and +# assemble the results into the global system. +# +# The following code creates the required `PatchTopology` for the problem at hand. We then +# take d-dimensional slices of it by the means of `PatchTriangulation` and `PatchBoundaryTriangulation`. +# These are the `Triangulation`s we will integrate our weakform over. + +ptopo = Geometry.PatchTopology(model) +Ωp = Geometry.PatchTriangulation(model,ptopo) # Patch volume triangulation +Γp = Geometry.PatchBoundaryTriangulation(model,ptopo) # Patch skeleton triangulation + +# ## FESpaces +# +# HDG uses three different finite element spaces: +# 1. A vector-valued space for the flux q (Q) +# 2. A scalar space for the solution u (V) +# 3. A scalar space for the interface variable m (M) +# +# We then define discontinuous finite element spaces of the approppriate order, locally $\mathbb{P}^k$. +# Note that only the skeletal space has Dirichlet boundary conditions. + +order = 1 # Polynomial order +reffe_Q = ReferenceFE(lagrangian, VectorValue{D, Float64}, order; space=:P) +reffe_V = ReferenceFE(lagrangian, Float64, order; space=:P) +reffe_M = ReferenceFE(lagrangian, Float64, order; space=:P) + +V = TestFESpace(Ω, reffe_V; conformity=:L2) # Discontinuous vector space +Q = TestFESpace(Ω, reffe_Q; conformity=:L2) # Discontinuous scalar space +M = TestFESpace(Γ, reffe_M; conformity=:L2, dirichlet_tags="boundary") # Interface space +N = TrialFESpace(M, u) + +# ## MultiField Structure +# +# Since we are doing static condensation, we need assemble by blocks. In particular, the +# `StaticCondensationOperator` expects the variables to be groupped in two blocks: +# - The eliminated variables (in this case, the volume variables q and u) +# - The retained variables (in this case, the interface variable m) +# We will assemble by blocks using the `BlockMultiFieldStyle` API. + +mfs = BlockMultiFieldStyle(2,(2,1)) # Special blocking for efficient static condensation +X = MultiFieldFESpace([V, Q, N]; style=mfs) + +# ## Weak Form and integration +# + +degree = 2*(order+1) # Integration degree +dΩp = Measure(Ωp,degree) # Volume measure, on the patch triangulation +dΓp = Measure(Γp,degree) # Surface measure, on the patch boundary triangulation + +τ = 1.0 # HDG stabilization parameter + +n = get_normal_vector(Γp) # Face normal vector +Πn(u) = u⋅n # Normal component +Π(u) = change_domain(u,Γp,DomainStyle(u)) # Project to skeleton + +a((uh,qh,sh),(vh,wh,lh)) = ∫( qh⋅wh - uh*(∇⋅wh) - qh⋅∇(vh) )dΩp + ∫(sh*Πn(wh))dΓp + + ∫((Πn(qh) + τ*(Π(uh) - sh))*(Π(vh) + lh))dΓp +l((vh,wh,lh)) = ∫( f*vh )dΩp + +# ## Static Condensation and Solution +# +# With all these ingredients, we can now build our statically-condensed operator ans +# solve the problem. Note that we are solving a scatically-condensed system. We can +# retrieve the internal `AffineFEOperator` from `op.sc_op`. + +op = MultiField.StaticCondensationOperator(ptopo,X,a,l) +uh, qh, sh = solve(op) + +dΩ = Measure(Ω,degree) +eh = uh - u +l2_uh = sqrt(sum(∫(eh⋅eh)*dΩ)) + +writevtk(Ω,"results",cellfields=["uh"=>uh,"qh"=>qh,"eh"=>eh]) + +# ## Going Further +# +# By modifying the stabilisation term, HDG can also work on polytopal meshes. An driver +# solving the same problem on a polytopal mesh is available [in the Gridap repository](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HDGPolytopalTests.jl). +# A tutorial for HHO on polytopal meshes is also available. +# diff --git a/src/poisson_hho.jl b/src/poisson_hho.jl new file mode 100644 index 0000000..e1967e4 --- /dev/null +++ b/src/poisson_hho.jl @@ -0,0 +1,282 @@ +# +# In this tutorial, we will learn how to implement a mixed-order Hybrid High-Order (HHO) method +# for solving the Poisson equation. HHO methods are a class of modern hybridizable finite element methods +# that provide optimal convergence rates while enabling static condensation for efficient solution. +# +# HHO is a mathematically complex method. This tutorial will **NOT** cover the method nor +# its mathematical foundations. You should be familiar with both before going into +# the tutorial itself. +# We also recommend going through the HDG tutorial first, which shares many of the +# same concepts but is simpler to understand. +# +# ## Problem statement +# +# We consider the Poisson equation with Dirichlet boundary conditions: +# +# ```math +# \begin{aligned} +# -\Delta u &= f \quad \text{in } \Omega \\ +# u &= g \quad \text{on } \partial\Omega +# \end{aligned} +# ``` +# +# where $\Omega$ is a bounded domain in $\mathbb{R}^2$, f is a source term, and $g$ is the prescribed boundary value. +# +# ## HHO discretization +# +# The HHO method introduces two types of unknowns: +# 1. Cell unknowns defined in the volume of each mesh cell +# 2. Face unknowns defined on the mesh facets +# +# This hybrid structure allows for efficient static condensation by eliminating the cell unknowns +# algebraically at the element level. +# +# We start by loading the required packages + +using Gridap +using Gridap.Geometry, Gridap.FESpaces, Gridap.MultiField +using Gridap.CellData, Gridap.Fields, Gridap.Helpers +using Gridap.ReferenceFEs +using Gridap.Arrays + +# ## Geometry +# +# We generate a 2-dimensional simplicial mesh from a Cartesian grid, then taking its Voronoi +# dual to create a polytopal mesh. + +u(x) = sin(2*π*x[1])*sin(2*π*x[2]) +f(x) = -Δ(u)(x) + +n = 10 +base_model = simplexify(CartesianDiscreteModel((0,1,0,1),(n,n))) +model = Geometry.voronoi(base_model) + +# From this mesh, we will require two triangulations where to define our HDG spaces: +# +# 1. A cell triangulation $\Omega$, for the volume variables +# 2. A face triangulation $\Gamma$, for the skeleton variables +# +# These are given by + +D = num_cell_dims(model) +Ω = Triangulation(ReferenceFE{D}, model) +Γ = Triangulation(ReferenceFE{D-1}, model) + +# ## FESpaces +# +# HHO uses two different finite element spaces: +# 1. A scalar space for the bulk variable uT (V) +# 2. A scalar space for the skeleton variable uF (M) +# +# We then define discontinuous finite element spaces of the approppriate order, locally $$\mathbb{P}^k$$. +# Because we are using a mixed-order scheme, the bulk space has a higher polynomial order +# than the skeleton space. +# Note that only the skeletal space has Dirichlet boundary conditions. + +order = 1 +V = FESpaces.PolytopalFESpace(Ω, Float64, order+1; space=:P) # Bulk space +M = FESpaces.PolytopalFESpace(Γ, Float64, order; space=:P, dirichlet_tags="boundary") # Skeleton space +N = TrialFESpace(M,u) + +# ## MultiField Structure +# +# Since we are doing static condensation, we need assemble by blocks. In particular, the +# `StaticCondensationOperator` expects the variables to be groupped in two blocks: +# - The eliminated variables (in this case, the volume variables q and u) +# - The retained variables (in this case, the interface variable m) +# We will assemble by blocks using the `BlockMultiFieldStyle` API. + +mfs = MultiField.BlockMultiFieldStyle(2,(1,1)) +X = MultiFieldFESpace([V, N];style=mfs) +Y = MultiFieldFESpace([V, M];style=mfs) + +# ## PatchTopology and PatchTriangulation +# +# A key aspect of hybrid methods is the use of static condensation, which is the +# elimination of cell unknowns to reduce the size of the global system. +# To achieve this, we need to be able to assemble and solve local problems on each cell, that +# involve +# - contributions from the cell itself +# - contributions from the cell faces +# To this end, Gridap provides a general framework for patch-assembly and solves. The idea +# is to define a patch decomposition of the mesh (in this case, a patch is a cell and its sourrounding +# faces). We can then gather contributions for each patch, solve the local problems, and +# assemble the results into the global system. +# +# The following code creates the required `PatchTopology` for the problem at hand. We then +# take d-dimensional slices of it by the means of `PatchTriangulation` and `PatchBoundaryTriangulation`. +# These are the `Triangulation`s we will integrate our weakform over. + +ptopo = Geometry.PatchTopology(model) +Ωp = Geometry.PatchTriangulation(model,ptopo) +Γp = Geometry.PatchBoundaryTriangulation(model,ptopo) + +qdegree = 2*(order+1) +dΩp = Measure(Ωp,qdegree) +dΓp = Measure(Γp,qdegree) + +# ## Local operators +# +# A key feature of HHO is the use of local solves to define local projections of our bulk and +# skeleton variables. Just like for static condensation, we will use patch assembly to +# gather the contributions from each patch and solve the local problems. +# The result is then an element of the space we are projecting to, given as a +# linear combination of the basis of that space. This whole abstraction is taken +# care of by the `LocalOperator` object. +# +# For the mixed-order Poisson problem, we require two local projections. +# +# ### L2 projection operator +# +# First, an L2 local projection operator onto the mesh faces. This is the simplest +# operator we can define, such that given $u$ in some undefined space, we find $Pu \in V$ st +# +# $$(Pu,q) = (u,q) \forall q \in V$$ +# +# This signature for the `LocalOperator` assumes that the rhs and lhs for each local problem +# are given by a single cell contribution (no patch assembly required). +# Note also the use of `FESpaceWithoutBCs`, which strips the boundary conditions from the +# space `V`. This is because we do not want to take into account boundary conditions +# when projecting onto the space. +# The local solve map is given by `LocalSolveMap`, which by default uses Julia's LU +# factorization to solve the local problem exactly. + +function projection_operator(V, Ω, dΩ) + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + mass(u,v) = ∫(u⋅Π(v,Ω))dΩ + V0 = FESpaces.FESpaceWithoutBCs(V) + P = LocalOperator( + LocalSolveMap(), V0, mass, mass; trian_out = Ω + ) + return P +end + +# ### Reconstruction operator +# +# Finally, we build the so-called reconstruction operator. This operator is highly tied to the +# ellipic projector, and projects our bulk-skeleton variable pair onto a bulk space of higher order. +# +# It's definition can be found in HHO literature, and requires solving a constrained +# local problem on each cell and its faces. We therefore use patch assembly, and impose +# our constraint using an additional space `Λ` as a local Lagrange multiplier. +# Naturally, we want to eliminate the multiplier and return a solution in the reconstructed +# space `L`. This is taken care of by the `LocalPenaltySolveMap`, and the kwarg `space_out = L` +# which overrides the default behavior of returning a solution in the test space `Y`. + +function reconstruction_operator(ptopo,order,X,Ω,Γp,dΩp,dΓp) + L = FESpaces.PolytopalFESpace(Ω, Float64, order+1; space=:P) + Λ = FESpaces.PolytopalFESpace(Ω, Float64, 0; space=:P) + + n = get_normal_vector(Γp) + Πn(v) = ∇(v)⋅n + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + lhs((u,λ),(v,μ)) = ∫( (∇(u)⋅∇(v)) + (μ*u) + (λ*v) )dΩp + rhs((uT,uF),(v,μ)) = ∫( (∇(uT)⋅∇(v)) + (uT*μ) )dΩp + ∫( (uF - Π(uT,Γp))*(Πn(v)) )dΓp + + Y = FESpaces.FESpaceWithoutBCs(X) + W = MultiFieldFESpace([L,Λ];style=BlockMultiFieldStyle()) + R = LocalOperator( + LocalPenaltySolveMap(), ptopo, W, Y, lhs, rhs; space_out = L + ) + return R +end + +PΓ = projection_operator(M, Γp, dΓp) +R = reconstruction_operator(ptopo,order,Y,Ωp,Γp,dΩp,dΓp) + +# ## Weakform and assembly +# +# We can now define: +# - The consistency term `a` +# - The stabilization term `s` +# - The rhs term `l` + +hTinv = CellField(1 ./ collect(get_array(∫(1)dΩp)), Ωp) + +function a(u,v) + Ru_Ω, Ru_Γ = R(u) + Rv_Ω, Rv_Γ = R(v) + return ∫(∇(Ru_Ω)⋅∇(Rv_Ω) + ∇(Ru_Γ)⋅∇(Rv_Ω) + ∇(Ru_Ω)⋅∇(Rv_Γ) + ∇(Ru_Γ)⋅∇(Rv_Γ))dΩp +end + +function s(u,v) + function SΓ(u) + u_Ω, u_Γ = u + return PΓ(u_Ω) - u_Γ + end + return ∫(hTinv * (SΓ(u)⋅SΓ(v)))dΓp +end + +l((vΩ,vΓ)) = ∫(f⋅vΩ)dΩp + +# ### Patch-FESpaces +# +# An additional difficulty in HHO methods is that our reconstructed functions `R(u)` are +# hard to assemble. They are defined on the cells, but depend on skeleton degrees of freedom. +# We therefore cannot assemble them using the original FESpace `N`. Instead, we will create \ +# view of the original space that is defined on the patches, and can be used to assemble +# the local contributions depending on `R`. It can be done by the means of `PatchFESpace`: + +Xp = FESpaces.PatchFESpace(X,ptopo) + +# ### Assembly without static condensation +# +# We can now proceed to evaluate and assemble all our contributions. Note that some of +# the contributions depend on the reconstructed variables, e.g `a(u,v)`, and need to +# be assembled using the `PatchFESpace` `Xp`, while others can be assembled using the original +# FESpace `X` (e.g. `s(u,v)` and `l(v)`). +# Assembly is therefore somewhat more complex than in the standard case. We need to +# collect the contributions for every (test,trial) pair, and then merge them into a single +# data structure that can be passed to the assembler. +# In the case of mixed-order HHO, we only have two combinations, (`Xp`, `Xp`) and (`X`, `X`), +# but the cross-terms appear also in the case of the original HHO formulation. + +global_assem = SparseMatrixAssembler(X,Y) + +function weakform() + u, v = get_trial_fe_basis(X), get_fe_basis(Y) + data = FESpaces.collect_and_merge_cell_matrix_and_vector( + (Xp, Xp, a(u,v), DomainContribution(), zero(Xp)), + (X, Y, s(u,v), l(v), zero(X)) + ) + assemble_matrix_and_vector(global_assem,data) +end + +A, b = weakform() +x = A \ b + +ui, ub = FEFunction(X,x) +eu = ui - u +l2u = sqrt(sum( ∫(eu * eu)dΩp)) +h1u = l2u + sqrt(sum( ∫(∇(eu) ⋅ ∇(eu))dΩp)) + +# ### Assembly with static condensation + +patch_assem = FESpaces.PatchAssembler(ptopo,X,Y) + +function patch_weakform() + u, v = get_trial_fe_basis(X), get_fe_basis(Y) + data = FESpaces.collect_and_merge_cell_matrix_and_vector(patch_assem, + (Xp, Xp, a(u,v), DomainContribution(), zero(Xp)), + (X, Y, s(u,v), l(v), zero(X)) + ) + return assemble_matrix_and_vector(patch_assem,data) +end + +op = MultiField.StaticCondensationOperator(X,patch_assem,patch_weakform()) +ui, ub = solve(op) + +eu = ui - u +l2u = sqrt(sum( ∫(eu * eu)dΩp)) +h1u = l2u + sqrt(sum( ∫(∇(eu) ⋅ ∇(eu))dΩp)) + +# ## Going further +# +# This tutorial has introduced the basic concepts of HHO methods using the simplest their +# simplest form, e.g. mixed-order HHO for the Poisson equation. +# More advanced drivers can be found with Gridap's tests. In particular: +# +# - [Poisson with original HHO formulation](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HHOPolytopalTests.jl) +# - [Incompressible Stokes](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HHOMixedStokesPolytopal.jl) +# - [Linear Elasticity](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HHOMixedElasticity.jl) +# diff --git a/src/stokes_blocks.jl b/src/stokes_blocks.jl new file mode 100644 index 0000000..9327e57 --- /dev/null +++ b/src/stokes_blocks.jl @@ -0,0 +1,220 @@ +# # Incompressible Stokes equations in a 2D/3D cavity +# +# This example solves the incompressible Stokes equations, given by +# +# ```math +# \left\lbrace +# \begin{aligned} +# -\Delta u - \nabla p &= f \quad \text{in} \quad \Omega, \\ +# \nabla \cdot u &= 0 \quad \text{in} \quad \Omega, \\ +# u &= \hat{x} \quad \text{in} \quad \Gamma_\text{top} \subset \partial \Omega, \\ +# u &= 0 \quad \text{in} \quad \partial \Omega \backslash \Gamma_\text{top} \\ +# \end{aligned} +# \right. +# ``` +# +# where $\Omega = [0,1]^d$. +# +# We use a mixed finite-element scheme, with $Q_k \times P_{k-1}^{-}$ elements for the velocity-pressure pair. +# +# To solve the linear system, we use a FGMRES solver preconditioned by a block-diagonal or +# block-triangular Shur-complement-based preconditioner. +# +# +# ## Block structure +# +# The discretized system has a natural 2×2 block structure: +# +# ```math +# \begin{bmatrix} +# A & B^T \\ +# B & 0 +# \end{bmatrix} +# \begin{bmatrix} +# u \\ +# p +# \end{bmatrix} = +# \begin{bmatrix} +# f \\ +# 0 +# \end{bmatrix} +# ``` +# +# where: +# - $A$: Vector Laplacian (velocity block) +# - $B$: Divergence operator +# - $B^T$: Gradient operator +# +# ## Solution strategy +# +# We use a FGMRES solver preconditioned by an upper block-triangular preconditioner: +# +# ```math +# P = \begin{bmatrix} +# A & B^T \\ +# 0 & -\hat{S} +# \end{bmatrix} +# ``` +# +# where $\hat{S}$ is an approximation of the Schur complement $S = BA^{-1}B^T$, which we +# will approximate using a pressure mass matrix. + +using LinearAlgebra +using BlockArrays + +using Gridap +using Gridap.MultiField + +using GridapSolvers +using GridapSolvers.LinearSolvers +using GridapSolvers.BlockSolvers: LinearSystemBlock, BiformBlock +using GridapSolvers.BlockSolvers: BlockDiagonalSolver, BlockTriangularSolver + +# ## Geometry and FESpaces +# +# See the basic Stokes tutorial for a detailed explanation. +# + +nc = (8,8) +Dc = length(nc) +domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) + +model = CartesianDiscreteModel(domain,nc) +labels = get_face_labeling(model) +if Dc == 2 + add_tag_from_tags!(labels,"top",[6]) + add_tag_from_tags!(labels,"walls",[1,2,3,4,5,7,8]) +else + add_tag_from_tags!(labels,"top",[22]) + add_tag_from_tags!(labels,"walls",[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25,26]) +end + +order = 2 +qdegree = 2*(order+1) +reffe_u = ReferenceFE(lagrangian,VectorValue{Dc,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +u_walls = (Dc==2) ? VectorValue(0.0,0.0) : VectorValue(0.0,0.0,0.0) +u_top = (Dc==2) ? VectorValue(1.0,0.0) : VectorValue(1.0,0.0,0.0) + +V = TestFESpace(model,reffe_u,dirichlet_tags=["walls","top"]); +U = TrialFESpace(V,[u_walls,u_top]); +Q = TestFESpace(model,reffe_p;conformity=:L2,constraint=:zeromean) + +# ## Block multi-field spaces +# +# Our first difference will come from how we define our multi-field spaces: +# Because we want to be able to use the block-structure of the linear system, +# we have to assemble our problem by blocks. The block structure of the resulting +# linear system is determined by the `BlockMultiFieldStyle` we use to define the +# multi-field spaces. +# +# A `BlockMultiFieldStyle` takes three arguments: +# - `N`: The number of blocks we want +# - `S`: An N-tuple with the number of fields in each block +# - `P`: A permutation of the fields in the multi-field space, which determines +# how fields are grouped into blocks. +# +# By default, we create as many blocks as there are fields in the multi-field space, +# and each block contains a single field with no permutation. + +mfs = BlockMultiFieldStyle(2,(1,1),(1,2)) +X = MultiFieldFESpace([U,Q];style=mfs) +Y = MultiFieldFESpace([V,Q];style=mfs) + +# ## Weak form and integration + +Ω = Triangulation(model) +dΩ = Measure(Ω,qdegree) + +α = 1.e1 +f = (Dc==2) ? VectorValue(1.0,1.0) : VectorValue(1.0,1.0,1.0) +a((u,p),(v,q)) = ∫(∇(v)⊙∇(u))dΩ - ∫(divergence(v)*p)dΩ - ∫(divergence(u)*q)dΩ +l((v,q)) = ∫(v⋅f)dΩ + +op = AffineFEOperator(a,l,X,Y) + +# ### Block structure of the linear system +# +# As per usual, we can extract the matrix and vector of the linear system from the operator. +# Notice now that unlike in previous examples, the matrix is a `BlockMatrix` type, from +# the `BlockArrays.jl` package, which allows us to work with block-structured matrices. + +A = get_matrix(op) +b = get_vector(op) + +# ## Block solvers +# +# We will now setup two types of block preconditioners for the Stokes system. In both cases, +# we will use the preconditioners to solve the linear system using a Flexible GMRES solver. +# The idea behind these preconditioners is the well-known property that the Schur complement +# of the velocity block can be well approximated by a scaled pressure mass matrix. Moreover, +# in our pressure discretization is discontinuous which means that the pressure mass matrix +# is block-diagonal and easily invertible. +# In this example, we will use an exact LU solver for the velocity block and a CG solver +# with Jacobi preconditioner for the pressure block. + +# ### Block diagonal preconditioner +# +# The simplest block preconditioner is the block-diagonal preconditioner. +# The only ingredients required are +# - the sub-solvers for each diagonal block and +# - the diagonal blocks we want to use. +# +# The sub-solvers are defined as follows: + +u_solver = LUSolver() +p_solver = CGSolver(JacobiLinearSolver();maxiter=20,atol=1e-14,rtol=1.e-6) + +# The block structure is defined using the block API. We provide different types of blocks, +# that might have different uses depending on the problem at hand. We will here use two of +# the most common block types: +# - The `LinearSystemBlock` defines a block that is taken directly (and aliased) form +# the linear system matrix `A`. We will use this for the velocity block. +# - For the pressure block, however, the pressure mass matrix is not directly available +# in the system matrix. Instead, we will have to integrate it using the Gridap API, as usual. +# These abstract concept is implemented in the `BiformBlock` type, which allows the +# user to define a block from a bilinear form. +# All in all, we define the block structure as follows: + +u_block = LinearSystemBlock() +p_block = BiformBlock((p,q) -> ∫(-(1.0/α)*p*q)dΩ,Q,Q) + +# With these ingredients, we can now define the block diagonal preconditioner as follows: + +PD = BlockDiagonalSolver([u_block,p_block],[u_solver,p_solver]) +solver_PD = FGMRESSolver(20,PD;atol=1e-10,rtol=1.e-12,verbose=true) + +uh, ph = solve(solver_PD, op) + +# ### Block upper-triangular preconditioner +# +# A slighly more elaborate preconditioner (but also more robust) is the +# block upper-triangular preconditioner. The ingredients are similar: +# - the sub-solvers for each diagonal block +# - the blocks that define the block structure, now including the off-diagonal blocks +# - the coefficients for the off-diagonal blocks, where zero coefficients +# indicate that the block is not used. +# +# We will also represent the off-diagonal blocks using the `LinearSystemBlock` type. + +sblocks = [ u_block LinearSystemBlock(); + LinearSystemBlock() p_block ] +coeffs = [1.0 1.0; + 0.0 1.0] +PU = BlockTriangularSolver(sblocks,[u_solver,p_solver],coeffs,:upper) +solver_PU = FGMRESSolver(20,PU;atol=1e-10,rtol=1.e-12,verbose=true) + +uh, ph = solve(solver_PU, op) + +# As you can see, the block upper-triangular preconditioner is quite better +# than the block diagonal one. + +# ## Going further +# +# If you want to see more examples of how to use the block solvers, +# you can check the documentation in [GridapSolvers.jl](https://gridap.github.io/GridapSolvers.jl/stable/), +# as well as it's `test/Applications` folder. +# +# There you will find more complicated examples, such as using a GMG solver to +# solve the velocity block. diff --git a/test/runtests.jl b/test/runtests.jl index 7456f93..1ad6452 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,7 @@ function get_title(files,filename) error("File $filename not found!") end -if (length(ARGS) != 0) +if !iszero(length(ARGS)) files = [get_title(Tutorials.files,filename)=>filename for filename in ARGS] else files = Tutorials.files