Alternative discretizations

.github/workflows/ci.yml

@@ -11,7 +11,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.10'
+          - '1.11'
         os:
           - ubuntu-latest
         arch:

.github/workflows/ci_mpi.yml

@@ -7,11 +7,12 @@ jobs:
     env:
       P4EST_ROOT_DIR: "/opt/p4est/2.8.5/"
       JULIA_PETSC_LIBRARY: "/opt/petsc/3.23.0/lib/libpetsc"
+      GMSHROOT: "/opt/gmsh-4.13-1-Linux64-sdk"
     strategy:
       fail-fast: false
       matrix:
         version:
-          - '1.10'
+          - '1.11'
         os:
           - ubuntu-latest
         arch:
@@ -21,6 +22,20 @@ jobs:
         run: |
           sudo apt-get update 
           sudo apt-get install -y wget gfortran g++ openmpi-bin libopenmpi-dev 
+      - name: Install gmsh-sdk
+        run: |
+          # Install gmsh-sdk from sources
+          CURR_DIR=$(pwd)
+          PACKAGE=gmsh
+          VERSION=4.13.1
+          INSTALL_ROOT=/opt
+          TAR_FILE=$PACKAGE-$VERSION-Linux64-sdk.tgz
+          URL="https://gmsh.info/bin/Linux"
+          ROOT_DIR=/tmp
+          wget -q $URL/$TAR_FILE -O $ROOT_DIR/$TAR_FILE
+          tar xzf $ROOT_DIR/$TAR_FILE -C $INSTALL_ROOT
+          rm -rf $ROOT_DIR/$TAR_FILE
+          cd $CURR_DIR
       - name: Install p4est
         run: |
           # Install p4est 2.2 from sources
@@ -71,7 +86,6 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v2
-
       - name: add MPIPreferences and external packages
         shell: julia --color=yes --project=. {0}
         run: |

Project.toml

@@ -25,18 +25,21 @@ Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
-gmsh_jll = "630162c2-fc9b-58b3-9910-8442a8a132e6"
 
 [compat]
 BSON = "0.3"
+BlockArrays = "1"
+DrWatson = "2"
 FileIO = "1"
 FillArrays = "1.13.0"
-Gridap = "0.19"
+FlameGraphs = "1"
+ForwardDiff = "1"
+Gridap = "0.19, 0.20"
 GridapDistributed = "0.4"
 GridapGmsh = "0.7.2"
 GridapP4est = "0.3.7"
 GridapPETSc = "0.5.0"
-GridapSolvers = "0.6"
+GridapSolvers = "0.6, 0.7"
 MPI = "0.20"
 PackageCompiler = "2"
 SparseMatricesCSR = "0.6.6"

compile/compile.jl

@@ -6,8 +6,8 @@ pkgs = Symbol[]
 append!(pkgs, [Symbol(v.name) for v in values(Pkg.dependencies()) if v.is_direct_dep],)
 
 create_sysimage(pkgs,
-  sysimage_path=joinpath(@__DIR__,"..","GridapMHD.so"),
-  precompile_execution_file=joinpath(@__DIR__,"warmup.jl"),
-  include_transitive_dependencies=false,
-  sysimage_build_args=`-O3 --check-bounds=no`
-  )
+  sysimage_path = joinpath(@__DIR__,"..","GridapMHD.so"),
+  precompile_execution_file = joinpath(@__DIR__,"warmup_h1h1.jl"),
+  include_transitive_dependencies = false,
+  sysimage_build_args = `-O3`
+)

compile/warmup_h1h1.jl

