Add matrix-free mode to SparseMatrixAssembler; drop dead H/D buffers

Project.toml

@@ -1,6 +1,6 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
-version = "0.13.1"
+version = "0.13.2"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
 
 [deps]

examples/app/src/MyApp.jl

@@ -19,16 +19,13 @@ function app_main(ARGS::Vector{String})
     #####################################
     # need to define some types
     #####################################
+    ET = ExodusDatabase{Int32, Int32, Int32, Float64}
     FT = AT.ExpressionFunction{Float64}
     SPT = FEC.CSCMatrix()
 
     #####################################
     # setup function space
     #####################################
-    # block_ids = FEC._setup_juliac_safe_block_to_ref_fe_id(sim.mesh)
-    # for n in 1:16
-    #     println(Core.stdout, "Block $n element type $(block_ids[n])")
-    # end
     V = FunctionSpace{true}(sim.mesh, H1Field, Lagrange)
     u = ScalarFunction(V, "u")
     dof = DofManager{false}(u)
@@ -44,53 +41,37 @@ function app_main(ARGS::Vector{String})
     end
 
     times = TimeStepper(0.0, 0.0, 1)
-    p = FEC.TypeStableParameters{FT}(sim.mesh, asm, physics, props, sim.ics, sim.dbcs, times)
-    # U = create_field(asm)
-    # Uu = create_unknowns(asm)
-
-    # # setting up ics
-    # # ics = InitialConditions{FT}(sim.mesh, dof, sim.ics)
-    # X = sim.mesh.nodal_coords
-    # # println(Core.stdout, "X size 1 = $(size(X, 1))")
-    # # println(Core.stdout, "X size 2 = $(size(X, 2))")
-    # # for block_id in 1:1
-    # #     ref_fe = FEC.block_reference_element(V, block_id)
-    # #     println(Core.stdout, ref_fe)
-    # # end
-
+    p = FEC.TypeStableParameters{FT}(
+        sim.mesh, asm,
+        physics, props,
+        sim.ics, sim.dbcs, sim.nbcs, sim.srcs, times
+    )
     FEC.initialize!(p)
-    # FEC._update_for_assembly!(p, dof, Uu)
 
-    # # lsolver = DirectLinearSolver(asm)
-    # # lsolver = x -> IterativeLinearSolver(x, :cg)
     preconditioner = I
     timer = TimerOutput()
     workspace = CgWorkspace(stiffness(asm), residual(asm))
     ΔUu = create_unknowns(asm)
     lsolver = IterativeLinearSolver(asm, preconditioner, workspace, timer, ΔUu)
     nlsolver = NewtonSolver(lsolver)
-
-    # println(Core.stdout, "Number of blocks = $(length(sim.mesh.element_types))")
-    # # N = length(sim.mesh.element_types)
-    # # # N = Val(N)
-    # # N = Val{N}()
-    # # ref_fe = _get_ref_fe(V, N)
-    # # indices = ntuple(x -> x, N)
-    # # indices = ntuple(x -> x, Val(N))
-    # # s = MyStruct{N}(N)::MyStruct{N}
     integrator = QuasiStaticIntegrator(nlsolver)
+
     println(sim.log_file.io, "Setup complete")
     println(sim.log_file.io, "Solving...")
     evolve!(integrator, p)
     println(sim.log_file.io, "Solve complete")
     println(Core.stdout, maximum(p.field.data))
     println(Core.stdout, minimum(p.field.data))
 
-    # pp = PostProcessor(sim.mesh, "juliac_output.e", u)
-
+    pp = PostProcessor{ET}(FEC.ExodusMesh, sim.mesh, "juliac_output.e")
+    FEC.add_function!(pp, u)
+    FEC.finalize_setup!(pp)
+    write_times(pp, 1, 0.0)
+    write_field(pp, 1, ("u",), p.field)
+    close(pp)
 end
 
 function @main(ARGS::Vector{String})
     app_main(ARGS)
     return 0
-end
+end

examples/expression/Project.toml

@@ -1,5 +1,6 @@
 name = "MyApp"
 
 [deps]
+DynamicExpressions = "a40a106e-89c9-4ca8-8020-a735e8728b6b"
 FiniteElementContainers = "d08262e4-672f-4e7f-a976-f2cea5767631"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

examples/expression/src/MyApp.jl

