Electromechanics

Project.toml

@@ -38,6 +38,9 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 
+[sources]
+OrdinaryDiffEqOperatorSplitting = {url = "https://github.com/SciML/OrdinaryDiffEqOperatorSplitting.jl.git", rev = "main"}
+
 [extensions]
 CuThunderboltExt = "CUDA"
 
@@ -55,7 +58,7 @@ Logging = "1.10"
 ModelingToolkit = "9"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqLowOrderRK = "1.2.0"
-OrdinaryDiffEqOperatorSplitting = "0.1"
+OrdinaryDiffEqOperatorSplitting = "0.2"
 Preferences = "1"
 SciMLBase = "2"
 UnPack = "1"

docs/Project.toml

@@ -14,3 +14,6 @@ OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Thunderbolt = "909927c2-98d5-4a67-bba9-79f03a9ad49b"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+
+[sources]
+OrdinaryDiffEqOperatorSplitting = {url = "https://github.com/SciML/OrdinaryDiffEqOperatorSplitting.jl.git", rev = "main"}

src/Thunderbolt.jl

@@ -223,6 +223,9 @@ export
     # Coupling
     Coupling,
     CoupledModel,
+    ElectroMechanicalCoupledModel,
+    ElectroMechanicalCoupler,
+    ElectroMechanicalSynchronizer,
     # Discretization
     semidiscretize,
     FiniteElementDiscretization,

src/discretization/fem.jl