@@ -0,0 +1,42 @@
+
+using GridapPETSc
+using SparseMatricesCSR
+using SparseArrays
+
+using GridapMHD: cavity
+
+np = (1,1,1)
+cavity(
+  nc = (4,4,4),
+  np = np,
+  backend = :mpi,
+  fluid_disc = :Qk_dPkm1,
+  current_disc = :H1,
+  solver = Dict(
+    :solver => :h1h1blocks,
+    :matrix_type    => SparseMatrixCSR{0,Float64,Int},
+    :vector_type    => Vector{Float64},
+    :block_solvers  => [:petsc_mumps,:petsc_cg_jacobi,:petsc_gmres_amg],
+    :petsc_options  => "-ksp_error_if_not_converged true -ksp_converged_reason",
+  ),
+  ζᵤ = 10.0,
+  solid = true
+)
+
+cavity(
+  nc = (8,8,8),
+  np = np,
+  backend = :mpi,
+  fluid_disc = :Qk_dPkm1,
+  current_disc = :H1,
+  solver = Dict(
+    :solver => :h1h1blocks,
+    :matrix_type    => SparseMatrixCSC{Float64,Int},
+    :vector_type    => Vector{Float64},
+    :block_solvers  => [:gmg,:petsc_cg_jacobi,:petsc_gmres_amg],
+    :petsc_options  => "-ksp_error_if_not_converged true -ksp_converged_reason",
+  ),
+  ranks_per_level = [np,np],
+  ζᵤ = 10.0,
+  solid = true
+)

src/Applications/cavity.jl

