Assembler bug fix. Addresses #315

src/AppTools.jl

@@ -15,6 +15,7 @@ import ..NeumannBC
 import ..PeriodicBC
 import ..RobinBC
 import ..Source
+import ..TimeStepper
 import ..UnstructuredMesh
 import ..element_blocks
 import ..element_ids
@@ -285,6 +286,7 @@ print_dict(l::LogFile, d::Dict{String, String}) = print_dict(l, d)
 # Input file
 #######################################################
 struct FunctionSettings{N, T <: Number}
+    dict::Dict{String, Any}
     scalar_expr_funcs::Dict{String, ScalarExpressionFunction{T}}
     vector_expr_funcs::Dict{String, VectorExpressionFunction{N, T}}
 
@@ -321,12 +323,15 @@ struct FunctionSettings{N, T <: Number}
                     @assert false "Unsupported function type $type"
                 end
             end
+        else
+            func_settings = Dict{String, Any}()
         end
-        return new{N, T}(scalar_functions, vector_functions)
+        return new{N, T}(func_settings, scalar_functions, vector_functions)
     end
 end
 
 struct BCSettings{N, T <: Number}
+    dict::Dict{String, Any}
     dirichlet::Vector{DirichletBC{ScalarExpressionFunction{T}}}
     neumann::Vector{NeumannBC{VectorExpressionFunction{N, T}}}
     periodic::Vector{PeriodicBC{ScalarExpressionFunction{T}}}
@@ -435,17 +440,10 @@ struct BCSettings{N, T <: Number}
             end
         end
 
-        new{N, T}(dbcs, nbcs, pbcs, rbcs, srcs)
+        new{N, T}(bc_settings, dbcs, nbcs, pbcs, rbcs, srcs)
     end
 end
 
-
-# struct FunctionSpaceSettings
-#     field_type::String
-#     polynomial_type::String
-# end
-
-# TODO currently only supports blocks
 struct ICSettings{T <: Number}
     ics::Vector{InitialCondition{ScalarExpressionFunction{T}}}
 
@@ -489,32 +487,65 @@ struct ICSettings{T <: Number}
     end
 end
 
+struct MaterialSettings
+    materials::Dict{String, Any}
+
+    function MaterialSettings(log_file, parser)
+        if haskey(parser, "materials")
+            materials = InputFileParser.get_nested_block(parser, "materials")
+        else
+            materials = Dict{String, Any}()
+        end
+        new(materials)
+    end
+end
+
 struct MeshSettings
     dimension::Int
     file_path::String
     file_type::String
 
-    function MeshSettings(log_file, data)
-        mesh_settings = data["mesh"]::Dict{String, Any}
+    function MeshSettings(log_file, parser)
+        mesh_settings = InputFileParser.get_nested_block(parser, "mesh")
         dimension = mesh_settings["dimension"]::Int
         file_path = mesh_settings["file path"]::String
         file_type = lowercase(mesh_settings["file type"]::String)
         new(dimension, file_path, file_type)
     end
 end
 
-# struct VariableSettings
-#     function_space::String
-#     name::String
-#     variable_type::String
-# end
+struct TimeSettings{T <: Number}
+    dict::Dict{String, Any}
+    time::Union{Nothing, TimeStepper{T}}
+
+    function TimeSettings(log_file, parser)
+        if haskey(parser, "time")
+            time_settings = InputFileParser.get_nested_block(parser, "time")
+            end_time = time_settings["end time"]::Float64
+            start_time = time_settings["start time"]::Float64
+            if haskey(time_settings, "number of time steps")
+                nt = time_settings["number of time steps"]::Int
+                Δt = (end_time - start_time) / nt
+            elseif haskey(time_settings, "time step")
+                Δt = time_settings["time step"]::Float64
+            else
+                @assert false "Requires either one of \"number of time steps\" or \"time step\" commands."
+            end
+            new{Float64}(time_settings, TimeStepper(start_time, end_time, start_time, Δt))
+        else
+            new{Float64}(Dict{String, Any}(), nothing)
+        end
+    end
+end
 
 struct InputSettings{N, T <: Number}
     bcs::BCSettings{N, T}
     functions::FunctionSettings{N, T}
     ics::ICSettings{T}
+    materials::MaterialSettings
     mesh::MeshSettings
     parser::InputFileParser.Parser
+    time::TimeSettings{T}
 end
 
 function InputSettings{N}(cli_args::CLIArgParser, log_file::LogFile, ::Type{T} = Float64) where {N, T <: Number}
