modifying physics to be block by block.

Author: cmhamel

examples/mechanics/Cube.jl

@@ -52,7 +52,9 @@ end
 function mechanics_test()
   mesh = UnstructuredMesh("cube.g")
   V = FunctionSpace(mesh, H1Field, Lagrange) 
-  physics = Mechanics(ThreeDimensional())
+  physics = (;
+    cube=Mechanics(ThreeDimensional())
+  )
 
   u = VectorFunction(V, :displ)
   asm = SparseMatrixAssembler(H1Field, u)

src/FiniteElementContainers.jl

@@ -143,6 +143,7 @@ include("physics/Physics.jl")
 include("PostProcessors.jl")
 
 #
+include("TimeSteppers.jl")
 include("Parameters.jl")
 include("solvers/Solvers.jl")
 include("integrators/Integrators.jl")

src/Parameters.jl

@@ -1,29 +1,3 @@
-# move to a seperate file
-abstract type AbstractTimeStepper{T} end
-current_time(t::AbstractTimeStepper) = sum(t.time_current)
-time_step(t::AbstractTimeStepper) = sum(t.Δt)
-
-struct TimeStepper{T} <: AbstractTimeStepper{T}
-  time_start::T
-  time_end::T
-  time_current::T
-  Δt::T
-end
-
-function TimeStepper(time_start_in::T, time_end_in::T, n_steps::Int) where T <: Number
-  time_start = zeros(1)
-  time_end = zeros(1)
-  time_current = zeros(1)
-  Δt = zeros(1)
-  Δt = zeros(1)
-  fill!(time_start, time_start_in)
-  fill!(time_end, time_end_in)
-  fill!(time_current, time_start_in)
-  fill!(Δt, (time_end_in - time_start_in) / n_steps)
-  return TimeStepper(time_start, time_end, time_current, Δt)
-  # return TimeStepper(time_start_in, time_end_in, )
-end
-
 abstract type AbstractParameters end
 
 struct Parameters{D, N, T, Phys, Props, S, V, H1} <: AbstractParameters
@@ -42,26 +16,13 @@ struct Parameters{D, N, T, Phys, Props, S, V, H1} <: AbstractParameters
   h1_field::H1
 end
 
-# TODO only works for H1Fields currently most likely
-# function Parameters(dof::DofManager, physics::AbstractPhysics)
-#   n_props = num_properties(physics)
-#   n_state = num_states(physics)
-
-#   # TODO
-#   fspace = dof.H1_vars[1].fspace
-
-#   for (name, ref_fe) in pairs(fspace.ref_fes)
-
-#   end
-# end
-
 # TODO 
 # 1. need to loop over bcs and vars in dof
 #    to make organize different bcs based on fspace type
 # 2. group all dbcs, nbcs of similar fspaces into a single struct?
-# 3. figure out how to handle function pointers on the GPU
+# 3. figure out how to handle function pointers on the GPU - done
 # 4. add different fspace types