@@ -41,7 +41,8 @@ function _cavity(;
   ν = 1.0,
   ρ = 1.0,
   σ = 1.0,
-  ζ = 0.0,   # Augmented Lagrangian weight  
+  ζᵤ = 0.0,   # Augmented Lagrangian weights 
+  ζⱼ = 0.0, 
   μ = 0,     # Stabilization weight
   B = VectorValue(0.0, 0.0, 10.0),
   f = VectorValue(0.0, 0.0, 0.0),
@@ -60,6 +61,7 @@ function _cavity(;
   verbose=true,
   vtk=true,
   convection=:newton,
+  solid = false,
   closed_cavity=true,
   adaptivity_method = 0,
   fluid_disc = ifelse(iszero(adaptivity_method),:Qk_dPkm1,:SV),
@@ -94,7 +96,7 @@ function _cavity(;
   params[:solver] = solver
 
   # Model
-  model = cavity_mesh(parts,params,nc,np,L,ranks_per_level,adaptivity_method)
+  model = cavity_mesh(parts,params,nc,np,L,ranks_per_level,adaptivity_method,solid)
 
   # Reduced quantities
   Re = u0 * L / ν
@@ -123,14 +125,24 @@ function _cavity(;
     :γ => γ,
     :f => f̄,
     :B => B̄,
-    :ζ => ζ,
+    :ζᵤ => ζᵤ,
+    :ζⱼ => ζⱼ,
     :convection => convection,
   )
+  if solid
+    params[:fluid][:domain] = "fluid"
+    params[:solid] = Dict(
+      :domain => "solid",
+      :σ => σ,
+      :ζⱼ => ζⱼ,
+    )
+  end
 
   # FESpaces and Boundary conditions
   uw = VectorValue(0.0, 0.0, 0.0)
   ul = VectorValue(1.0, 0.0, 0.0)
   ji = VectorValue(0.0, 0.0, 0.0)
+  params[:x0] = initial_value
   params[:fespaces] = Dict{Symbol,Any}(
     :order_u => order,
     :order_j => order_j,
@@ -139,8 +151,6 @@ function _cavity(;
     :current_disc => current_disc,
   )
 
-  params[:x0] = initial_value
-
   if closed_cavity
     params[:bcs] = Dict{Symbol,Any}(
       :u => Dict(:tags => ["cavity", "lid"], :values => [uw, ul]),
@@ -149,10 +159,13 @@ function _cavity(;
     params[:fespaces][:p_constraint] = :zeromean
   else
     params[:bcs] = Dict{Symbol,Any}(
-      :u => Dict(:tags => ["wall", "lid"], :values => [uw, ul]), # Bottom wall is Newman
+      :u => Dict(:tags => ["wall", "lid"], :values => [uw, ul]), # Bottom is Neumann
       :j => Dict(:tags => "insulating", :values => ji),
     )
   end
+  if current_disc == :H1
+    params[:bcs][:φ] = Dict(:tags => "lid", :values => 0.0)
+  end
 
   if μ > 0
     params[:bcs][:stabilization] = Dict(:μ=>μ)
@@ -173,13 +186,23 @@ function _cavity(;
   if vtk
     tic!(t, barrier=true)
     Ω = Interior(model)
-    ūh, p̄h, j̄h, φ̄h = xh
-    uh = u0 * ūh
-    ph = (ρ * u0^2) * p̄h
-    jh = (σ * u0 * B0) * j̄h
-    φh = (u0 * B0 * L) * φ̄h
-    div_jh = ∇·jh
-    div_uh = ∇·uh
+    if current_disc != :H1
+      ūh, p̄h, j̄h, φ̄h = xh
+      uh = u0 * ūh
+      ph = (ρ * u0^2) * p̄h
+      jh = (σ * u0 * B0) * j̄h
+      φh = (u0 * B0 * L) * φ̄h
+      div_jh = ∇·jh
+      div_uh = ∇·uh
+    else
+      ūh, p̄h, φ̄h = xh
+      uh = u0 * ūh
+      ph = (ρ * u0^2) * p̄h
+      φh = (u0 * B0 * L) * φ̄h
+      jh = σ * (uh × B - ∇(φh))
+      div_jh = σ*((∇×uh)⋅B - Δ(φh)) 
+      div_uh = ∇·uh
+    end
     writevtk(
       Ω, joinpath(path,title), order=max(order,order_j), 
       cellfields=["uh" => uh, "ph" => ph, "jh" => jh, "phi" => φh, "div_jh" => div_jh, "div_uh" => div_uh],
@@ -191,7 +214,7 @@ function _cavity(;
   info[:ncells]  = num_cells(model)
   info[:ndofs_u] = length(get_free_dof_values(ūh))
   info[:ndofs_p] = length(get_free_dof_values(p̄h))
-  info[:ndofs_j] = length(get_free_dof_values(j̄h))
+  info[:ndofs_j] = (current_disc != :H1) ? length(get_free_dof_values(j̄h)) : 0
   info[:ndofs_φ] = length(get_free_dof_values(φ̄h))
   info[:ndofs]   = sum([info[:ndofs_u], info[:ndofs_p], info[:ndofs_j], info[:ndofs_φ]])
   info[:Re]      = Re
@@ -200,44 +223,66 @@ function _cavity(;
   return info, t
 end
 
-function add_cavity_tags!(model::CartesianDiscreteModel)
-  labels = get_face_labeling(model)
-  add_cavity_tags!(labels)
+function add_cavity_tags!(model::GridapDistributed.DistributedDiscreteModel, tw, L)
+  map(local_views(model)) do model
+    add_cavity_tags!(model, tw, L)
+  end
 end
 
-function add_cavity_tags!(model::GridapDistributed.DistributedDiscreteModel)
-  map(add_cavity_tags!, local_views(model))
+function add_cavity_tags!(mh::MultilevelTools.ModelHierarchy, tw, L)
+  map(mh) do mhl
+    m = get_model(mhl)
+    add_cavity_tags!(m, tw, L)
+    mred = get_model_before_redist(mhl)
+    if ! isnothing(mred) && mred !== m
+      add_cavity_tags!(mred, tw, L)
+    end
+  end
 end
 
-function add_cavity_tags!(labels)
-  Γl = [22]
-  Γb = [21]
-  Γw = vcat(collect(1:20),collect(23:26))
-  add_tag_from_tags!(labels, "lid", Γl)
-  add_tag_from_tags!(labels, "bottom", Γb)
-  add_tag_from_tags!(labels, "wall", Γw)
-  add_tag_from_tags!(labels, "cavity", ["wall","bottom"])
-  add_tag_from_tags!(labels, "insulating", "boundary")
+function add_cavity_tags!(model::DiscreteModel, tw, L)
+  topo = get_grid_topology(model)
+  labels = get_face_labeling(model)
+  add_tag_from_tags!(labels, "top", [22])
+  add_tag_from_tags!(labels, "bottom", [21])
+  add_tag_from_tags!(labels, "sides", vcat(collect(1:20),collect(23:26)))
+
+  if isnothing(tw) # Only fluid
+    add_tag_from_tags!(labels, "lid", ["top"])
+    add_tag_from_tags!(labels, "cavity", ["sides","bottom"])
+    add_tag_from_tags!(labels, "insulating", "boundary")
+  else # Solid and fluid
+    # The solid is the first layer of cells touching either the bottom or the sides
+    solid_mask(x) = (x[1] < tw[1]) || (x[1] > L - tw[1]) || (x[2] < tw[2]) || (x[2] > L - tw[2]) || (x[3] < tw[3])
+    cell_to_issolid = lazy_map(solid_mask,lazy_map(mean, get_cell_coordinates(model)))
+    cell_to_color = Gridap.Arrays.collect1d(map(x -> ifelse(x,1,2), cell_to_issolid))
+    merge!(labels,Geometry.face_labeling_from_cell_tags(topo,cell_to_color,["solid","fluid"]))
+
+    Geometry.add_tag_from_tags_intersection!(labels,"cavity",["solid","fluid"])
+    Geometry.add_tag_from_tags_setdiff!(labels,"lid",["top"],["solid"])
+    add_tag_from_tags!(labels, "insulating", "boundary")
+  end
 end
 
-function cavity_mesh(parts,params,nc::Int,np,L,ranks_per_level,simplexify)
-  return cavity_mesh(parts,params,(nc,nc,nc),np,L,ranks_per_level,simplexify)
+function cavity_mesh(parts,params,nc::Int,np,L,ranks_per_level,adaptivity_method,solid)
+  return cavity_mesh(parts,params,(nc,nc,nc),np,L,ranks_per_level,adaptivity_method,solid)
 end
 
-function cavity_mesh(parts,params,nc::Tuple,np::Int,L,ranks_per_level,simplexify)
-  return cavity_mesh(parts,params,nc,(np,1,1),L,ranks_per_level,simplexify)
+function cavity_mesh(parts,params,nc::Tuple,np::Int,L,ranks_per_level,adaptivity_method,solid)
+  return cavity_mesh(parts,params,nc,(np,1,1),L,ranks_per_level,adaptivity_method,solid)
 end
 
-function cavity_mesh(parts,params,nc::Tuple,np::Tuple,L,ranks_per_level,adaptivity_method)
+function cavity_mesh(parts,params,nc::Tuple,np::Tuple,L,ranks_per_level,adaptivity_method,solid)
   domain = (0.0,L,0.0,L,0.0,L)
+  tw = solid ? L ./ nc : nothing
   if isnothing(ranks_per_level) # Single grid
     model = CartesianDiscreteModel(parts,np,domain,nc)
-    add_cavity_tags!(model)
+    add_cavity_tags!(model,tw,L)
   else # Multigrid
     mh = CartesianModelHierarchy(
-      parts,np_per_level,domain,nc;
-      add_labels! = add_cavity_tags!,
+      parts,ranks_per_level,domain,nc
     )
+    add_cavity_tags!(mh,tw,L)
     params[:multigrid] = Dict{Symbol,Any}(
       :mh => mh,
       :num_refs_coarse => 0,

src/Applications/expansion.jl

@@ -54,7 +54,8 @@ function _expansion(;
   Ha = 1.0,
   cw = 0.028,
   τ  = 100,
-  ζ  = 0.0,   # Augmented Lagrangian weight  
+  ζᵤ=0.0, # Augmented Lagrangian weights  
+  ζⱼ=0.0, # 
   μ  = 0,     # Stabilization weight
   order = 2,
   order_j = order,
@@ -136,7 +137,8 @@ function _expansion(;
     :γ => γ,
     :f => VectorValue(0.0,0.0,0.0),
     :B => VectorValue(0.0,1.0,0.0),
-    :ζ => ζ,
+    :ζᵤ => ζᵤ,
+    :ζⱼ => ζⱼ,
     :convection => convection,
   )