diff --git a/ext/FiniteElementContainersAdaptExt.jl b/ext/FiniteElementContainersAdaptExt.jl
index 256e5bc..aab1876 100644
--- a/ext/FiniteElementContainersAdaptExt.jl
+++ b/ext/FiniteElementContainersAdaptExt.jl
@@ -152,7 +152,16 @@ function Adapt.adapt_structure(to, var::VectorFunction)
   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),
diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl
index 4fbb3e0..3712a2c 100644
--- a/src/FiniteElementContainers.jl
+++ b/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
diff --git a/src/Parameters.jl b/src/Parameters.jl
index 319099d..2ec021e 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -63,10 +63,10 @@ end
 # 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 @@ function Parameters(
   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 @@ function Parameters(
   end
 
   p = Parameters(
-    dbcs, nbcs, 
+    dirichlet_bcs,
+    neumann_bcs, 
     times,
     physics, 
     properties, 
@@ -113,8 +114,13 @@ function Parameters(
   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)
@@ -122,4 +128,9 @@ function update_dofs!(asm::SparseMatrixAssembler, p::Parameters)
   update_dofs!(asm, p.dirichlet_bcs)
   Base.resize!(p.h1_dbcs, length(asm.dof.H1_bc_dofs))
   return nothing
-end
\ No newline at end of file
+end
+
+function update_time!(p::Parameters)
+  fill!(p.times.time_current, current_time(p.times) + time_step(p.times))
+  return nothing
+end
diff --git a/src/integrators/Integrators.jl b/src/integrators/Integrators.jl
index 73de0f5..f1bc01d 100644
--- a/src/integrators/Integrators.jl
+++ b/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")
diff --git a/src/integrators/QuasiStaticIntegrator.jl b/src/integrators/QuasiStaticIntegrator.jl
new file mode 100644
index 0000000..1d7b0b6
--- /dev/null
+++ b/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
diff --git a/src/integrators/StaticIntegrator.jl b/src/integrators/StaticIntegrator.jl
deleted file mode 100644
index 6aa374e..0000000
--- a/src/integrators/StaticIntegrator.jl
+++ /dev/null
@@ -1,9 +0,0 @@
-# dummy placeholder for static problems to have a consistent interface
-
-struct StaticIntegrator <: AbstractStaticIntegrator{Sol}
-  U::Sol
-end
-
-function StaticIntegrator()
-
-end
diff --git a/src/solvers/DirectLinearSolver.jl b/src/solvers/DirectLinearSolver.jl
new file mode 100644
index 0000000..529adba
--- /dev/null
+++ b/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
diff --git a/src/solvers/IterativeLinearSolver.jl b/src/solvers/IterativeLinearSolver.jl
new file mode 100644
index 0000000..c794776
--- /dev/null
+++ b/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
diff --git a/src/solvers/NewtonSolver.jl b/src/solvers/NewtonSolver.jl
new file mode 100644
index 0000000..f78db2f
--- /dev/null
+++ b/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
diff --git a/src/solvers/Solvers.jl b/src/solvers/Solvers.jl
index c789840..d94ee15 100644
--- a/src/solvers/Solvers.jl
+++ b/src/solvers/Solvers.jl
@@ -1,12 +1,19 @@
 # abstract type AbstractPreconditioner end
 abstract type AbstractSolver end
 
-abstract type AbstractLinearSolver{A <: AbstractAssembler, P, Inc} <: AbstractSolver end
-abstract type AbstractNonLinearSolver{L <: AbstractLinearSolver} <: AbstractSolver end
+abstract type AbstractLinearSolver{
+  A <: AbstractAssembler, 
+  P, 
+  U <: AbstractArray{<:Number, 1}
+} <: AbstractSolver end
+abstract type AbstractNonLinearSolver{
+  L <: AbstractLinearSolver
+} <: AbstractSolver end
 
 # traditional solver interface
 function preconditioner end
 function residual end
+function solve! end
 function tangent end
 function update_preconditioner! end
 
@@ -24,7 +31,7 @@ KA.@kernel function _update_bcs_kernel!(bc, p)
   p.h1_field[dof] = val
 end
 
-function _update_bcs!(type::Type{H1Field}, bc, p, backend::KA.Backend)
+function _update_bcs!(::Type{H1Field}, bc, p, backend::KA.Backend)
   kernel! = _update_bcs_kernel!(backend)
   kernel!(bc, p, ndrange=length(bc.vals))
   return nothing
@@ -51,98 +58,6 @@ function update_bcs!(::Type{H1Field}, solver::AbstractNonLinearSolver, Uu, p)
   return nothing
 end
 
-# optimization based solver interface
-function gradient end
-function hessian end
-function value end
-
-function solve! end
-
-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)
-  assemble!(solver.assembler, p.physics, p.h1_field, :residual_and_stiffness)
-  R = residual(solver.assembler)
-  K = stiffness(solver.assembler)
-  # TODO make a KA kernel for a copy here
-  # 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, +)
-  return nothing
-end
-
-struct IterativeSolver{
-  A,
-  P,
-  Inc,
-  S
-} <: AbstractLinearSolver{A, P, Inc}
-  assembler::A
-  preconditioner::P
-  ΔUu::Inc
-  solver::S
-end
-
-function IterativeSolver(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 IterativeSolver(assembler, preconditioner, ΔUu, solver)
-end
-
-# TODO specialize for operator like assemblers
-function solve!(solver::IterativeSolver, Uu, p)
-  # assemble!(solver.assembler, p.physics, p.h1_field, :residual_and_stiffness)
-  assemble!(solver.assembler, p.physics, p.h1_field, :residual)
-  assemble!(solver.assembler, p.physics, p.h1_field, :stiffness)
-  Krylov.solve!(solver.solver, stiffness(solver.assembler), residual(solver.assembler))
-  ΔUu = -Krylov.solution(solver.solver)
-  # solver.ΔUu .= ΔUu 
-  copyto!(solver.ΔUu, ΔUu)
-  update_field_unknowns!(p.h1_field, solver.assembler.dof, solver.ΔUu, +)
-  return nothing
-end
-
-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
+include("DirectLinearSolver.jl")
+include("IterativeLinearSolver.jl")
+include("NewtonSolver.jl")
diff --git a/test/mechanics/TestMechanics.jl b/test/mechanics/TestMechanics.jl
index f4c52e5..a6015bc 100644
--- a/test/mechanics/TestMechanics.jl
+++ b/test/mechanics/TestMechanics.jl
@@ -8,7 +8,7 @@ mesh_file = "./mechanics/mechanics.g"
 output_file = "./mechanics/mechanics.e"
 
 fixed(_, _) = 0.
-displace(_, _) = 1.e-3
+displace(_, t) = 1.e-3 * t
 
 struct Mechanics{Form} <: AbstractPhysics{2, 0, 0}
   formulation::Form
@@ -67,13 +67,16 @@ function mechanics_test()
   ]
 
   # pre-setup some scratch arrays
-  p = create_parameters(asm, physics, dbcs)
-  Uu = create_unknowns(asm)
+  times = TimeStepper(0., 1., 1)
+  p = create_parameters(asm, physics; dirichlet_bcs=dbcs, times=times)
+  # Uu = create_unknowns(asm)
 
-  solver = NewtonSolver(IterativeSolver(asm, :CgSolver))
-  update_bcs!(H1Field, solver, Uu, p)
+  solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
+  integrator = QuasiStaticIntegrator(solver)
+  evolve!(integrator, p)
+  # update_bcs!(H1Field, solver, Uu, p)
 
-  FiniteElementContainers.solve!(solver, Uu, p)
+  # FiniteElementContainers.solve!(solver, Uu, p)
 
   pp = PostProcessor(mesh, output_file, u)
   write_times(pp, 1, 0.0)
diff --git a/test/mechanics/TestMechanicsCUDA.jl b/test/mechanics/TestMechanicsCUDA.jl
index 0c1ed9b..29d8575 100644
--- a/test/mechanics/TestMechanicsCUDA.jl
+++ b/test/mechanics/TestMechanicsCUDA.jl
@@ -10,7 +10,7 @@ mesh_file = "./test/mechanics/mechanics.g"
 output_file = "./test/mechanics/mechanics.e"
 
 fixed(_, _) = 0.
-displace(_, _) = 1.e-3
+displace(_, t) = 1.e-3 * t
 
 struct Mechanics{Form} <: AbstractPhysics{2, 0, 0}
   formulation::Form
@@ -69,7 +69,8 @@ end
   ]
 
   # pre-setup some scratch arrays
