Skip to content

Commit

Permalink
First steps towards continuation
Browse files Browse the repository at this point in the history
  • Loading branch information
JordiManyer committed Oct 31, 2024
1 parent 3a4780a commit 8d3a181
Show file tree
Hide file tree
Showing 5 changed files with 104 additions and 93 deletions.
3 changes: 1 addition & 2 deletions src/Applications/channel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -173,10 +173,9 @@ function _channel(;

toc!(t,"pre_process")


params[:solver][:niter] = niter
if initial_value == :inlet
params[:solver][:initial_values] = Dict(
params[:x0] = Dict(
:u=>ū,:j=>j_zero,:p=>0.0,=>0.0)
end
if params[:solver][:solver] == :petsc
Expand Down
4 changes: 2 additions & 2 deletions src/Applications/expansion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ function _expansion(;
params[:fespaces] = Dict{Symbol,Any}(
:order_u => order,
:order_j => order_j,
:rt_scaling => rt_scaling ? 1.0/_get_mesh_size(model) : nothing
:rt_scaling => rt_scaling ? 1.0/get_mesh_size(model) : nothing
)

params[:fluid] = Dict{Symbol,Any}(
Expand Down Expand Up @@ -176,7 +176,7 @@ function _expansion(;

j_zero = VectorValue(0.0,0.0,0.0)
if initial_value == :inlet
params[:solver][:initial_values] = Dict(
params[:x0] = Dict(
:u=>u_in,:j=>j_zero,:p=>0.0,=>0.0
)
end
Expand Down
28 changes: 11 additions & 17 deletions src/Applications/transient.jl
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,6 @@ function _transient(;
solver = :julia,
verbose = true,
man_solution = nothing,
initial_value_type = :zero,
max_error = 0.0,
time_solver = :theta,
θ = 0.5,
Expand Down Expand Up @@ -150,30 +149,25 @@ function _transient(;
end

# Setup ODE solver
ode_solver_params = Dict(=>θ)

params[:transient] = Dict{Symbol,Any}(
:solver => ode_solver,
:t0 => t0,
:tf => tf,
:Δt => Δt,
)
params[:transient][:solver] = Dict{Symbol,Any}(
:solver => time_solver,
=> θ
)
if is_manufactured
initial_values = Dict(
params[:x0] = Dict{Symbol,Any}(
:u => u(t0),
:p => p(t0),
:j => j(t0),
=> φ(t0),
)
initial_value_type = :value
else
initial_values = nothing
end

params[:ode] = Dict(
:solver => time_solver,
:t0 => t0,
:tf => tf,
:Δt => Δt,
:solver_params => ode_solver_params,
:U0 => initial_value_type,
:initial_values => initial_values,
)

toc!(t,"pre_process")

# Lazy ODE solver
Expand Down
68 changes: 28 additions & 40 deletions src/Main.jl
Original file line number Diff line number Diff line change
Expand Up @@ -153,12 +153,12 @@ function main(_params::Dict;output::Dict=Dict{Symbol,Any}())
toc!(t,"solve")
else
op = _fe_operator(U,V,params)
xh = _allocate_solution(op,params)
xh = initial_guess(op,params)
if params[:solve]
solver = _solver(op,params)
toc!(t,"solver_setup")
tic!(t;barrier=true)
xh,cache = _solve(xh,solver,op,params)
xh, cache = _solve(xh,solver,op,params)
solver_postpro = params[:solver][:solver_postpro]
solver_postpro(cache,output)
toc!(t,"solve")
Expand All @@ -183,13 +183,12 @@ end

function _solver(op,params)
solver = _solver(Val(params[:solver][:solver]),op,params)
if params[:transient]
if !isnothing(params[:transient])
solver = _ode_solver(solver,params)
end
return solver
end

#_solver(::Val{:julia},op,params) = NLSolver(show_trace=true,method=:newton)
function _solver(::Val{:julia},op,params)
verbose = i_am_main(get_parts(params[:model]))
GridapSolvers.NewtonSolver(
Expand All @@ -200,16 +199,16 @@ _solver(::Val{:petsc},op,params) = PETScNonlinearSolver()
_solver(::Val{:li2019},op,params) = Li2019Solver(op,params)
_solver(::Val{:badia2024},op,params) = Badia2024Solver(op,params)

_ode_solver(solver,params) = _ode_solver(Val(params[:ode][:solver]),solver,params)
_ode_solver(solver,params) = _ode_solver(Val(params[:transient][:solver][:solver]),solver,params)

function _ode_solver(::Val{:theta},solver,params)
Δt = params[:ode][:Δt]
θ = params[:ode][:solver_params][]
Δt = params[:transient][:Δt]
θ = params[:transient][:solver][]
ThetaMethod(solver,Δt,θ)
end

function _ode_solver(::Val{:forward},solver,params)
Δt = params[:ode][:Δt]
Δt = params[:transient][:Δt]
ForwardEuler(solver,Δt)
end

Expand Down Expand Up @@ -305,18 +304,18 @@ end
function _trial_fe_space(V,v_bc,transient)
if isnothing(v_bc) || (isa(v_bc,Number) && iszero(v_bc))
return V
elseif !transient
return TrialFESpace(V,v_bc)
else
elseif !isnothing(transient)
return TransientTrialFESpace(V,v_bc)
else
return TrialFESpace(V,v_bc)
end
end

# FEOperator

function _fe_operator(U,V,params)
mfs = _multi_field_style(params)
if !params[:transient]
if isnothing(params[:transient])
_fe_operator(mfs,U,V,params)
else
_ode_fe_operator(mfs,U,V,params)
Expand Down Expand Up @@ -410,43 +409,32 @@ end

# Solve

_solve(xh,solver,op,params) = _solve(Val(params[:transient]),xh,solver,op,params)

_solve(::Val{false},xh,solver,op,params) = solve!(xh,solver,op)
_solve(xh,solver,op::FEOperator,params) = solve!(xh,solver,op)

function _solve(::Val{true},xh,solver,op,params)
xh0 = initial_value(op,params)
t0,tf = time_interval(params)
function _solve(xh,solver,op::TransientFEOperator,params)
t0, tf = params[:transient][:t0], params[:transient][:tf]
cache = nothing
solve(solver,op,t0,tf,xh0), cache
solve(solver,op,t0,tf,xh), cache
end

initial_value(op,params) = initial_value(Val(params[:ode][:U0]),op,params)
initial_value(::Val{:zero},op,params) = zero(get_trial(op))
initial_value(::Val{:solve},op,params) = @notimplemented
function initial_value(::Val{:value},op,params)
t0 = params[:ode][:t0]
U0 = get_trial(op)(t0)
v = params[:ode][:initial_values]
interpolate([v[:u],v[:p],v[:j],v[]],U0)
end

time_interval(params) = ( params[:ode][:t0], params[:ode][:tf] )

_allocate_solution(op,params) = _allocate_solution(op,params,params[:solver][:initial_values])

function _allocate_solution(op::FEOperator,params,x0::Dict)
x0 = [x0[:u],x0[:p],x0[:j],x0[]]
function initial_guess(op::FEOperator,params)
U = get_trial(op)
interpolate(x0,U)
initial_guess(params[:x0],U,op,params)
end

function _allocate_solution(op::FEOperator,params,::Nothing)
zero(get_trial(op))
function initial_guess(op::TransientFEOperator,params)
t0 = params[:transient][:t0]
U0 = get_trial(op)(t0)
initial_guess(params[:x0],U0,op,params)
end

function _allocate_solution(op::TransientFEOperator,args...)
nothing
initial_guess(x0::Symbol,trial,op,params) = initial_guess(Val(x0),trial,op,params)
initial_guess(::Val{:zero},trial,op,params) = zero(trial)
initial_guess(::Val{:solve},trial,op,params) = @notimplemented

function initial_guess(x0::Dict,trial,op,params)
vals = [x0[:u],x0[:p],x0[:j],x0[]]
interpolate(vals,trial)
end

# Mesh sizes
Expand Down
94 changes: 62 additions & 32 deletions src/parameters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -84,54 +84,55 @@ It also checks the validity of the main parameter dictionary `params`.
"""
function add_default_params(_params)
mandatory = Dict(
:ptimer=>false,
:debug=>false,
:solve=>true,
:res_assemble=>false,
:jac_assemble=>false,
:model=>true,
:fluid=>true,
:solid=>false,
:bcs=>true,
:fespaces=>false,
:solver=>true,
:multigrid=>false,
:check_valid=>false,
:ode=>false,
:transient=>false,
:ptimer => false,
:debug => false,
:solve => true,
:res_assemble => false,
:jac_assemble => false,
:model => true,
:fluid => true,
:solid => false,
:bcs => true,
:fespaces => false,
:solver => true,
:multigrid => false,
:check_valid => false,
:transient => false,
:x0 => false,
)
_check_mandatory(_params,mandatory,"")
optional = Dict(
:ptimer=>default_ptimer(_params[:model]),
:debug=>false,
:res_assemble=>false,
:jac_assemble=>false,
:solid=>nothing,
:fespaces=>nothing,
:multigrid=>nothing,
:check_valid=>true,
:ode=>nothing,
:transient=>default_transient(_params),
:ptimer => default_ptimer(_params[:model]),
:debug => false,
:res_assemble => false,
:jac_assemble => false,
:solid => nothing,
:fespaces => nothing,
:multigrid => nothing,
:check_valid => true,
:transient => nothing,
:x0 => :zero,
)
params = _add_optional(_params,mandatory,optional,_params,"")
_check_unused(params,mandatory,params,"")
# Process sub-params
params[:fluid] = params_fluid(params)
if params[:solid] !== nothing
if !isnothing(params[:solid])
params[:solid] = params_solid(params)
end
params[:bcs] = params_bcs(params)
params[:fespaces] = params_fespaces(params)
params[:solver] = params_solver(params)
params[:multigrid] = params_multigrid(params)
if !isnothing(params[:transient])
params[:transient] = params_transient(params)
end
params
end

default_ptimer(model) = PTimer(DebugArray(LinearIndices((1,))))
default_ptimer(model::GridapDistributed.DistributedDiscreteModel) = PTimer(get_parts(model))

default_transient(params) = haskey(params,:ode) && !isnothing(params[:ode])

"""
Valid keys for `params[:solver]` are the following:
Expand Down Expand Up @@ -172,7 +173,6 @@ function default_solver_params(::Val{:julia})
:solver_postpro => ((cache,info) -> nothing),
:niter => 10,
:rtol => 1e-5,
:initial_values => nothing,
)
end

Expand All @@ -185,7 +185,6 @@ function default_solver_params(::Val{:petsc})
:petsc_options => "-snes_monitor -ksp_error_if_not_converged true -ksp_converged_reason -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_7 0",
:niter => 100,
:rtol => 1e-5,
:initial_values => nothing,
)
end

Expand All @@ -199,7 +198,6 @@ function default_solver_params(::Val{:li2019})
:block_solvers => [:petsc_mumps,:petsc_gmres_schwarz,:petsc_cg_jacobi,:petsc_cg_jacobi],
:niter => 80,
:rtol => 1e-5,
:initial_values => nothing,
)
end

Expand All @@ -213,7 +211,6 @@ function default_solver_params(::Val{:badia2024})
:block_solvers => [:petsc_mumps,:petsc_cg_jacobi,:petsc_cg_jacobi],
:niter => 80,
:rtol => 1e-5,
:initial_values => nothing,
)
end

Expand Down Expand Up @@ -606,3 +603,36 @@ function params_bcs_stabilization(params::Dict{Symbol,Any})
optional = Dict(:domain=>params[:fluid][:domain])
_check_mandatory_and_add_optional_weak(params[:bcs][:stabilization],mandatory,optional,params,"[:bcs][:stabilization]")
end

function params_transient(params::Dict{Symbol,Any})
mandatory = Dict(
:t0 => true,
:tf => true,
:Δt => true,
:solver => false,
)
optional = Dict(
:solver => :theta,
)
transient = _check_mandatory_and_add_optional_weak(params[:transient],mandatory,optional,params,"[:transient]")
if isa(transient[:solver],Symbol)
transient[:solver] = default_transient_solver_params(transient[:solver])
else
optional_solver = default_solver_params(transient[:solver][:solver])
transient[:solver] = _add_optional_weak(transient[:solver],optional_solver,params,"[:solver]")
end
return transient
end

function default_solver_params(::Val{:theta})
return Dict(
:solver => :theta,
=> 0.5,
)
end

function default_solver_params(::Val{:forward})
return Dict(
:solver => :forward,
)
end

0 comments on commit 8d3a181

Please sign in to comment.