moving properties down to block level.

Author: cmhamel

ext/FiniteElementContainersAdaptExt.jl

@@ -176,7 +176,7 @@ function Adapt.adapt_structure(to, p::FiniteElementContainers.Parameters)
     adapt(to, p.properties),
     adapt(to, p.state_old),
     adapt(to, p.state_new),
-    adapt(to, p.h1_dbcs),
+    adapt(to, p.h1_coords),
     adapt(to, p.h1_field)
   )
 end

src/Parameters.jl

@@ -1,7 +1,10 @@
 abstract type AbstractParameters end
 
 # TODO need to break up bcs to different field types
-struct Parameters{D, N, T, Phys, Props, S, V, H1} <: AbstractParameters
+struct Parameters{
+  D, N, T, Phys, Props, S, 
+  H1Coords, H1
+} <: AbstractParameters
   # actual parameter fields
   # TODO add boundary condition stuff and time stepping stuff
   dirichlet_bcs::D
@@ -13,7 +16,7 @@ struct Parameters{D, N, T, Phys, Props, S, V, H1} <: AbstractParameters
   state_old::S
   state_new::S
   # scratch fields
-  h1_dbcs::V # TODO remove this
+  h1_coords::H1Coords
   h1_field::H1
 end
 
@@ -30,7 +33,7 @@ function Parameters(
   neumann_bcs, 
   times
 )
-  h1_dbcs = create_bcs(assembler, H1Field)
+  h1_coords = assembler.dof.H1_vars[1].fspace.coords
   h1_field = create_field(assembler, H1Field)
 
   # TODO
