adding initial stab at integrator interface.

ext/FiniteElementContainersAdaptExt.jl

@@ -152,7 +152,16 @@
   fspace = adapt(to, var.fspace)
   return VectorFunction{syms, typeof(fspace)}(fspace)
 end
-
+
+# Integrators
+function Adapt.adapt_structure(to, integrator::QuasiStaticIntegrator)
+  return QuasiStaticIntegrator(
+    adapt(to, integrator.solution),
+    adapt(to, integrator.solver)
+  )
+end
+
+# parameters
 function Adapt.adapt_structure(to, p::FiniteElementContainers.Parameters)
   return FiniteElementContainers.Parameters(
     adapt(to, p.dirichlet_bcs),

src/FiniteElementContainers.jl

@@ -61,6 +61,11 @@ export SymmetricTensorFunction
 export TensorFunction
 export VectorFunction
 
+# Integrators
+export AbstractIntegrator
+export QuasiStaticIntegrator
+export evolve!
+
 # Meshes
 export FileMesh
 export UnstructuredMesh
@@ -84,6 +89,7 @@ export sideset_names
 
 # Parameters
 export Parameters
+export TimeStepper
 export create_parameters
 
 # Physics
@@ -100,7 +106,7 @@ export write_times
 # export AbstractPreconditioner
 # export AbstractSolver
 export DirectLinearSolver
-export IterativeSolver
+export IterativeLinearSolver
 export NewtonSolver
 export solve!
 
@@ -139,5 +145,6 @@ include("PostProcessors.jl")
 #
 include("Parameters.jl")
 include("solvers/Solvers.jl")
+include("integrators/Integrators.jl")
 
 end # module

src/Parameters.jl

@@ -63,10 +63,10 @@
 # 4. add different fspace types
 # 5. convert vectors of dbcs/nbcs into namedtuples
 function Parameters(
-  assembler, physics, 
-  dbcs=nothing, 
-  nbcs=nothing, 
-  times=nothing
+  assembler, physics,
+  dirichlet_bcs, 
+  neumann_bcs, 
+  times
 )
   h1_dbcs = create_bcs(assembler, H1Field)
   h1_field = create_field(assembler, H1Field)
@@ -76,22 +76,22 @@
   state_old = nothing
   state_new = nothing
 
-  if dbcs !== nothing
-    syms = map(x -> Symbol("dirichlet_bc_$x"), 1:length(dbcs))
+  if dirichlet_bcs !== nothing
+    syms = map(x -> Symbol("dirichlet_bc_$x"), 1:length(dirichlet_bcs))
     # dbcs = NamedTuple{tuple(syms...)}(tuple(dbcs...))
     # dbcs = DirichletBCContainer(dbcs, size(assembler.dof.H1_vars[1].fspace.coords, 1))
-    dbcs = DirichletBCContainer.((assembler.dof,), dbcs)
-    temp_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, dbcs)
+    dirichlet_bcs = DirichletBCContainer.((assembler.dof,), dirichlet_bcs)
+    temp_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, dirichlet_bcs)
     temp_dofs = unique(sort(temp_dofs))
-    dbcs = NamedTuple{tuple(syms...)}(tuple(dbcs...))
+    dirichlet_bcs = NamedTuple{tuple(syms...)}(tuple(dirichlet_bcs...))
     # TODO eventually do something different here
     # update_dofs!(assembler.dof, dbcs.bookkeeping.dofs)
     update_dofs!(assembler.dof, temp_dofs)
   end
 
-  if nbcs !== nothing
-    syms = map(x -> Symbol("neumann_bc_$x"), 1:length(nbcs))
-    dbcs = NamedTuple{tuple(syms...)}(tuple(nbcs...))
+  if neumann_bcs !== nothing
+    syms = map(x -> Symbol("neumann_bc_$x"), 1:length(neumann_bcs))
+    neumann_bcs = NamedTuple{tuple(syms...)}(tuple(neumann_bcs...))
   end
 
   # dummy time stepper for a static problem