@@ -219,3 +219,19 @@ function semidiscretize(model::QuasiStaticModel, discretization::FiniteElementDi
 
     return semidiscrete_problem
 end
+
+function Thunderbolt.semidiscretize(coupled_model::ElectroMechanicalCoupledModel, discretizations::Tuple, meshes::Tuple{<:Thunderbolt.AbstractGrid, <:Thunderbolt.AbstractGrid})
+    ep_sdp = semidiscretize(coupled_model.ep_model, discretizations[1], meshes[1])
+    mech_sdp = semidiscretize(coupled_model.structural_model, discretizations[2], meshes[2])
+    dh_mech_displacement = mech_sdp.dh
+    dh_ep_φₘ = ep_sdp.functions[1].dh
+    dh_mech_φₘ = create_dh_helper(get_grid(dh_mech_displacement), discretizations[1], :φₘ)
+    ep_mech_transfer_op = Thunderbolt.NodalIntergridInterpolation(dh_ep_φₘ, dh_mech_φₘ)
+    mech_ep_transfer_op = Thunderbolt.NodalIntergridInterpolation(dh_mech_φₘ, dh_ep_φₘ)
+    sync = ElectroMechanicalSynchronizer(ep_mech_transfer_op, mech_ep_transfer_op)
+    return GenericSplitFunction(
+        (ep_sdp, mech_sdp),
+        (Thunderbolt.OS.get_solution_indices(ep_sdp, 2), (last(Thunderbolt.OS.get_solution_indices(ep_sdp, 2))+1):((last(Thunderbolt.OS.get_solution_indices(ep_sdp, 2))) + solution_size(mech_sdp))),
+        (OS.NoExternalSynchronization(), sync)
+    )
+end

src/ferrite-addons/transfer_operators.jl

@@ -50,7 +50,7 @@ struct NodalIntergridInterpolation{PH <: PointEvalHandler, DH1 <: AbstractDofHan
     field_name_from::Symbol
     field_name_to::Symbol
 
-    function NodalIntergridInterpolation(dh_from::DofHandler{sdim}, dh_to::DofHandler{sdim}, field_name_from::Symbol, field_name_to::Symbol; subdomains_from = 1:length(dh_from.subdofhandlers), subdomains_to = 1:length(dh_to.subdofhandlers)) where sdim
+    function NodalIntergridInterpolation(dh_from::DofHandler{sdim}, dh_to::DofHandler{sdim}, field_name_from::Symbol, field_name_to::Symbol; subdomains_from = 1:length(dh_from.subdofhandlers), subdomains_to = 1:length(dh_to.subdofhandlers), search_nneighbors = 3, strategy = Ferrite.NewtonLineSearchPointFinder(residual_tolerance=1e-6)) where sdim
         @assert field_name_from ∈ Ferrite.getfieldnames(dh_from)
         @assert field_name_to ∈ Ferrite.getfieldnames(dh_to)
 
@@ -94,7 +94,7 @@ struct NodalIntergridInterpolation{PH <: PointEvalHandler, DH1 <: AbstractDofHan
             _compute_dof_nodes_barrier!(nodes, sdh, Ferrite.dof_range(sdh, field_name_to), gip, dof_to_node_map, ref_coords)
         end
 
-        ph = PointEvalHandler(Ferrite.get_grid(dh_from), nodes; warn=false)
+        ph = PointEvalHandler(Ferrite.get_grid(dh_from), nodes; search_nneighbors = search_nneighbors, strategy = strategy, warn=false)
 
         n_missing = sum(x -> x === nothing, ph.cells)
         n_missing == 0 || @warn "Constructing the interpolation for $field_name_from to $field_name_to failed. $n_missing (out of $(length(ph.cells))) points not found."

src/modeling/coupler/electromechanics.jl

@@ -0,0 +1,130 @@
+function create_dh_helper(mesh, discretization, sym)
+    ipc = _get_interpolation_from_discretization(discretization, sym)
+    dh = DofHandler(mesh)
+    for name in discretization.subdomains
+        add_subdomain!(dh, name, [ApproximationDescriptor(sym, ipc)])
+    end
+    close!(dh)
+    return dh
+end
+
+struct ElectroMechanicalCoupler <: Thunderbolt.VolumeCoupler 
+    ep_index::Int
+    mech_index::Int
+end
+
+struct ElectroMechanicalCoupledModel{SM #=<: QuasiStaticModel =#, EM, CT <: ElectroMechanicalCoupler}
+    ep_model::EM
+    structural_model::SM
+    coupler::CT
+end
+
+struct ElectroMechanicalSynchronizerCache{T<:Real, Dim, ∇uType<:DenseDataRange, I4Type<:DenseDataRange}
+    Ca_mech::Vector{T}
+    u_mech::Vector{Vec{Dim, T}}
+    ∇u_mech::∇uType
+    I4_mech::I4Type
+    I4_vecs::Vector{Vector{T}}
+end
+
+struct ElectroMechanicalSynchronizer{NodalTransferOp1, NodalTransferOp2, T, Dim}
+    transfer_op_EP_Mech::NodalTransferOp1
+    transfer_op_Mech_EP::NodalTransferOp2
+    cache::ElectroMechanicalSynchronizerCache{T, Dim}
+end
+
+function ElectroMechanicalSynchronizer(transfer_op_EP_Mech, transfer_op_Mech_EP)
+    qrc=Thunderbolt.QuadratureRuleCollection(1)
+    dh = transfer_op_EP_Mech.dh_to # mechanical dh
+    dim = Ferrite.getrefdim(Ferrite.getfieldinterpolation(dh, (1,1)))
+    Ca_mech = zeros(ndofs(dh))
+    u_mech  = zeros(Vec{dim, Float64}, ndofs(dh))
+    ∇u_mech = Thunderbolt.construct_qvector(Vector{Tensor{2, dim, Float64}}, Vector{Int64}, get_grid(dh), qrc)
+    I4_mech = Thunderbolt.construct_qvector(Vector{Float64}, Vector{Int64}, get_grid(dh), qrc)
+    I4_vecs = [zeros(getnquadpoints(getquadraturerule(qrc, cell))) for cell in getcells(get_grid(dh))] 
+    ElectroMechanicalSynchronizer(
+        transfer_op_EP_Mech,
+        transfer_op_Mech_EP,
+        ElectroMechanicalSynchronizerCache(
+            Ca_mech ,
+            u_mech  ,
+            ∇u_mech ,
+            I4_mech ,
+            I4_vecs
+        )
+    )
+end
+
+function get_calcium_state_idx end
+
+function OS.forward_sync_external!(outer_integrator::OS.OperatorSplittingIntegrator, inner_integrator::DiffEqBase.DEIntegrator, sync::ElectroMechanicalSynchronizer)
+    # Tying holds a buffer for the 3D problem with some meta information about the 0D problem
+    calcium_field_coeff = inner_integrator.cache.inner_solver_cache.op.integrator.volume_model.material_model.contraction_model.calcium_field
+    calcium_state_idx = get_calcium_state_idx(outer_integrator.subintegrator_tree[1][2].f.ode)
+    isnothing(calcium_state_idx) && return nothing # TODO: error?
+    npoints = outer_integrator.subintegrator_tree[1][2].f.npoints
+    calcium_in_EP = @view outer_integrator.subintegrator_tree[1][2].u[npoints * (calcium_state_idx-1)+1:npoints * calcium_state_idx]
+    dh_mech = sync.transfer_op_EP_Mech.dh_to
+    Ca_mech = sync.cache.Ca_mech
+    Thunderbolt.transfer!(Ca_mech, sync.transfer_op_EP_Mech, calcium_in_EP)
+    for sdh in dh_mech.subdofhandlers
+        for cell in CellIterator(sdh)
+            cell_idx = cellid(cell)
+            for dof_idx in 1:length(celldofs(cell))
+                calcium_field_coeff.elementwise_data[dof_idx, cell_idx] = Ca_mech[celldofs(cell)[dof_idx]]
+            end
+        end
+    end
+end
+
+function OS.backward_sync_external!(outer_integrator::OS.OperatorSplittingIntegrator, inner_integrator::DiffEqBase.DEIntegrator, sync::ElectroMechanicalSynchronizer)    # Tying holds a buffer for the 3D problem with some meta information about the 0D problem
+    u_mech_flat = inner_integrator.u
+    u_mech = sync.cache.u_mech
+    dim = Ferrite.getrefdim(Ferrite.getfieldinterpolation(inner_integrator.f.dh, (1,1)))
+    for i in 1:length(u_mech)
+        u_mech[i] = Tensors.Vec(ntuple(j -> u_mech_flat[(i-1)*dim+j], Val(dim)))
+    end
+    dh = inner_integrator.f.dh
+    qrc=Thunderbolt.QuadratureRuleCollection(1)
+    ∇u_mech = sync.cache.∇u_mech
+    ∇u_mech.data .= Ref(zero(eltype(∇u_mech.data)))
+    integrator = Thunderbolt.BilinearDiffusionIntegrator(
+        ConstantCoefficient(diagm(Tensor{2,dim}, 1.0)),
+        qrc, # Allow e.g. mass lumping for explicit integrators.
+        :displacement,
+    )
+    Thunderbolt.compute_quadrature_fluxes!(∇u_mech, sync.transfer_op_Mech_EP.dh_from, u_mech, :φₘ, integrator)
+    fiber_coeff = inner_integrator.f.integrator.volume_model.material_model.microstructure_model.fiber_coefficient
+    ∇u_mech.data .= Ref(diagm(Tensor{2,dim}, 1.0)) .+ ∇u_mech.data
+    I4_mech = sync.cache.I4_mech
+    proj = L2Projector(dh.grid)
+    for sdh in dh.subdofhandlers
+        fiber_coeff_cache = Thunderbolt.setup_coefficient_cache(fiber_coeff, Thunderbolt.getquadraturerule(qrc,sdh), sdh)
+        qr = Thunderbolt.getquadraturerule(qrc,sdh)
+        add!(proj, sdh.cellset, Thunderbolt.getinterpolation(LagrangeCollection{1}(), getcells(dh.grid, first(sdh.cellset))) ; qr_rhs = qr)
+        for cell in CellIterator(sdh)
+            cell_idx = cellid(cell)
+            I4_mech_cell = Thunderbolt.get_data_for_index(I4_mech, cell_idx)
+            F_cell = Thunderbolt.get_data_for_index(∇u_mech, cell_idx)
+            for qp in QuadratureIterator(qr)
+                f = evaluate_coefficient(fiber_coeff_cache, cell, qp, 0.0)
+                I4_mech_cell[qp.i] = (F_cell[qp.i] ⋅ f) ⋅ (F_cell[qp.i] ⋅ f)
+            end
+        end
+    end
+    I4_vecs = sync.cache.I4_vecs
+    cell_idx = 1
+    for sdh in dh.subdofhandlers
+        for cell in CellIterator(sdh)
+            cell_idx = cellid(cell)
+            I4_mech_cell = Thunderbolt.get_data_for_index(I4_mech, cell_idx)
+            I4_vecs[cell_idx] .= I4_mech_cell
+            cell_idx += 1
+        end
+    end
+    close!(proj)
+    projected = project(proj, I4_vecs)
+    I4_in_EP = outer_integrator.subintegrator_tree[1][2].f.ode.I4
+    Thunderbolt.transfer!(I4_in_EP, sync.transfer_op_Mech_EP, (projected))    
+end
+

src/modeling/multiphysics.jl

@@ -1,3 +1,4 @@
 include("coupler/interface.jl")
 include("coupler/tying.jl")
 include("coupler/fsi.jl")
+include("coupler/electromechanics.jl")