@@ -1,7 +1,11 @@
 import FiniteElementContainers: Expressions as Exp
 import FiniteElementContainers.Expressions: Parser
+using DynamicExpressions
 using StaticArrays
 
+const operators = OperatorEnum(1 => (sin, cos, exp), 2 => (+, -, *, /))
+
+
 function test_method()
     # constant
     println(Core.stdout, "Number")
@@ -54,6 +58,14 @@ function test_method()
     func = Exp.ExpressionFunction{Float64}(string, ["x", "y"])
     val = func([1.0, 2.0])
     println(Core.stdout, "func = $val")
+
+    # operators = OperatorEnum(1 => (sin, cos, exp), 2 => (+, -, *, /))
+    variable_names = ["x", "y"]
+    parsed_expr = parse_expression(
+        :(sin(2.0 * x + exp(y + 5.0))); 
+        operators = operators,
+        variable_names = variable_names
+    )
 end
 
 function @main(ARGS)

src/Fields.jl

@@ -102,7 +102,7 @@ struct Connectivity{
     nelems::Vector{T}
     offsets::Vector{T}
 
-    function Connectivity(mats::Vector{<:AbstractArray{<:Integer, 2}})
+    function Connectivity(mats::Vector{<:AbstractMatrix{<:Integer}})
         nblocks = length(mats)
         nepes = map(x -> size(x, 1), mats)
         nelems = map(x -> size(x, 2), mats)
@@ -135,6 +135,8 @@ function Adapt.adapt_structure(to, conn::Connectivity{T, D}) where {T, D}
     )
 end
 
+Base.eltype(::Connectivity{T, D}) where {T, D} = T
+
 # NOT GPU safe
 function connectivity(conn::Connectivity, b::Int)
     nepe = conn.nepes[b]

src/FunctionSpaces.jl

@@ -88,18 +88,15 @@ end
 
 # Lagrange elements
 # defaulting to fully integrated elements for now