-# 5. convert vectors of dbcs/nbcs into namedtuples
+# 5. convert vectors of dbcs/nbcs into namedtuples - done
 function Parameters(
   assembler, physics,
   dirichlet_bcs, 
@@ -76,6 +37,16 @@ function Parameters(
   state_old = nothing
   state_new = nothing
 
+  # for mixed spaces we'll need to do this more carefully
+  if isa(physics, AbstractPhysics)
+    syms = keys(values(assembler.dof.H1_vars)[1].fspace.elem_conns)
+    physics = map(x -> physics, syms)
+    physics = NamedTuple{tuple(syms...)}(tuple(physics...))
+  else
+    @assert isa(physics, NamedTuple)
+    # TODO re-arrange physics tuple to match fspaces when appropriate
+  end
+
   if dirichlet_bcs !== nothing
     syms = map(x -> Symbol("dirichlet_bc_$x"), 1:length(dirichlet_bcs))
     # dbcs = NamedTuple{tuple(syms...)}(tuple(dbcs...))

src/TimeSteppers.jl

@@ -0,0 +1,24 @@
+abstract type AbstractTimeStepper{T} end
+current_time(t::AbstractTimeStepper) = sum(t.time_current)
+time_step(t::AbstractTimeStepper) = sum(t.Δt)
+
+struct TimeStepper{T} <: AbstractTimeStepper{T}
+  time_start::T
+  time_end::T
+  time_current::T
+  Δt::T
+end
+
+function TimeStepper(time_start_in::T, time_end_in::T, n_steps::Int) where T <: Number
+  time_start = zeros(1)
+  time_end = zeros(1)
+  time_current = zeros(1)
+  Δt = zeros(1)
+  Δt = zeros(1)
+  fill!(time_start, time_start_in)
+  fill!(time_end, time_end_in)
+  fill!(time_current, time_start_in)
+  fill!(Δt, (time_end_in - time_start_in) / n_steps)
+  return TimeStepper(time_start, time_end, time_current, Δt)
+  # return TimeStepper(time_start_in, time_end_in, )
+end

src/assemblers/Assemblers.jl

@@ -38,26 +38,6 @@ function _assemble_element!(global_val::H1Field, local_val, conn, e, b)
   return nothing
 end
 
-# """
-# $(TYPEDSIGNATURES)
-# Assembly method for a block labelled as block_id.
-
-# This is a top level method that simply dispatches to CPU/GPU specializations
-# based on the backend of assembler and also checks for consistent backends.
-
-# This method is common acroos all value types, e.g. energy, residual, stiffness, etc.
-# to maintain a common method call syntax. The difference in dispatch is handled by the 
-# sym argument. 
-# """
-# function _assemble_block!(assembler, physics, sym, ref_fe, U, X, conns, block_id)
-#   backend = KA.get_backend(assembler)
-#   # TODO add get_backend method of ref_fe
-#   @assert backend == KA.get_backend(U)
-#   @assert backend == KA.get_backend(X)
-#   @assert backend == KA.get_backend(conns)
-#   _assemble_block!(assembler, physics, sym, ref_fe, U, X, conns, block_id, backend)
-# end
-
 _assemble_block_method_from_sym(::Val{:mass}) = _assemble_block_mass!
 _assemble_block_method_from_sym(::Val{:residual}) = _assemble_block_residual!
 _assemble_block_method_from_sym(::Val{:residual_and_stiffness}) = _assemble_block_residual_and_stiffness!
@@ -76,29 +56,23 @@ end
 $(TYPEDSIGNATURES)
 Top level assembly method for ```H1Field``` that loops over blocks and dispatches
 to appropriate kernels based on sym.
+
+TODO need to make sure at setup time that physics and elem_conns have the same
+values order. Otherwise, shenanigans.
 """
 function assemble!(assembler, physics, U::H1Field, sym)
   val_sym = Val{sym}()
   _assemble_block_method! = _assemble_block_method_from_sym(val_sym)
   _zero_storage(assembler, val_sym)
   fspace = assembler.dof.H1_vars[1].fspace
   X = fspace.coords
-  for (b, conns) in enumerate(values(fspace.elem_conns))
+  for (b, (conns, block_physics)) in enumerate(zip(values(fspace.elem_conns), values(physics)))
     ref_fe = values(fspace.ref_fes)[b]
     backend = _check_backends(assembler, U, X, conns)
-    _assemble_block_method!(assembler, physics, ref_fe, U, X, conns, b, backend)
+    _assemble_block_method!(assembler, block_physics, ref_fe, U, X, conns, b, backend)
   end
 end
 
-# wrapper that uses Uu and p
-# TODO only works on H1 fields right now
-function assemble!(assembler, physics, Uu, p, sym)
-  update_field_bcs!(p.h1_field, assembler.dof, p.h1_bcs, p.t)
-  update_field_unknowns!(p.h1_field, assembler.dof, Uu)
-  assemble!(assembler, physics, p.h1_field, sym)
-  return nothing
-end
-
 """
 $(TYPEDSIGNATURES)
 """
@@ -117,16 +91,6 @@ $(TYPEDSIGNATURES)
 """
 create_unknowns(asm::AbstractAssembler, type::Type{<:AbstractField}) = create_unknowns(asm.dof, type)
 
-"""
-$(TYPEDSIGNATURES)
-"""
-function _element_level_coordinates(X::H1Field, ref_fe, conns, e)
-  NDim = size(X, 1)
-  NNPE = ReferenceFiniteElements.num_vertices(ref_fe)
-  x_el = @views SMatrix{NDim, NNPE, Float64, NDim * NNPE}(X[:, conns[:, e]])
-  return x_el
-end
-
 """
 $(TYPEDSIGNATURES)
 """

src/assemblers/CPUGeneral.jl

@@ -38,7 +38,7 @@ function _assemble_block_residual!(assembler, physics, ref_fe, U, X, conns, bloc
   NNPE = ReferenceFiniteElements.num_vertices(ref_fe)
   NxNDof = NNPE * ND
   for e in axes(conns, 2)
-    x_el = _element_level_coordinates(X, ref_fe, conns, e)
+    x_el = _element_level_fields(X, ref_fe, conns, e)
     u_el = _element_level_fields(U, ref_fe, conns, e)
     R_el = zeros(SVector{NxNDof, eltype(assembler.residual_storage)})
 

src/assemblers/GPUGeneral.jl

@@ -11,7 +11,7 @@ KA.@kernel function _assemble_block_residual_kernel!(assembler, physics, ref_fe,
   NNPE = ReferenceFiniteElements.num_vertices(ref_fe)
   NxNDof = NNPE * ND
 
-  x_el = _element_level_coordinates(X, ref_fe, conns, E)
+  x_el = _element_level_fields(X, ref_fe, conns, E)
   u_el = _element_level_fields(U, ref_fe, conns, E)
   R_el = zeros(SVector{NxNDof, Float64})
 
@@ -52,7 +52,7 @@ KA.@kernel function _assemble_block_stiffness_kernel!(assembler, physics, ref_fe
   NNPE = ReferenceFiniteElements.num_vertices(ref_fe)
   NxNDof = NNPE * ND
 
-  x_el = _element_level_coordinates(X, ref_fe, conns, E)
+  x_el = _element_level_fields(X, ref_fe, conns, E)
   u_el = _element_level_fields(U, ref_fe, conns, E)
   K_el = zeros(SMatrix{NxNDof, NxNDof, Float64, NxNDof * NxNDof})
 

src/assemblers/SparseMatrixAssembler.jl

@@ -103,7 +103,7 @@ function _assemble_block_mass!(assembler, physics, ref_fe, U, X, conns, block_id
   NNPE = ReferenceFiniteElements.num_vertices(ref_fe)
   NxNDof = NNPE * ND
   for e in axes(conns, 2)
-    x_el = _element_level_coordinates(X, ref_fe, conns, e)
+    x_el = _element_level_fields(X, ref_fe, conns, e)
     u_el = _element_level_fields(U, ref_fe, conns, e)
     M_el = zeros(SMatrix{NxNDof, NxNDof, eltype(assembler.mass_storage), NxNDof * NxNDof})
 
@@ -130,7 +130,7 @@ function _assemble_block_stiffness!(assembler, physics, ref_fe, U, X, conns, blo
   NNPE = ReferenceFiniteElements.num_vertices(ref_fe)
   NxNDof = NNPE * ND
   for e in axes(conns, 2)
-    x_el = _element_level_coordinates(X, ref_fe, conns, e)
+    x_el = _element_level_fields(X, ref_fe, conns, e)
     u_el = _element_level_fields(U, ref_fe, conns, e)
     K_el = zeros(SMatrix{NxNDof, NxNDof, eltype(assembler.stiffness_storage), NxNDof * NxNDof})
 
@@ -158,7 +158,7 @@ function _assemble_block_residual_and_stiffness!(assembler, physics, ref_fe, U,
   NNPE = ReferenceFiniteElements.num_vertices(ref_fe)
   NxNDof = NNPE * ND
   for e in axes(conns, 2)
-    x_el = _element_level_coordinates(X, ref_fe, conns, e)
+    x_el = _element_level_fields(X, ref_fe, conns, e)
     u_el = _element_level_fields(U, ref_fe, conns, e)
     R_el = zeros(SVector{NxNDof, eltype(assembler.residual_storage)})
     K_el = zeros(SMatrix{NxNDof, NxNDof, eltype(assembler.stiffness_storage), NxNDof * NxNDof})

test/mechanics/TestMechanicsCUDA.jl

@@ -73,7 +73,7 @@ end
   p = create_parameters(asm, physics; dirichlet_bcs=dbcs, times=times)
 
   U = create_field(asm, H1Field)
-  @time assemble!(asm, physics, U, :stiffness)
+  @time assemble!(asm, p.physics, U, :stiffness)
   K = stiffness(asm)
 
   # move to device

test/poisson/TestPoissonCUDA.jl

@@ -55,7 +55,7 @@ end
   # need to assemble once before moving to GPU
   # TODO try to wrap this in the |> gpu call
   U = create_field(asm, H1Field)
-  @time assemble!(asm, physics, U, :stiffness)
+  @time assemble!(asm, p.physics, U, :stiffness)
   K = stiffness(asm)
 
   # device movement