From e9e7c9adc46ce687a11572d225b9b87286607304 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 28 Dec 2023 18:36:11 +1100 Subject: [PATCH] Started NonlinearSolvers --- src/BlockSolvers/BlockSolverInterfaces.jl | 70 +++++++++++---------- src/BlockSolvers/BlockSolvers.jl | 2 + src/GridapSolvers.jl | 2 + src/NonlinearSolvers/NewtonRaphsonSolver.jl | 70 +++++++++++++++++++++ src/NonlinearSolvers/NonlinearSolvers.jl | 19 ++++++ 5 files changed, 131 insertions(+), 32 deletions(-) create mode 100644 src/NonlinearSolvers/NewtonRaphsonSolver.jl create mode 100644 src/NonlinearSolvers/NonlinearSolvers.jl diff --git a/src/BlockSolvers/BlockSolverInterfaces.jl b/src/BlockSolvers/BlockSolverInterfaces.jl index 38bc28ca..0708b076 100644 --- a/src/BlockSolvers/BlockSolverInterfaces.jl +++ b/src/BlockSolvers/BlockSolverInterfaces.jl @@ -1,8 +1,34 @@ abstract type SolverBlock end + abstract type LinearSolverBlock <: SolverBlock end abstract type NonlinearSolverBlock <: SolverBlock end +is_nonlinear(::LinearSolverBlock) = false +is_nonlinear(::NonlinearSolverBlock) = true + +function instantiate_block_cache(block::LinearSolverBlock,mat::AbstractMatrix) + @abstractmethod +end +function instantiate_block_cache(block::NonlinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) + @abstractmethod +end +function instantiate_block_cache(block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) + instantiate_block_cache(block,mat) +end + +function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix) + return cache +end +function update_block_cache!(cache,block::NonlinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) + @abstractmethod +end +function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) + update_block_cache!(cache,block,mat) +end + +# MatrixBlock + struct MatrixBlock{A} <: LinearSolverBlock mat :: A function MatrixBlock(mat::AbstractMatrix) @@ -11,9 +37,18 @@ struct MatrixBlock{A} <: LinearSolverBlock end end +instantiate_block_cache(block::MatrixBlock,::AbstractMatrix) = block.mat + +# SystemBlocks + struct LinearSystemBlock <: LinearSolverBlock end struct NonlinearSystemBlock <: NonlinearSolverBlock end +instantiate_block_cache(block::LinearSystemBlock,mat::AbstractMatrix) = mat +instantiate_block_cache(block::NonlinearSystemBlock,mat::AbstractMatrix,::AbstractVector) = mat +update_block_cache!(cache,block::NonlinearSystemBlock,mat::AbstractMatrix,::AbstractVector) = mat + +# BiformBlock/TriformBlock struct BiformBlock <: LinearSolverBlock f :: Function trial :: FESpace @@ -40,50 +75,21 @@ struct TriformBlock <: NonlinearSolverBlock end end -# Instantiate blocks - -function instantiate_block_cache(block::LinearSolverBlock,mat::AbstractMatrix) - @abstractmethod -end -function instantiate_block_cache(block::NonlinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) - @abstractmethod -end -function instantiate_block_cache(block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) - instantiate_block_cache(block,mat) -end - -function instantiate_block_cache(block::MatrixBlock,mat::AbstractMatrix) - return block.mat -end function instantiate_block_cache(block::BiformBlock,mat::AbstractMatrix) return assemble_matrix(block.f,block.assem,block.trial,block.test) end -instantiate_block_cache(block::LinearSystemBlock,mat::AbstractMatrix) = mat function instantiate_block_cache(block::TriformBlock,mat::AbstractMatrix,x::AbstractVector) uh = FEFunction(block.trial,x) f(u,v) = block.f(uh,u,v) return assemble_matrix(f,block.assem,block.trial,block.test) end -instantiate_block_cache(block::NonlinearSystemBlock,mat::AbstractMatrix,x::AbstractVector) = mat - -# Update blocks - -function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix) - return cache -end -function update_block_cache!(cache,block::NonlinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) - @abstractmethod -end -function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) - update_block_cache!(cache,block,mat) -end function update_block_cache!(cache,block::TriformBlock,mat::AbstractMatrix,x::AbstractVector) uh = FEFunction(block.trial,x) f(u,v) = block.f(uh,u,v) assemble_matrix!(mat,f,block.assem,block.trial,block.test) end -function update_block_cache!(cache,block::NonlinearSystemBlock,mat::AbstractMatrix,x::AbstractVector) - return cache -end + +# CompositeBlock +# How do we deal with different sparsity patterns? Not trivial... diff --git a/src/BlockSolvers/BlockSolvers.jl b/src/BlockSolvers/BlockSolvers.jl index 2a7d0fa9..921f01b7 100644 --- a/src/BlockSolvers/BlockSolvers.jl +++ b/src/BlockSolvers/BlockSolvers.jl @@ -21,6 +21,8 @@ module BlockSolvers export BlockFEOperator + export MatrixBlock, LinearSystemBlock, NonlinearSystemBlock, BiformBlock, TriformBlock + export BlockDiagonalSolver export BlockTriangularSolver end diff --git a/src/GridapSolvers.jl b/src/GridapSolvers.jl index 5316d068..516819f2 100644 --- a/src/GridapSolvers.jl +++ b/src/GridapSolvers.jl @@ -5,12 +5,14 @@ module GridapSolvers include("BlockSolvers/BlockSolvers.jl") include("LinearSolvers/LinearSolvers.jl") include("PatchBasedSmoothers/PatchBasedSmoothers.jl") + include("NonlinearSolvers/NonlinearSolvers.jl") using GridapSolvers.SolverInterfaces using GridapSolvers.MultilevelTools using GridapSolvers.BlockSolvers using GridapSolvers.LinearSolvers using GridapSolvers.PatchBasedSmoothers + using GridapSolvers.NonlinearSolvers # MultilevelTools export get_parts, generate_level_parts, generate_subparts diff --git a/src/NonlinearSolvers/NewtonRaphsonSolver.jl b/src/NonlinearSolvers/NewtonRaphsonSolver.jl new file mode 100644 index 00000000..08e7bdac --- /dev/null +++ b/src/NonlinearSolvers/NewtonRaphsonSolver.jl @@ -0,0 +1,70 @@ + +# TODO: This should be called NewtonRaphsonSolver, but it would clash with Gridap. +struct NewtonSolver <: Algebra.NonlinearSolver + ls ::Algebra.LinearSolver + log::ConvergenceLog{Float64} +end + +function NewtonSolver(ls;maxiter=100,atol=1e-12,rtol=1.e-6,verbose=0,name="Newton-Raphson") + tols = SolverTolerances{Float64}(;maxiter=maxiter,atol=atol,rtol=rtol) + log = ConvergenceLog(name,tols;verbose=verbose) + return NewtonSolver(ls,log) +end + +struct NewtonCache + A::AbstractMatrix + b::AbstractVector + dx::AbstractVector + ns::NumericalSetup +end + +function Algebra.solve!(x::AbstractVector,nls::NewtonSolver,op::NonlinearOperator,cache::Nothing) + b = residual(op, x) + A = jacobian(op, x) + dx = similar(b) + ss = symbolic_setup(nls.ls, A) + ns = numerical_setup(ss,A) + _solve_nr!(x,A,b,dx,ns,nls,op) + return NewtonCache(A,b,dx,ns) +end + +function Algebra.solve!(x::AbstractVector,nls::NewtonSolver,op::NonlinearOperator,cache::NewtonCache) + A,b,dx,ns = cache.A, cache.b, cache.dx, cache.ns + residual!(b, op, x) + jacobian!(A, op, x) + numerical_setup!(ns,A) + _solve_nr!(x,A,b,dx,ns,nls,op) + return cache +end + +function _solve_nr!(x,A,b,dx,ns,nls,op) + log = nls.log + + # Check for convergence on the initial residual + res = norm(b) + done = init!(log,res) + + # Newton-like iterations + while !done + + # Solve linearized problem + rmul!(b,-1) + solve!(dx,ns,b) + x .+= dx + + # Check convergence for the current residual + residual!(b, op, x) + res = norm(b) + done = update!(log,res) + + if !done + # Update jacobian and solver + jacobian!(A, op, x) + numerical_setup!(ns,A) + end + + end + + finalize!(log,res) + return x +end diff --git a/src/NonlinearSolvers/NonlinearSolvers.jl b/src/NonlinearSolvers/NonlinearSolvers.jl new file mode 100644 index 00000000..0ee37902 --- /dev/null +++ b/src/NonlinearSolvers/NonlinearSolvers.jl @@ -0,0 +1,19 @@ +module NonlinearSolvers + using LinearAlgebra + using SparseArrays + using SparseMatricesCSR + using BlockArrays + using IterativeSolvers + + using Gridap + using Gridap.Helpers, Gridap.Algebra, Gridap.CellData, Gridap.Arrays, Gridap.FESpaces, Gridap.MultiField + using PartitionedArrays + using GridapDistributed + + using GridapSolvers.MultilevelTools + using GridapSolvers.SolverInterfaces + + include("NewtonRaphsonSolver.jl") + export NewtonSolver + +end \ No newline at end of file