-function _setup_juliac_safe_ref_fes(::Type{Lagrange}, q_type::Type{QT}) where {QT <: AbstractQuadratureType}
-  ref_fes = (
-    ReferenceFE(Hex{Lagrange, 1}(), QT(2, 2)),  # HEX8 for Lagrange
-    ReferenceFE(Quad{Lagrange, 1}(), QT(2, 2)), # QUAD4 for Lagrange
-    ReferenceFE(Quad{Lagrange, 2}(), QT(2, 2)), # QUAD9 for Lagrange
-    ReferenceFE(Tri{Lagrange, 1}(), QT(2, 2)),  # Tri3 for Lagrange
-    ReferenceFE(Tri{Lagrange, 2}(), QT(2, 2)),  # Tri3 for Lagrange
-    ReferenceFE(Tet{Lagrange, 1}(), QT(2, 2)),  # Tri3 for Lagrange
-    ReferenceFE(Tet{Lagrange, 2}(), QT(2, 2))   # Tri3 for Lagrange
-  )
-  return ref_fes
-end
+const _juliac_safe_ref_fes = (
+  ReferenceFE(Hex{Lagrange, 1}(), GaussLobattoLegendre(2, 2)),  # HEX8 for Lagrange
+  ReferenceFE(Quad{Lagrange, 1}(), GaussLobattoLegendre(2, 2)), # QUAD4 for Lagrange
+  ReferenceFE(Quad{Lagrange, 2}(), GaussLobattoLegendre(2, 2)), # QUAD9 for Lagrange
+  ReferenceFE(Tri{Lagrange, 1}(), GaussLobattoLegendre(2, 2)),  # Tri3 for Lagrange
+  ReferenceFE(Tri{Lagrange, 2}(), GaussLobattoLegendre(2, 2)),  # Tri3 for Lagrange
+  ReferenceFE(Tet{Lagrange, 1}(), GaussLobattoLegendre(2, 2)),  # Tri3 for Lagrange
+  ReferenceFE(Tet{Lagrange, 2}(), GaussLobattoLegendre(2, 2))   # Tri3 for Lagrange
+)
 
 """
 default code path that sets up ref fes as a namedtuple
@@ -167,10 +164,10 @@ struct FunctionSpace{
   RefFEs
 } <: AbstractFunctionSpace
   block_names::Vector{String}
-  # block_to_ref_fe_id::Vector{IT} # only used in juliac builds
   block_to_ref_fe_id::BTRE
   coords::Coords
   elem_conns::Connectivity{IT, IV}
+  # need to remove this mabye?
   elem_id_maps::Vector{Vector{IT}} # TODO create new type for ID map similar to connectivity
   ref_fes::RefFEs
 
@@ -209,13 +206,15 @@ function FunctionSpace{is_juliac_safe}(
     elseif p_degree == 0
       @assert false "TODO 0 order elements"
     else
+      # TODO all of this needs TLC
       ref_fes = _setup_ref_fes(mesh, interp_type, p_degree, q_type, q_degree)
       coords, conns = create_higher_order_mesh(mesh, H1Field, interp_type, p_degree)
       conns = Connectivity([val for val in values(conns)])
     end
   else
     if is_juliac_safe
-      ref_fes = _setup_juliac_safe_ref_fes(interp_type, q_type)
+      # ref_fes = _setup_juliac_safe_ref_fes(interp_type, q_type)
+      ref_fes = _juliac_safe_ref_fes
     else
       ref_fes = _setup_ref_fes(mesh, interp_type, nothing, q_type, q_degree)
     end
@@ -268,7 +267,8 @@ function Adapt.adapt_structure(to, fspace::FunctionSpace)
     fspace.block_names, fspace.block_to_ref_fe_id,
     adapt(to, fspace.coords),
     adapt(to, fspace.elem_conns), 
-    fspace.elem_id_maps,
+    # fspace.elem_id_maps,
+    map(x -> adapt(to, x), fspace.elem_id_maps),
     map(x -> adapt(to, x), fspace.ref_fes)
   )
 end
@@ -282,23 +282,54 @@ function Base.show(io::IO, fspace::FunctionSpace)
   end
 end
 
-function _is_juliac_safe(::FunctionSpace{B, I, V, C, R}) where {B, I, V, C, R}
+function _is_juliac_safe(::FunctionSpace{B, I, V, BTRE, C, R}) where {B, I, V, BTRE, C, R}
   return B
 end
 
 function block_entity_size(fspace::FunctionSpace, b::Int)
   return (num_entities_per_element(fspace, b), num_elements(fspace, b))
 end
 
-function block_reference_element(fspace::FunctionSpace{false, I, V, C, R}, block_id::Int) where {I, V, C, R}
+function block_reference_element(fspace::FunctionSpace{false, I, V, BTRE, C, R}, block_id::Int) where {I, V, BTRE, C, R}
   return fspace.ref_fes[block_id]
 end
 
-function block_reference_element(fspace::FunctionSpace{true, I, V, C, R}, block_id::Int) where {I, V, C, R}
-  ref_fe_id = fspace.block_to_ref_fe_id[block_id]
-  return fspace.ref_fes[ref_fe_id]
+@generated function block_reference_element(
+  fspace::FunctionSpace{true, IT, IV, BTRE, C, R},
+  block_id::Int
+) where {IT, IV, BTRE, C, R}
+
+  n_refs = length(BTRE.parameters)
+
+  cases = map(1:n_refs) do j
+      quote
+          if block_id == $j
+              return fspace.ref_fes[$j]
+          end
+      end
+  end
+
+  quote
+      id = fspace.block_to_ref_fe_id[block_id]  # runtime value (OK)
+
+      if id == -1
+          error("Inactive block ", block_id)
+      end
+
+      $(cases...)  # ONLY depends on compile-time j
+
+      error("Invalid ref_fe id ", id)
+  end
 end
 
+# function block_reference_element(
+#   fspace::FunctionSpace{true, IT, IV, BTRE, C, R},
+#   block_id::Int
+# ) where {IT, IV, BTRE, C, R}
+#   index = fspace.block_to_ref_fe_id[block_id]
+#   return _block_reference_element(fspace, Val{index}())
+# end
+
 # this does not work behind juliac
 function block_quadrature_size(fspace::FunctionSpace{false}, b::Int)
   ref_fe = block_reference_element(fspace, b)

src/Parameters.jl

@@ -111,26 +111,6 @@ function Parameters(
   end
 
   # setup state variables
-  # state_old = Array{Float64, 3}[]
-  # state_new = Array{Float64, 3}[]
-  # for (b, val) in enumerate(values(physics))
-  #   # create state variables for this block physics
-  #   NS = num_states(val)
-  #   NQ, NE = block_quadrature_size(fspace, b)
-
-  #   state_old_temp = zeros(NS, NQ, NE)
-  #   state_new_temp = zeros(NS, NQ, NE)
-  #   for e in 1:NE
-  #     for q in 1:NQ
-  #       state_old_temp[:, q, e] = create_initial_state(val)
-  #       state_new_temp[:, q, e] = create_initial_state(val)
-  #     end
-  #   end
-  #   push!(state_old, state_old_temp)
-  #   push!(state_new, state_new_temp)
-  # end
-  # state_old = L2Field(state_old)
-  # state_new = L2Field(state_new)
   state_old, state_new = _setup_state_variables(fspace, physics)
 
   # scratch
@@ -214,20 +194,17 @@ struct TypeStableParameters{
   RT     <: Number,
   IV     <: AbstractVector{IT},
   RV     <: AbstractVector{RT},
-  # RM1    <: AbstractMatrix,
-  # RM2    <: AbstractMatrix,
-  # RM3    <: AbstractMatrix,
-  # RM4    <: AbstractMatrix,
+  RM     <: AbstractMatrix{<:SVector},
   Phys,
   Props,
   Coords <: AbstractField,
   Field  <: AbstractField 
 } <: AbstractParameters
   ics::InitialConditions{Vector{InitialConditionFunction{FuncT}}, IV, RV}
   dirichlet_bcs::DirichletBCs{Vector{DirichletBCFunction{FuncT, FuncT, FuncT}}, IV, RV}
-  # neumann_bcs::NeumannBCs{NBCFuncs, IT, IV, RM1}
+  neumann_bcs::NeumannBCs{Vector{NeumannBCFunction{FuncT}}, IT, IV, RM}
   # robin_bcs::RobinBCs{RBCFuncs, IT, IV, RM2, RM3}
-  # sources::Sources{SRCFuncs, RM4}
+  sources::Sources{Vector{SourceFunction{FuncT}}, RM}
   times::TimeStepper{RT}
   physics::Phys
   properties::Props
@@ -239,11 +216,14 @@ struct TypeStableParameters{
   # scratch fields
   hvp_scratch_field::Field
 
-  function TypeStableParameters{F}(mesh, assembler, physics, props, ics, dbcs, times) where F
+  function TypeStableParameters{F}(mesh, assembler, physics, props, ics, dbcs, nbcs, srcs, times) where F
     dof = assembler.dof
+    ND = size(dof, 1)
     fspace = function_space(dof)
     ics = InitialConditions{F}(mesh, dof, ics)
     dbcs = DirichletBCs{F}(mesh, dof, dbcs)
+    nbcs = NeumannBCs{F}(mesh, dof, nbcs)
+    srcs = Sources{F}(mesh, dof, srcs)
 
     state_old, state_new = _setup_state_variables(fspace, physics)
 
@@ -256,27 +236,13 @@ struct TypeStableParameters{
     update_dofs!(assembler, dbcs)
 
     new{
-      F, Int, Float64, Vector{Int}, Vector{Float64},
+      F, Int, Float64, Vector{Int}, Vector{Float64}, Matrix{SVector{ND, Float64}},
       typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field)
     }(
-      ics, dbcs, times, 
+      ics, dbcs, nbcs, srcs,
+      times, 
       physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field
     )
-
-    # TODO need to settle on expected format for nbcs, rbcs, and sources...
-    # should we always expect a staticarrays even for a 1 dof problem?
-    # RM1 = Matrix{Float64}
-    # RM2 = Matrix{Float64}
-    # RM3 = Matrix{Float64}
-    # RM4 = Matrix{Float64}
-
-    # new{
-    #   F, Int, Float64, Vector{Int}, Vector{Float64}, RM1, RM2, RM3, RM4,
-    #   typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field)
-    # }(
-    #   ics, dbcs, times, 
-    #   physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field
-    # )
   end
 end
 
@@ -360,8 +326,8 @@ function update_bc_values!(p::TypeStableParameters, assembler)
   X = coordinates(p)
   t = current_time(p)
   update_bc_values!(p.dirichlet_bcs, X, t)
-  # update_bc_values!(p.neumann_bcs, assembler, X, t)
-  # update_source_values!(p.sources, assembler, X, t)
+  update_bc_values!(p.neumann_bcs, assembler, X, t)
+  update_source_values!(p.sources, assembler, X, t)
 
   # TODO how to handle Robin BCs?
   # currently assembly methods handle updating the field