Started NonlinearSolvers

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...

src/BlockSolvers/BlockSolvers.jl

@@ -21,6 +21,8 @@ module BlockSolvers
 
   export BlockFEOperator
 
+  export MatrixBlock, LinearSystemBlock, NonlinearSystemBlock, BiformBlock, TriformBlock
+
   export BlockDiagonalSolver
   export BlockTriangularSolver
 end

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

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

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