Updated GMG drivers

src/Applications/hunt.jl

@@ -54,6 +54,8 @@ function _hunt(;
   σw1=0.1,
   σw2=10.0,
   tw=0.0,
+  order = 2,
+  order_j = order,
   nsums = 10,
   vtk=true,
   title = "test",
@@ -63,12 +65,14 @@ function _hunt(;
   jac_assemble = false,
   solve = true,
   solver = :julia,
+  formulation = :mhd,
+  rt_scaling = false,
   verbose = true,
   BL_adapted = true,
   kmap_x = 1,
   kmap_y = 1,
   ranks_per_level = nothing
-  )
+)
 
   info = Dict{Symbol,Any}()
   params = Dict{Symbol,Any}(
@@ -105,12 +109,21 @@ function _hunt(;
   N = Ha^2/Re
   f̄ = (L/(ρ*u0^2))*VectorValue(f)
   B̄ = (1/B0)*VectorValue(B)
-  α = 1.0
-  β = 1.0/Re
-  γ = N
   σ̄1 = σw1/σ
   σ̄2 = σw2/σ
 
+  if formulation == :cfd # Option 1 (CFD)
+    α = 1.0
+    β = 1.0/Re
+    γ = N
+  elseif formulation == :mhd # Option 2 (MHD) is chosen in the experimental article
+    α = (1.0/N)
+    β = (1.0/Ha^2)
+    γ = 1.0
+  else
+    error("Unknown formulation")
+  end
+
   # DiscreteModel in terms of reduced quantities
 
   model = hunt_mesh(parts,params,nc,rank_partition,L,tw,Ha,kmap_x,kmap_y,BL_adapted,ranks_per_level)
@@ -119,30 +132,28 @@ function _hunt(;
     writevtk(model,"hunt_model")
   end
 
+  params[:fluid] = Dict{Symbol,Any}(
+    :domain=>nothing,
+    :α=>α,
+    :β=>β,
+    :γ=>γ,
+    :f=>f̄,
+    :B=>B̄,
+    :ζ=>ζ,
+  )
+
   if tw > 0.0
     σ_Ω = solid_conductivity(σ̄1,σ̄2,Ω,get_cell_gids(model),get_face_labeling(model))
     params[:solid] = Dict(:domain=>"solid",:σ=>σ_Ω)
-    params[:fluid] = Dict(
-      :domain=>"fluid",
-      :α=>α,
-      :β=>β,
-      :γ=>γ,
-      :B=>B̄,
-      :f=>f̄,
-      :ζ=>ζ,
-     )
-  else
-    params[:fluid] = Dict(
-      :domain=>nothing,
-      :α=>α,
-      :β=>β,
-      :γ=>γ,
-      :f=>f̄,
-      :B=>B̄,
-      :ζ=>ζ,
-    )
+    params[:fluid][:domain] = "fluid"
   end
 
+  params[:fespaces] = Dict{Symbol,Any}(
+    :order_u => order,
+    :order_j => order_j,
+    :rt_scaling => rt_scaling ? 1.0/get_mesh_size(model) : nothing
+  )
+
   # Boundary conditions
 
   params[:bcs] = Dict(
@@ -255,13 +266,13 @@ end
 function hunt_mesh(
   parts,params,
   nc::Tuple,np::Tuple,L::Real,tw::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool,
-  ranks_per_level)
+  ranks_per_level
+)
   if isnothing(ranks_per_level) # Single grid
     model = Meshers.hunt_generate_base_mesh(parts,np,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted)
     params[:model] = model
   else # Multigrid
-    base_model = Meshers.hunt_generate_base_mesh(nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted)
-    mh = Meshers.generate_mesh_hierarchy(parts,base_model,0,ranks_per_level)
+    mh = Meshers.hunt_generate_mesh_hierarchy(parts,ranks_per_level,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted)
     params[:multigrid] = Dict{Symbol,Any}(
       :mh => mh,
       :num_refs_coarse => 0,

src/Fixes.jl

@@ -8,4 +8,4 @@
 function Geometry.get_polytopes(model::GridapDistributed.DistributedDiscreteModel)
   polys = map(get_polytopes,local_views(model))
   return PartitionedArrays.getany(polys)
-end
+end

src/Main.jl

@@ -269,7 +269,7 @@ function _fe_space(::Val{:j},params)
   uses_mg = space_uses_multigrid(params[:solver])[3]
   model = uses_mg ? params[:multigrid][:mh] : params[:model]
 
-  phi = rt_scaling(model,params[:fespaces])
+  phi = rt_scaling(params[:model],params[:fespaces])
   reffe_j = ReferenceFE(raviart_thomas,Float64,k-1;basis_type=:jacobi,phi=phi)
   params[:fespaces][:reffe_j] = reffe_j
 
@@ -356,10 +356,11 @@ const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.Distribu
 
 _interior(model,domain::DiscreteModelTypes) = Interior(domain)
 _interior(model,domain::TriangulationTypes) = domain
-_interior(model,domain::Nothing) = Triangulation(model) # This should be removed, but Gridap needs fixes
+_interior(model,domain::Nothing) = Triangulation(model)
 _interior(model,domain) = Interior(model,tags=domain)
 
 _boundary(model,domain::TriangulationTypes) = domain
+_boundary(model,domain::Nothing) = Boundary(model)
 _boundary(model,domain) = Boundary(model,tags=domain)
 
 _skeleton(model,domain::TriangulationTypes) = SkeletonTriangulation(domain)
@@ -375,18 +376,25 @@ function _setup_trians!(params)
   solid = params[:solid]
   Ωs = !isnothing(solid) ? _interior(params[:model],solid[:domain]) : nothing
 
-  if !uses_multigrid(params[:solver])
+  #if !uses_multigrid(params[:solver])
     params[:Ω]  = Ω
     params[:Ωf] = Ωf
     params[:Ωs] = Ωs
-  else
-    params[:multigrid][:Ω]  = Ω
-    params[:multigrid][:Ωf] = Ωf
-    params[:multigrid][:Ωs] = Ωs
-    params[:Ω]  = Ω[1]
-    params[:Ωf] = Ωf[1]
-    params[:Ωs] = Ωs[1]
+  #else
+  #  params[:multigrid][:Ω]  = Ω
+  #  params[:multigrid][:Ωf] = Ωf
+  #  params[:multigrid][:Ωs] = Ωs
+  #  params[:Ω]  = Ω[1]
+  #  params[:Ωf] = Ωf[1]
+  #  params[:Ωs] = Ωs[1]
+  #end
+
+  if uses_multigrid(params[:solver])
+    params[:multigrid][:Ω]  = params[:multigrid][:mh]
+    params[:multigrid][:Ωf] = params[:multigrid][:mh]
+    params[:multigrid][:Ωs] = params[:multigrid][:mh]
   end
+
 end
 
 # Random vector generation
@@ -441,21 +449,29 @@ function _allocate_solution(op::TransientFEOperator,args...)
   nothing
 end
 
-function _get_cell_size(t::Triangulation)
+# Mesh sizes
+
+get_cell_size(t::TriangulationTypes) = CellField(_get_cell_size(t),t)
+
+get_cell_size(m::DiscreteModelTypes) = get_cell_size(Triangulation(m))
+_get_cell_size(m::DiscreteModelTypes) = _get_cell_size(Triangulation(m))
+
+function _get_cell_size(t::Triangulation) :: Vector{Float64}
+  if iszero(num_cells(t))
+    return Float64[]
+  end
   meas = get_cell_measure(t)
   d = num_dims(t)
-  map(m->m^(1/d),meas)
+  return collect(Float64, meas .^ (1/d))
 end
 
 function _get_cell_size(t::GridapDistributed.DistributedTriangulation)
   map(_get_cell_size,local_views(t))
 end
 
-_get_cell_size(m::DiscreteModelTypes) = _get_cell_size(Triangulation(m))
-
-_get_mesh_size(m::DiscreteModel) = minimum(_get_cell_size(m))
+get_mesh_size(m::DiscreteModel) = minimum(_get_cell_size(m))
 
-function _get_mesh_size(m::GridapDistributed.DistributedDiscreteModel)
-  h = map(_get_mesh_size,local_views(m))
+function get_mesh_size(m::GridapDistributed.DistributedDiscreteModel)
+  h = map(get_mesh_size,local_views(m))
   return reduce(min,h;init=one(eltype(h)))
 end

src/Meshers/hunt_mesher.jl

@@ -53,6 +53,10 @@ end
 
 function hunt_add_tags!(model,L::Real,tw::Real)
   labels = get_face_labeling(model)
+  hunt_add_tags!(labels,model,L,tw)
+end
+
+function hunt_add_tags!(labels,model,L::Real,tw::Real)
   if tw > 0.0 ## add solid tags
     # When model is part of a distributed model we are using
     # that the maximum entity (9) is the same in all parts, 
@@ -128,3 +132,18 @@ function hunt_generate_base_mesh(
   hunt_add_tags!(model,L,tw)
   return model
 end
+
+function hunt_generate_mesh_hierarchy(
+  parts,np_per_level,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted
+)
+  Lt = L+tw
+  _nc = (nc[1],nc[2],3)
+  domain = (-1.0,1.0,-1.0,1.0,0.0,0.1)
+  CartesianModelHierarchy(
+    parts,np_per_level,domain,_nc;
+    nrefs = (2,2,1),
+    isperiodic = (false,false,true),
+    map = hunt_stretch_map(Lt,Ha,kmap_x,kmap_y,BL_adapted),
+    add_labels! = labels -> hunt_add_tags!(labels,nothing,L,tw) # TODO: will not work for solid
+  )
+end