diff --git a/Project.toml b/Project.toml index 5ed3a364..24b641da 100644 --- a/Project.toml +++ b/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" diff --git a/docs/Project.toml b/docs/Project.toml index eee86dd6..eeb7bbc2 100644 --- a/docs/Project.toml +++ b/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"} diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl index 817832fa..e29c8541 100644 --- a/src/Thunderbolt.jl +++ b/src/Thunderbolt.jl @@ -223,6 +223,9 @@ export # Coupling Coupling, CoupledModel, + ElectroMechanicalCoupledModel, + ElectroMechanicalCoupler, + ElectroMechanicalSynchronizer, # Discretization semidiscretize, FiniteElementDiscretization, diff --git a/src/discretization/fem.jl b/src/discretization/fem.jl index f5363a47..6d67b2e0 100644 --- a/src/discretization/fem.jl +++ b/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 diff --git a/src/ferrite-addons/transfer_operators.jl b/src/ferrite-addons/transfer_operators.jl index 2e0dacfc..49a8be1b 100644 --- a/src/ferrite-addons/transfer_operators.jl +++ b/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." diff --git a/src/modeling/coupler/electromechanics.jl b/src/modeling/coupler/electromechanics.jl new file mode 100644 index 00000000..b8e100aa --- /dev/null +++ b/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 + diff --git a/src/modeling/multiphysics.jl b/src/modeling/multiphysics.jl index f13893e0..8c64bbee 100644 --- a/src/modeling/multiphysics.jl +++ b/src/modeling/multiphysics.jl @@ -1,3 +1,4 @@ include("coupler/interface.jl") include("coupler/tying.jl") include("coupler/fsi.jl") +include("coupler/electromechanics.jl")