Added Picard iteration

src/Applications/cavity.jl

@@ -59,11 +59,12 @@ function _cavity(;
   ranks_per_level=nothing,
   verbose=true,
   vtk=true,
-  convection=true,
+  convection=:newton,
   closed_cavity=true
 )
   @assert formulation ∈ [:cfd,:mhd]
   @assert initial_value ∈ [:zero,:solve]
+  @assert convection ∈ [:newton,:picard,:none]
 
   info = Dict{Symbol,Any}()
   params = Dict{Symbol,Any}(

src/Applications/channel.jl

@@ -66,14 +66,15 @@ function _channel(;
   γB = 0.45,
   bl_orders=(1,1,1),
   initial_value = :zero,
-  convection = true,
+  convection = :newton,
   np_per_level = nothing,
   rt_scaling = false,
   formulation = :mhd
 )
   @assert inlet ∈ [:parabolic,:shercliff,:constant]
   @assert initial_value ∈ [:zero,:inlet,:solve]
   @assert formulation ∈ [:cfd,:mhd]
+  @assert convection ∈ [:newton,:picard,:none]
 
   info = Dict{Symbol,Any}()
   params = Dict{Symbol,Any}(

src/Applications/expansion.jl

@@ -62,7 +62,7 @@ function _expansion(;
   initial_value = :zero,
   solid_coupling = :none,
   niter = nothing,
-  convection = :true,
+  convection = :newton,
   rt_scaling = false,
   savelines = false,
   petsc_options = "",
@@ -71,6 +71,7 @@ function _expansion(;
   @assert inlet ∈ [:parabolic,:shercliff,:constant]
   @assert initial_value ∈ [:zero,:inlet,:solve]
   @assert formulation ∈ [:cfd,:mhd]
+  @assert convection ∈ [:newton,:picard,:none]
 
   info   = Dict{Symbol,Any}()
   params = Dict{Symbol,Any}(

src/Main.jl

@@ -440,12 +440,14 @@ initial_guess(::Val{:zero},trial,op,params) = zero(trial)
 
 function initial_guess(::Val{:solve},trial,op,params)
   @notimplementedif isa(op,TransientFEOperator)
-  @assert params[:fluid][:convection] "Convection must be enabled to use initial guess :solve"
-  params[:fluid][:convection] = false
+  @assert has_convection(params) "Convection must be enabled to use initial guess :solve"
+
+  convection = params[:fluid][:convection]
+  params[:fluid][:convection] = :none
   xh = initial_guess(:zero,trial,op,params)
   solver = _solver(op,params)
   xh, cache = _solve(xh,solver,op,params)
-  params[:fluid][:convection] = true
+  params[:fluid][:convection] = convection
   return xh
 end
 

src/Solvers/gmg.jl

@@ -17,8 +17,8 @@ function gmg_solver(::Val{(:u,:j)},params)
   nlevs = num_levels(mh)
   k = params[:fespaces][:k]
   qdegree = map(lev -> 2*k+1,1:nlevs)
-  is_nonlinear = params[:fluid][:convection]
-
+  is_nonlinear = has_convection(params)
+  
   reffe_p = params[:fespaces][:reffe_p]
   Πp = MultilevelTools.LocalProjectionMap(divergence,reffe_p,2*k)
 

src/parameters.jl

@@ -368,22 +368,31 @@ Valid keys for `params[:fluid]` are the following.
 """
 function params_fluid(params::Dict{Symbol,Any})
   mandatory = Dict(
-   :domain=>true,
-   :α=>true,
-   :β=>true,
-   :γ=>true,
-   :B=>true,
-   :f=>false,
-   :σ=>false,
-   :ζ=>false,
-   :g=>false,
-   :convection=>false,
+   :domain => true,
+   :α => true,
+   :β => true,
+   :γ => true,
+   :B => true,
+   :f => false,
+   :σ => false,
+   :ζ => false,
+   :g => false,
+   :convection => false,
+  )
+  optional = Dict(
+    :σ => 1.0,
+    :f => VectorValue(0,0,0),
+    :ζ => 0.0,
+    :g => VectorValue(0,0,0),
+    :convection => :newton
   )
-  optional = Dict(:σ=>1.0,:f=>VectorValue(0,0,0),:ζ=>0.0,:g=>VectorValue(0,0,0),:convection=>true)
   fluid = _check_mandatory_and_add_optional(params[:fluid],mandatory,optional,params,"[:fluid]")
-  fluid
+  @assert fluid[:convection] in (:none,:newton,:picard)
+  return fluid
 end
 
+has_convection(params) = haskey(params[:fluid],:convection) && params[:fluid][:convection] != :none
+
 """
 Valid keys for `params[:solid]` are the following
 

src/weakforms.jl

@@ -68,20 +68,24 @@ function _weak_form(params,k)
   end
 
   function dc(x,dx,dy)
-    r = dc_mhd(x,dx,dy,α,dΩf)
+    if fluid[:convection] == :picard
+      r = p_dc_mhd(x,dx,dy,α,dΩf)
+    else
+      r = n_dc_mhd(x,dx,dy,α,dΩf)
+    end
     r
   end
 
   function res(x,dy)
     r = a(x,dy) - ℓ(dy)
-    if fluid[:convection]
+    if has_convection(params)
       r = r + c(x,dy)
     end
     r
   end
   function jac(x,dx,dy)
     r = a(dx,dy)
-    if fluid[:convection]
+    if has_convection(params)
       r = r + dc(x,dx,dy)
     end
     r
@@ -145,7 +149,11 @@ function _ode_weak_form(params,k)
   end
 
   function dc(x,dx,dy)
-    r = dc_mhd(x,dx,dy,α,dΩf)
+    if fluid[:convection] == :picard
+      r = p_dc_mhd(x,dx,dy,α,dΩf)
+    else
+      r = n_dc_mhd(x,dx,dy,α,dΩf)
+    end
     r
   end
 
@@ -285,20 +293,34 @@ function ℓ_mhd(dy,f,dΩ)
 end
 ℓ_mhd_u(v_u,f,dΩ) = ∫( v_u⋅f )*dΩ
 
+# Convection
+
 function c_mhd(x,dy,α,dΩ)
   u, p, j, φ = x
   v_u, v_p, v_j, v_φ = dy
   c_mhd_u_u(u,v_u,α,dΩ)
 end
 c_mhd_u_u(u,v_u,α,dΩ) = ∫( α*v_u⋅(conv∘(u,∇(u))) ) * dΩ
 
-function dc_mhd(x,dx,dy,α,dΩ)
+dc_mhd(x,dx,dy,α,dΩ) = n_dc_mhd(x,dx,dy,α,dΩ) # Default to full Newton iteration
+
+# Convection derivative: Newton iteration
+function n_dc_mhd(x,dx,dy,α,dΩ)
+  u, p, j, φ = x
+  du , dp , dj , dφ  = dx
+  v_u, v_p, v_j, v_φ = dy
+  n_dc_mhd_u_u(u,du,v_u,α,dΩ)
+end
+n_dc_mhd_u_u(u,du,v_u,α,dΩ) = ∫( α*v_u⋅( dconv∘(du,∇(du),u,∇(u)) ) ) * dΩ
+
+# Convection derivative: Picard iteration
+function p_dc_mhd(x,dx,dy,α,dΩ)
   u, p, j, φ = x
   du , dp , dj , dφ  = dx
   v_u, v_p, v_j, v_φ = dy
-  dc_mhd_u_u(u,du,v_u,α,dΩ)
+  p_dc_mhd_u_u(u,du,v_u,α,dΩ)
 end
-dc_mhd_u_u(u,du,v_u,α,dΩ) = ∫( α*v_u⋅( dconv∘(du,∇(du),u,∇(u)) ) ) * dΩ
+p_dc_mhd_u_u(u,du,v_u,α,dΩ) = ∫( α*v_u⋅( conv∘(u,∇(du)) ) ) * dΩ
 
 # Augmented lagrangian
 

test/dev/bugfix.jl

@@ -41,7 +41,7 @@ params[:fluid] = Dict{Symbol,Any}(
   :f => VectorValue(0.0,0.0,0.0),
   :B => VectorValue(0.0,1.0,0.0),
   :ζ => 1.0,
-  :convection => true,
+  :convection => :newton,
 )
 u_in = GridapMHD.u_inlet(:parabolic,1.0,4.0,1.0)
 params[:bcs] = Dict{Symbol,Any}(