@@ -532,35 +563,11 @@ function InputSettings{N}(cli_args::CLIArgParser, log_file::LogFile, ::Type{T} =
     functions = FunctionSettings{N, T}(log_file, parser)
     bcs = BCSettings{N, T}(log_file, parser, functions)
     ics = ICSettings{T}(log_file, parser, functions)
+    materials = MaterialSettings(log_file, parser)
     mesh = MeshSettings(log_file, parser)
-    return InputSettings{N, T}(bcs, functions, ics, mesh, parser)
-end
-
-# function _parse_function_space_settings(log_file, data)
-#     settings = data["function_spaces"]::Dict{String, Any}
-#     fspaces = Dict{String, FunctionSpaceSettings}()
-#     for (k, v) in settings
-#         name = k::String
-#         temp = v::Dict{String, Any}
-#         field_type = temp["field_type"]::String
-#         polynomial_type = temp["polynomial_type"]::String
-#         fspaces[name] = FunctionSpaceSettings(field_type, polynomial_type)
-#     end
-#     return fspaces
-# end
-
-# function _parse_variable_settings(log_file, data)
-#     var_settings = data["variables"]::Dict{String, Any}
-#     vars = Dict{String, VariableSettings}()
-#     for (k, v) in var_settings
-#         name = k::String
-#         temp = v::Dict{String, Any}
-#         fspace = temp["function_space"]::String
-#         type = temp["type"]::String
-#         vars[name] = VariableSettings(fspace, name, type)
-#     end
-#     return vars
-# end
+    time = TimeSettings(log_file, parser)
+    return InputSettings{N, T}(bcs, functions, ics, materials, mesh, parser, time)
+end
 
 #######################################################
 # MeshIO strongly typed helpers
@@ -656,6 +663,7 @@ end
 struct Simulation{D, N, T <: Number, IO, Mesh}
     dbcs::Vector{DirichletBC{ScalarExpressionFunction{T}}}
     ics::Vector{InitialCondition{ScalarExpressionFunction{T}}}
+    input_settings::InputSettings{N, T}
     log_file::LogFile{IO}
     mesh::Mesh
     nbcs::Vector{NeumannBC{VectorExpressionFunction{N, T}}}
@@ -676,12 +684,13 @@ struct Simulation{D, N, T <: Number, IO, Mesh}
         #     println(log_file.io, fspace::FunctionSpace)
         # end
         if length(settings.ics.ics) > 1
-            T = eltype(settings.ics[1])
+            T = eltype(settings.ics.ics[1])
         else
             T = Float64
         end
         new{D, N, T, IO, typeof(mesh)}(
-            settings.bcs.dirichlet, settings.ics.ics, log_file, mesh,
+            settings.bcs.dirichlet, settings.ics.ics, 
+            settings, log_file, mesh,
             settings.bcs.neumann, settings.bcs.periodic,
             settings.bcs.robin, settings.bcs.source
         )

src/Parameters.jl

@@ -253,6 +253,36 @@ struct TypeStableParameters{
       physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field
     )
   end
+
+  function TypeStableParameters{SF, VF}(
+    mesh, assembler, physics, props, state_old, state_new,
+    ics, dbcs, nbcs, pbcs, srcs, times
+  ) where {SF, VF}
+    dof = assembler.dof
+    ND = size(dof, 1)
+    ics = InitialConditions{SF}(mesh, dof, ics)
+    dbcs = DirichletBCs{SF}(mesh, dof, dbcs)
+    nbcs = NeumannBCs{VF}(mesh, dof, nbcs)
+    pbcs = PeriodicBCs{SF}(mesh, dof, pbcs)
+    srcs = Sources{VF}(mesh, dof, srcs)
+
+    coords = mesh.nodal_coords
+    field = create_field(assembler)
+    field_old = create_field(assembler)
+    hvp_scratch_field = create_field(assembler)
+
+    # update assembler, where should this really live?
+    update_dofs!(assembler, dbcs, pbcs)
+
+    new{
+      SF, VF, Int, Float64, Vector{Int}, Vector{Float64}, Matrix{SVector{ND, Float64}},
+      typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field)
+    }(
+      ics, dbcs, nbcs, pbcs, srcs,
+      times, 
+      physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field
+    )
+  end
 end
 
 function create_parameters(

src/assemblers/Assemblers.jl

@@ -250,14 +250,31 @@ function _quadrature_level_state(state::AbstractArray{<:Number, 3}, q::Int, e::I
   return state_q
 end
 
-function _sparse_matrix(asm::AbstractAssembler, storage)
+function _sparse_matrix_mass(asm::AbstractAssembler, coo_storage)
   type = _sparse_matrix_type(asm)
   if isa(type, COOMatrix)
-    return _coo_matrix(asm.matrix_pattern, asm.dof, storage)
+    return _coo_matrix(asm.matrix_pattern, asm.dof, coo_storage)
   elseif isa(type, CSCMatrix)
-    return _csc_matrix(asm.matrix_pattern, asm.dof, storage)
+    csc_storage = asm.matrix_pattern.cscnzval_mass
+    return _csc_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage)
   elseif isa(type, CSRMatrix)
-    return _csr_matrix(asm.matrix_pattern, asm.dof, storage)
+    csc_storage = asm.matrix_pattern.cscnzval_mass
+    return _csr_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage)
+  else
+    @assert false "Should never happen. Handle this case better. Type is $(_sparse_matrix_type(asm))"
+  end
+end
+
+function _sparse_matrix_stiffness(asm::AbstractAssembler, coo_storage)
+  type = _sparse_matrix_type(asm)
+  if isa(type, COOMatrix)
+    return _coo_matrix(asm.matrix_pattern, asm.dof, coo_storage)
+  elseif isa(type, CSCMatrix)
+    csc_storage = asm.matrix_pattern.cscnzval_stiffness
+    return _csc_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage)
+  elseif isa(type, CSRMatrix)
+    csc_storage = asm.matrix_pattern.cscnzval_stiffness
+    return _csr_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage)
   else
     @assert false "Should never happen. Handle this case better. Type is $(_sparse_matrix_type(asm))"
   end