-  p = create_parameters(asm, physics, dbcs)
+  times = TimeStepper(0., 1., 1)
+  p = create_parameters(asm, physics; dirichlet_bcs=dbcs, times=times)
 
   U = create_field(asm, H1Field)
   @time assemble!(asm, physics, U, :stiffness)
@@ -79,13 +80,15 @@ end
   p_gpu = p |> gpu
   asm_gpu = asm |> gpu
 
-  solver = NewtonSolver(IterativeSolver(asm_gpu, :CgSolver))
-  Uu = create_unknowns(asm_gpu)
+  solver = NewtonSolver(IterativeLinearSolver(asm_gpu, :CgSolver))
+  integrator = QuasiStaticIntegrator(solver)
+  evolve!(integrator, p_gpu)
+  # Uu = create_unknowns(asm_gpu)
 
-  update_bcs!(H1Field, solver, Uu, p_gpu)
+  # update_bcs!(H1Field, solver, Uu, p_gpu)
 
-  @time FiniteElementContainers.solve!(solver, Uu, p_gpu)
-  @time FiniteElementContainers.solve!(solver, Uu, p_gpu)
+  # @time FiniteElementContainers.solve!(solver, Uu, p_gpu)
+  # @time FiniteElementContainers.solve!(solver, Uu, p_gpu)
 
   p = p_gpu |> cpu
 