@@ -100,7 +100,8 @@
   end
 
   p = Parameters(
-    dbcs, nbcs, 
+    dirichlet_bcs,
+    neumann_bcs, 
     times,
     physics, 
     properties, 
@@ -113,13 +114,23 @@
   return p
 end
 
-function create_parameters(assembler, physics, dbcs=nothing, nbcs=nothing)
-  return Parameters(assembler, physics, dbcs, nbcs)
+function create_parameters(
+  assembler, physics; 
+  dirichlet_bcs=nothing, 
+  neumann_bcs=nothing,
+  times=nothing
+)
+  return Parameters(assembler, physics, dirichlet_bcs, neumann_bcs, times)
 end
 
 function update_dofs!(asm::SparseMatrixAssembler, p::Parameters)
   # temp_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, values(p.dirichlet_bcs))
   update_dofs!(asm, p.dirichlet_bcs)
   Base.resize!(p.h1_dbcs, length(asm.dof.H1_bc_dofs))
   return nothing
-end
+end
+
+function update_time!(p::Parameters)
+  fill!(p.times.time_current, current_time(p.times) + time_step(p.times))
+  return nothing
+end

src/integrators/Integrators.jl

@@ -1,6 +1,14 @@
-abstract type AbstractIntegrator{Solver} end
+abstract type AbstractIntegrator{
+  Solution <: AbstractArray{<:Number, 1},
+  Solver <: AbstractNonLinearSolver
+} end
 # abstract type AbstractFirstOrderIntegrator <: AbstractIntegrator end
 # abstract type AbstractSecondOrderIntegrator <: AbstractIntegrator end
-abstract type AbstractStaticIntegrator{Solver, Solution} <: AbstractIntegrator{Solver, Solution} end
+abstract type AbstractStaticIntegrator{
+  Solution,
+  Solver
+} <: AbstractIntegrator{Solution, Solver} end
 
-include("StaticIntegrator.jl")
+function evolve! end
+
+include("QuasiStaticIntegrator.jl")

src/integrators/QuasiStaticIntegrator.jl

@@ -0,0 +1,21 @@
+struct QuasiStaticIntegrator{
+  Solution, 
+  Solver
+} <: AbstractStaticIntegrator{Solution, Solver}
+  solver::Solver
+  solution::Solution
+end
+
+function QuasiStaticIntegrator(solver::AbstractNonLinearSolver)
+  solution = similar(solver.linear_solver.ΔUu)
+  fill!(solution, zero(eltype(solution)))
+  return QuasiStaticIntegrator(solver, solution)
+end
+
+function evolve!(integrator::QuasiStaticIntegrator, p)
+  update_time!(p)
+  update_bcs!(H1Field, integrator.solver, integrator.solution, p)
+  solve!(integrator.solver, integrator.solution, p)
+  # copyto!(integrator.solution, Uu)
+  return nothing
+end

src/solvers/DirectLinearSolver.jl

@@ -0,0 +1,31 @@
+struct DirectLinearSolver{
+  A <: SparseMatrixAssembler, 
+  P, 
+  Inc <: AbstractArray{<:Number, 1}
+} <: AbstractLinearSolver{A, P, Inc}
+  assembler::A
+  preconditioner::P
+  # TODO add some tolerances
+  # what's the best way to do this with general solvers?
+  ΔUu::Inc
+end
+
+function DirectLinearSolver(assembler::SparseMatrixAssembler)
+  preconditioner = I
+  ΔUu = similar(assembler.residual_unknowns)
+  fill!(ΔUu, zero(eltype(ΔUu)))
+  return DirectLinearSolver(assembler, preconditioner, ΔUu)
+end 
+
+function solve!(solver::DirectLinearSolver, Uu, p)
+  update_field_unknowns!(p.h1_field, solver.assembler.dof, Uu)
+  assemble!(solver.assembler, p.physics, p.h1_field, :residual_and_stiffness)
+  R = residual(solver.assembler)
+  K = stiffness(solver.assembler)
+  # TODO specialize to backend solvers if they exists
+  # solver.ΔUu .= -K \ R
+  copyto!(solver.ΔUu, -K \ R)
+  # update_field_unknowns!(p.h1_field, solver.assembler.dof, solver.ΔUu, +)
+  map!((x, y) -> x + y, Uu, Uu, solver.ΔUu)
+  return nothing
+end

src/solvers/IterativeLinearSolver.jl

@@ -0,0 +1,38 @@
+struct IterativeLinearSolver{
+  A <: AbstractAssembler,
+  P,
+  S,
+  U <: AbstractArray{<:Number, 1}
+} <: AbstractLinearSolver{A, P, U}
+  assembler::A
+  preconditioner::P
+  solver::S
+  ΔUu::U
+end
+
+function IterativeLinearSolver(assembler::SparseMatrixAssembler, solver_sym)
+  # TODO
+  preconditioner = I
+  # ΔUu = KA.zeros(KA.get_backend)
+  ΔUu = similar(assembler.residual_unknowns)
+  fill!(ΔUu, zero(eltype(ΔUu)))
+  n = length(ΔUu)
+  solver = eval(solver_sym)(n, n, typeof(ΔUu))
+  return IterativeLinearSolver(assembler, preconditioner, solver, ΔUu)
+end
+
+# TODO specialize for operator like assemblers
+function solve!(solver::IterativeLinearSolver, Uu, p)
+  # update unknown dofs
+  update_field_unknowns!(p.h1_field, solver.assembler.dof, Uu)
+  # assemble relevant fields
+  assemble!(solver.assembler, p.physics, p.h1_field, :residual)
+  assemble!(solver.assembler, p.physics, p.h1_field, :stiffness)
+  # solve and fetch solution
+  Krylov.solve!(solver.solver, stiffness(solver.assembler), residual(solver.assembler))
+  ΔUu = -Krylov.solution(solver.solver)
+  # make necessary copies and updates
+  copyto!(solver.ΔUu, ΔUu)
+  map!((x, y) -> x + y, Uu, Uu, ΔUu)
+  return nothing
+end

src/solvers/NewtonSolver.jl

@@ -0,0 +1,21 @@
+struct NewtonSolver{L, T} <: AbstractNonLinearSolver{L}
+  linear_solver::L
+  #
+  max_iters::Int
+  abs_increment_tol::T
+  abs_residual_tol::T
+end
+
+function NewtonSolver(linear_solver)
+  return NewtonSolver(linear_solver, 10, 1e-12, 1e-12)
+end
+
+function solve!(solver::NewtonSolver, Uu, p)
+  for n in 1:solver.max_iters
+    solve!(solver.linear_solver, Uu, p)
+    @show norm(solver.linear_solver.ΔUu) norm(residual(solver.linear_solver.assembler))
+    if norm(solver.linear_solver.ΔUu) < 1e-12 || norm(residual(solver.linear_solver.assembler)) < 1e-12
+      break
+    end
+  end
+end