From 54d8e04eaf45fb5d2ed08ce2f64de1027d8acf7f Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 7 Nov 2024 00:08:35 +1100 Subject: [PATCH] Added Picard iteration --- src/Applications/cavity.jl | 3 ++- src/Applications/channel.jl | 3 ++- src/Applications/expansion.jl | 3 ++- src/Main.jl | 8 +++++--- src/Solvers/gmg.jl | 4 ++-- src/parameters.jl | 33 ++++++++++++++++++++------------ src/weakforms.jl | 36 ++++++++++++++++++++++++++++------- test/dev/bugfix.jl | 2 +- 8 files changed, 64 insertions(+), 28 deletions(-) diff --git a/src/Applications/cavity.jl b/src/Applications/cavity.jl index a38547a..1944948 100644 --- a/src/Applications/cavity.jl +++ b/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}( diff --git a/src/Applications/channel.jl b/src/Applications/channel.jl index 34ad77c..e07c4a0 100644 --- a/src/Applications/channel.jl +++ b/src/Applications/channel.jl @@ -66,7 +66,7 @@ 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 @@ -74,6 +74,7 @@ function _channel(; @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}( diff --git a/src/Applications/expansion.jl b/src/Applications/expansion.jl index df1175b..041fe07 100644 --- a/src/Applications/expansion.jl +++ b/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}( diff --git a/src/Main.jl b/src/Main.jl index ec71d0a..39f40a0 100644 --- a/src/Main.jl +++ b/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 diff --git a/src/Solvers/gmg.jl b/src/Solvers/gmg.jl index 76b75d1..008bf66 100644 --- a/src/Solvers/gmg.jl +++ b/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) diff --git a/src/parameters.jl b/src/parameters.jl index 6eff6f2..76e3f5f 100644 --- a/src/parameters.jl +++ b/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 diff --git a/src/weakforms.jl b/src/weakforms.jl index 73a6308..97ce594 100644 --- a/src/weakforms.jl +++ b/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,6 +293,8 @@ 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 @@ -292,13 +302,25 @@ function c_mhd(x,dy,α,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 diff --git a/test/dev/bugfix.jl b/test/dev/bugfix.jl index 2059627..08b39d1 100644 --- a/test/dev/bugfix.jl +++ b/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}(