@@ -47,8 +50,14 @@ function Parameters(
   end
 
   # state_old = Array{Float64, 3}[]
+  properties = []
   state_old = L2QuadratureField[]
   for (key, val) in pairs(physics)
+    # create properties for this block physics
+    # TODO specialize to allow for element level properties
+    push!(properties, create_properties(val))
+
+    # create state variables for this block physics
     NS = num_states(val)
     NQ = ReferenceFiniteElements.num_quadrature_points(
       getfield(values(assembler.dof.H1_vars)[1].fspace.ref_fes, key)
@@ -57,10 +66,21 @@ function Parameters(
       getfield(values(assembler.dof.H1_vars)[1].fspace.elem_conns, key),
       2
     )
+
     syms = tuple(map(x -> Symbol("state_variable_$x"), 1:NS)...)
-    state_old_temp = L2QuadratureField(zeros(NS, NQ, NE), syms)
+    state_old_temp = zeros(NS, NQ, NE)
+    for e in 1:NE
+      for q in 1:NQ
+        state_old_temp[:, q, e] = create_initial_state(val)
+      end
+    end
+    state_old_temp = L2QuadratureField(state_old_temp, syms)
+
     push!(state_old, state_old_temp)
   end
+
+  properties = NamedTuple{keys(physics)}(tuple(properties...))
+
   state_new = copy(state_old)
   state_old = NamedTuple{keys(physics)}(tuple(state_old...))
   state_new = NamedTuple{keys(physics)}(tuple(state_new...))
@@ -95,7 +115,7 @@ function Parameters(
     physics, 
     properties, 
     state_old, state_new, 
-    h1_dbcs, h1_field
+    h1_coords, h1_field
   )
 
   update_dofs!(assembler, p)
@@ -118,9 +138,7 @@ function create_parameters(
 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
 

src/assemblers/Assemblers.jl

@@ -46,60 +46,65 @@ function _check_backends(assembler, U, X, state_old, state_new, conns)
   @assert backend == KA.get_backend(conns)
   @assert backend == KA.get_backend(state_old)
   @assert backend == KA.get_backend(state_new)
+  # props will be complicated...
+  # TODO
   return backend
 end
 
 function assemble!(assembler, ::Type{H1Field}, p, val_sym::Val{:mass})
-  _zero_storage(assembler, val_sym)
   fspace = assembler.dof.H1_vars[1].fspace
-  for (b, (conns, block_physics, state_old, state_new)) in enumerate(zip(
+  _zero_storage(assembler, val_sym)
+  for (b, (conns, block_physics, state_old, state_new, props)) in enumerate(zip(
     values(fspace.elem_conns), 
     values(p.physics),
-    values(p.state_old), values(p.state_new)
+    values(p.state_old), values(p.state_new),
+    values(p.properties)
   ))
     ref_fe = values(fspace.ref_fes)[b]
-    backend = _check_backends(assembler, p.h1_field, fspace.coords, state_old, state_new, conns)
+    backend = _check_backends(assembler, p.h1_field, p.h1_coords, state_old, state_new, conns)
     _assemble_block_mass!(
       assembler, block_physics, ref_fe, 
-      p.h1_field, fspace.coords, state_old, state_new,
+      p.h1_field, p.h1_coords, state_old, state_new, props,
       conns, b, 
       backend
     )
   end
 end
 
 function assemble!(assembler, ::Type{H1Field}, p, val_sym::Val{:residual})
-  _zero_storage(assembler, val_sym)
   fspace = assembler.dof.H1_vars[1].fspace
-  for (b, (conns, block_physics, state_old, state_new)) in enumerate(zip(
+  _zero_storage(assembler, val_sym)
+  for (b, (conns, block_physics, state_old, state_new, props)) in enumerate(zip(
     values(fspace.elem_conns), 
     values(p.physics),
-    values(p.state_old), values(p.state_new)
+    values(p.state_old), values(p.state_new),
+    values(p.properties)
   ))
     ref_fe = values(fspace.ref_fes)[b]
-    backend = _check_backends(assembler, p.h1_field, fspace.coords, state_old, state_new, conns)
+    backend = _check_backends(assembler, p.h1_field, p.h1_coords, state_old, state_new, conns)
     _assemble_block_residual!(
       assembler, block_physics, ref_fe, 
-      p.h1_field, fspace.coords, state_old, state_new,
+      p.h1_field, p.h1_coords, state_old, state_new, props,
       conns, b, 
       backend
     )
   end
 end
 
 function assemble!(assembler, ::Type{H1Field}, p, val_sym::Val{:residual_and_stiffness})
-  _zero_storage(assembler, val_sym)
   fspace = assembler.dof.H1_vars[1].fspace
-  for (b, (conns, block_physics, state_old, state_new)) in enumerate(zip(
+  _zero_storage(assembler, val_sym)
+  for (b, (conns, block_physics, state_old, state_new, props)) in enumerate(zip(
     values(fspace.elem_conns), 
     values(p.physics),
-    values(p.state_old), values(p.state_new)
+    values(p.state_old), values(p.state_new),
+    values(p.properties)
   ))
     ref_fe = values(fspace.ref_fes)[b]
-    backend = _check_backends(assembler, p.h1_field, fspace.coords, state_old, state_new, conns)
+    backend = _check_backends(assembler, p.h1_field, p.h1_coords, state_old, state_new, conns)
     _assemble_block_residual_and_stiffness!(
       assembler, block_physics, ref_fe, 
-      p.h1_field, fspace.coords, state_old, state_new,
+      p.h1_field, p.h1_coords, state_old, state_new, props,
       conns, b, 
       backend
     )
@@ -109,16 +114,17 @@ end
 function assemble!(assembler, ::Type{H1Field}, p, val_sym::Val{:stiffness})
   _zero_storage(assembler, val_sym)
   fspace = assembler.dof.H1_vars[1].fspace
-  for (b, (conns, block_physics, state_old, state_new)) in enumerate(zip(
+  for (b, (conns, block_physics, state_old, state_new, props)) in enumerate(zip(
     values(fspace.elem_conns), 
     values(p.physics),
-    values(p.state_old), values(p.state_new)
+    values(p.state_old), values(p.state_new),
+    values(p.properties)
   ))
     ref_fe = values(fspace.ref_fes)[b]
-    backend = _check_backends(assembler, p.h1_field, fspace.coords, state_old, state_new, conns)
+    backend = _check_backends(assembler, p.h1_field, p.h1_coords, state_old, state_new, conns)
     _assemble_block_stiffness!(
       assembler, block_physics, ref_fe, 
-      p.h1_field, fspace.coords, state_old, state_new,
+      p.h1_field, p.h1_coords, state_old, state_new, props,
       conns, b, 
       backend
     )

src/assemblers/CPUGeneral.jl

@@ -35,7 +35,7 @@ TODO remove Float64 typing below for eventual unitful use
 """
 function _assemble_block_residual!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id, ::KA.CPU
 )
   ND = size(U, 1)

src/assemblers/GPUGeneral.jl

@@ -6,7 +6,7 @@ Kernel for residual block assembly
 """
 KA.@kernel function _assemble_block_residual_kernel!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id
 )
   E = KA.@index(Global)
@@ -45,7 +45,7 @@ TODO add state variables and physics properties
 """
 function _assemble_block_residual!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id, backend::KA.Backend
 )
   kernel! = _assemble_block_residual_kernel!(backend)
@@ -59,7 +59,7 @@ end
 
 KA.@kernel function _assemble_block_stiffness_kernel!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id
 )
   E = KA.@index(Global)
@@ -98,7 +98,7 @@ TODO add state variables and physics properties
 """
 function _assemble_block_stiffness!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id, backend::KA.Backend
 )
   kernel! = _assemble_block_stiffness_kernel!(backend)

src/assemblers/SparseMatrixAssembler.jl

@@ -100,7 +100,7 @@ TODO add state variables and physics properties
 """
 function _assemble_block_mass!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id, ::KA.CPU
 )
   ND = size(U, 1)
@@ -131,7 +131,7 @@ TODO add state variables and physics properties
 """
 function _assemble_block_stiffness!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id, ::KA.CPU
 )
   ND = size(U, 1)
@@ -163,7 +163,7 @@ TODO remove Float64 typing below for eventual unitful use
 """
 function _assemble_block_residual_and_stiffness!(
   assembler, physics, ref_fe, 
-  U, X, state_old, state_new,
+  U, X, state_old, state_new, props,
   conns, block_id, ::KA.CPU
 )
   ND = size(U, 1)

src/integrators/Integrators.jl

@@ -4,11 +4,36 @@ abstract type AbstractIntegrator{
 } end
 # abstract type AbstractFirstOrderIntegrator <: AbstractIntegrator end
 # abstract type AbstractSecondOrderIntegrator <: AbstractIntegrator end
+abstract type AbstractSecondOrderIntegrator{
+  Solution,
+  Solver
+} <: AbstractIntegrator{Solution, Solver} end
 abstract type AbstractStaticIntegrator{
   Solution,
   Solver
 } <: AbstractIntegrator{Solution, Solver} end
 
 function evolve! end
 
+# TODO maybe move these methods below to BCs
+function _update_bcs!(bc, U, ::KA.CPU)
+  for (dof, val) in zip(bc.bookkeeping.dofs, bc.vals)
+    U[dof] = val
+  end
+  return nothing
+end
+
+KA.@kernel function _update_bcs_kernel!(bc, U)
+  I = KA.@index(Global)
+  dof = bc.bookkeeping.dofs[I]
+  val = bc.vals[I]
+  U[dof] = val
+end
+
+function _update_bcs!(bc, U, backend::KA.Backend)
+  kernel! = _update_bcs_kernel!(backend)
+  kernel!(bc, U, ndrange=length(bc.vals))
+  return nothing
+end
+
 include("QuasiStaticIntegrator.jl")

src/integrators/QuasiStaticIntegrator.jl

@@ -14,7 +14,20 @@ end
 
 function evolve!(integrator::QuasiStaticIntegrator, p)
   update_time!(p)
-  update_bcs!(H1Field, integrator.solver, integrator.solution, p)
+  # update_bcs!(H1Field, integrator.solver, p)
+  update_bcs!(integrator, p)
   solve!(integrator.solver, integrator.solution, p)
   return nothing
 end
+
+function update_bcs!(::QuasiStaticIntegrator, p::Parameters)
+  X = p.h1_coords # TODO won't work for mixed bcs
+  t = current_time(p.times)
+
+  update_bc_values!(p.dirichlet_bcs, X, t)
+
+  for bc in values(p.dirichlet_bcs)
+    _update_bcs!(bc, p.h1_field, KA.get_backend(bc))
+  end
+  return nothing
+end

src/physics/Physics.jl

@@ -4,38 +4,34 @@ num_fields(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS} = NF
 num_properties(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS} = NP
 num_states(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS} = NS
 
-# this won't work
-# function create_physics(name, vars...)
-#   # @show field_vars
-#   # @show prop_vars
-#   H1_vars = _filter_field_type(vars, H1Field)
-#   # TODO setup Hdiv and Hcurl logic as well
-#   L2_element_vars = _filter_field_type(vars, L2ElementField)
-#   L2_quadrature_vars = _filter_field_type(vars, L2QuadratureField)
-  
-#   H1_var_names = ()
-#   for var in H1_vars
-#     H1_var_names = (H1_var_names..., names(var)...)
-#   end
-#   H1_var_names = NamedTuple{H1_var_names}(1:length(H1_var_names))
+function create_properties(physics::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS}
+  @assert false "You need to implement the create_properties method for physics $(physics) or 
+  type $(typeof(physics))!"
+end
 
-#   L2_element_var_names = ()
-#   for var in L2_element_vars
-#     L2_element_var_names = (L2_element_var_names..., names(var)...)
-#   end
-#   L2_element_var_names = NamedTuple{L2_element_var_names}(1:length(L2_element_var_names))
+function create_properties(::AbstractPhysics{NF, 0, NS}) where {NF, NS}
+  return SVector{0, Float64}()
+end
 
-#   L2_quadrature_var_names = ()
-#   for var in L2_quadrature_vars
-#     L2_quadrature_var_names = (L2_quadrature_var_names..., names(var)...)
-#   end
-#   L2_quadrature_var_names = NamedTuple{L2_quadrature_var_names}(1:length(L2_quadrature_var_names))
-  
-# end
+function create_initial_state(::AbstractPhysics{NF, NP, 0}) where {NF, NP}
+  return SVector{0, Float64}()
+end
 
-# function initialize_state(::AbstractPhysics{NF, NP, 0}) where {NF, NP}
-#   return 
-# end
+# Can we make something that makes interfacing with kernels easier?
+# How can we make something like this work nicely with AD?
+# struct PhysicsQuadratureState{T, ND, NN, NF, NP, NS, NDxNF, NNxND, NNxNNxND}
+#   # u::SVector{NF, T}
+#   # ∇u::SMatrix{NF, ND, T, NDxNF}
+#   # props::SVector{NP, T}
+#   # state_old::SVector{NS, T}
+#   # interpolants and gauss weight at quadrature point
+#   N::SVector{NN, T}
+#   ∇N_ξ::SMatrix{NN, ND, T, NNxND}
+#   ∇∇N_ξ::SArray{Tuple{NN, ND, ND}, T, 3, NNxNNxND}
+#   JxW::T
+#   # element level fields
+#   u_el::SMatrix{}
+# end 
 
 # physics like methods
 function damping end