@@ -314,7 +331,7 @@ function mass(asm::AbstractAssembler)
     return _zero_sparse_matrix(asm)
   end
 
-  M = _sparse_matrix(asm, asm.mass_storage)
+  M = _sparse_matrix_mass(asm, asm.mass_storage)
 
   if _is_condensed(asm.dof)
     _adjust_matrix_entries_for_constraints!(M, asm.constraint_storage)
@@ -361,7 +378,7 @@ function stiffness(asm::AbstractAssembler)
     return _zero_sparse_matrix(asm)
   end
 
-  K = _sparse_matrix(asm, asm.stiffness_storage)
+  K = _sparse_matrix_stiffness(asm, asm.stiffness_storage)
 
   if _is_condensed(asm.dof)
     _adjust_matrix_entries_for_constraints!(K, asm.constraint_storage)

src/assemblers/Diagonal.jl

@@ -14,11 +14,26 @@ gives the true diagonal, whereas the row-sum approximation `K·1` can
 be zero at interior nodes of uniform meshes.
 """
 function assemble_diagonal!(
-  assembler, func::F, Uu, p
+  assembler::AbstractAssembler, func::F, u, p
 ) where F <: Function
-  storage = assembler.residual_storage
+  assemble_diagonal!(
+    assembler.residual_storage,
+    assembler.vector_pattern, assembler.dof,
+    func, u, p;
+    use_inplace_methods = _use_inplace_methods(assembler)
+  )
+end
+
+# function assemble_diagonal!(
+#   assembler, func::F, Uu, p
+# ) where F <: Function
+  # storage = assembler.residual_storage
+function assemble_diagonal!(
+  storage, pattern, dof, func, Uu, p;
+  use_inplace_methods::Bool = false
+)
   fill!(storage, zero(eltype(storage)))
-  dof = assembler.dof
+  # dof = assembler.dof
   fspace = function_space(dof)
   X = coordinates(p)
   t = current_time(p)
@@ -29,16 +44,28 @@ function assemble_diagonal!(
   return_type = AssembledDiagonal()
   conns = fspace.elem_conns
   foreach_block(fspace, p) do physics, props, ref_fe, b
-    _assemble_block!(
-      storage,
-      conns.data, conns.offsets[b],
-      func,
-      physics, ref_fe,
-      X, t, Δt,
-      U, U_old,
-      block_view(p.state_old, b), block_view(p.state_new, b), props,
-      return_type
-    )
+    if use_inplace_methods
+      _assemble_block!(
+        storage,
+        func,
+        physics,
+        t, Δt,
+        props,
+        block_view(p.state_old, b), block_view(p.state_new, b),
+        conns.data, conns.offsets[b], ref_fe, X, U, U_old
+      )
+    else
+      _assemble_block!(
+        storage,
+        conns.data, conns.offsets[b],
+        func,
+        physics, ref_fe,
+        X, t, Δt,
+        U, U_old,
+        block_view(p.state_old, b), block_view(p.state_new, b), props,
+        return_type
+      )
+    end
   end
   return nothing
 end

src/assemblers/SparseMatrixAssembler.jl

@@ -140,15 +140,16 @@ function _empty_matrix_pattern(dof::DofManager)
     similar(dof.unknown_dofs, 0),  # Is
     similar(dof.unknown_dofs, 0),  # Js
     similar(dof.unknown_dofs, 0),  # unknown_dofs
-    Int[],                          # block_start_indices
-    [0],                            # max_entries  (single zero — see create_assembler_cache)
+    Int[],                         # block_start_indices
+    [0],                           # max_entries  (single zero — see create_assembler_cache)
     similar(dof.unknown_dofs, 0),  # klasttouch
     similar(dof.unknown_dofs, 0),  # csrrowptr
     similar(dof.unknown_dofs, 0),  # csrcolval
-    Float64[],                      # csrnzval
+    Float64[],                     # csrnzval
     similar(dof.unknown_dofs, 0),  # csccolptr
     similar(dof.unknown_dofs, 0),  # cscrowval
-    Float64[],                      # cscnzval
+    Float64[],                     # cscnzval_mass
+    Float64[],                     # cscnzval_stiffness
     similar(dof.unknown_dofs, 0)   # permutation
   )
 end