diff --git a/test/poisson/TestPoisson.jl b/test/poisson/TestPoisson.jl
index 5cf66e6..c485042 100644
--- a/test/poisson/TestPoisson.jl
+++ b/test/poisson/TestPoisson.jl
@@ -43,19 +43,14 @@ function poisson()
     DirichletBC(:u, :sset_3, bc_func),
     DirichletBC(:u, :sset_4, bc_func),
   ]
-  # update_dofs!(asm, dbcs)
 
-  # # pre-setup some scratch arrays
-  p = create_parameters(asm, physics, dbcs)
-  Uu = create_unknowns(asm)
+  # setup the parameters
+  p = create_parameters(asm, physics; dirichlet_bcs=dbcs)
 
-  # # # solver = NewtonSolver(DirectLinearSolver(asm))
-  solver = NewtonSolver(IterativeSolver(asm, :CgSolver))
-  # # # this call below updates the BCs...
-  # # # maybe we should change that name?
-  update_bcs!(H1Field, solver, Uu, p)
-
-  FiniteElementContainers.solve!(solver, Uu, p)
+  # setup solver and integrator
+  solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
+  integrator = QuasiStaticIntegrator(solver)
+  evolve!(integrator, p)
 
   write_times(pp, 1, 0.0)
   write_field(pp, 1, p.h1_field)
diff --git a/test/poisson/TestPoissonCUDA.jl b/test/poisson/TestPoissonCUDA.jl
index 1701c75..73f061b 100644
--- a/test/poisson/TestPoissonCUDA.jl
+++ b/test/poisson/TestPoissonCUDA.jl
@@ -48,12 +48,10 @@ end
     DirichletBC(:u, :sset_3, bc_func),
     DirichletBC(:u, :sset_4, bc_func_2),
   ]
-  # TODO this one will be tough to do on the GPU
-  # update_dofs!(asm, dbcs)
 
   # create parameters on CPU
   # TODO make a better constructor
-  p = create_parameters(asm, physics, dbcs)
+  p = create_parameters(asm, physics; dirichlet_bcs=dbcs)
   # need to assemble once before moving to GPU
   # TODO try to wrap this in the |> gpu call
   U = create_field(asm, H1Field)
@@ -64,11 +62,9 @@ end
   p_gpu = p |> gpu
   asm_gpu = asm |> gpu
 
-  solver = NewtonSolver(IterativeSolver(asm_gpu, :CgSolver))
-  Uu = create_unknowns(asm_gpu)
-  update_bcs!(H1Field, solver, Uu, p_gpu)
-  @time FiniteElementContainers.solve!(solver, Uu, p_gpu)
-  @time FiniteElementContainers.solve!(solver, Uu, p_gpu)
+  solver = NewtonSolver(IterativeLinearSolver(asm_gpu, :CgSolver))
+  integrator = QuasiStaticIntegrator(solver)
+  evolve!(integrator, p_gpu)
 
   # transfer to cpu to post-process
   p = p_gpu |> cpu