From 8ee2d7fc71a2e7b0c7a75315cfb8d91d29e0ff55 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 22 Jul 2025 13:08:36 +0200 Subject: [PATCH 1/6] Merge branch 'ahm/electromechanics' --- src/Thunderbolt.jl | 3 + src/discretization/fem.jl | 16 ++++ src/ferrite-addons/transfer_operators.jl | 4 +- src/modeling/coupler/electromechanics.jl | 93 ++++++++++++++++++++++++ src/modeling/multiphysics.jl | 1 + 5 files changed, 115 insertions(+), 2 deletions(-) create mode 100644 src/modeling/coupler/electromechanics.jl diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl index 817832fae..e29c85419 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 f5363a474..6d67b2e07 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 2e0dacfcc..49a8be1bf 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 000000000..0c504ed57 --- /dev/null +++ b/src/modeling/coupler/electromechanics.jl @@ -0,0 +1,93 @@ +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 ElectroMechanicalSynchronizer{NodalTransferOp1, NodalTransferOp2} + transfer_op_EP_Mech::NodalTransferOp1 + transfer_op_Mech_EP::NodalTransferOp2 +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 = 4 + npoints = outer_integrator.subintegrator_tree[1][2].f.npoints + calcium_in_EP = 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 = zeros(ndofs(dh_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 = [Tensors.Vec((u_mech_flat[i:i+2]...)) for i in 1:3:length(u_mech_flat)] + dh = inner_integrator.f.dh + qrc=Thunderbolt.QuadratureRuleCollection(1) + dim = Ferrite.getrefdim(Ferrite.getfieldinterpolation(inner_integrator.f.dh, (1,1))) + ∇u_mech = Thunderbolt.construct_qvector(Vector{Tensor{2, dim}}, Vector{Int64}, dh.grid, qrc) + I=SymmetricTensor{2,3}((1.0,0.0,0.0,1.0,0.0,1.0)) + integrator = Thunderbolt.BilinearDiffusionIntegrator( + ConstantCoefficient( I), + # Thunderbolt.NodalQuadratureRuleCollection(Thunderbolt.LagrangeCollection{1}()), # Allow e.g. mass lumping for explicit integrators. + 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 + F = deepcopy(∇u_mech) + F.data .= Ref(I) .+ F.data + I4_mech = Thunderbolt.construct_qvector(Vector{Float64}, Vector{Int64}, dh.grid, qrc) + 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(F, 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_vec_of_vecs = Vector{Float64}[] + 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) + push!(I4_vec_of_vecs, I4_mech_cell) + end + end + close!(proj) + projected = project(proj, I4_vec_of_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 f13893e09..8c64bbee0 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") From afdaa323882d0db9f7861743d0c605e5c326ca2d Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 22 Jul 2025 16:23:05 +0200 Subject: [PATCH 2/6] Caching intermediate arrays --- src/modeling/coupler/electromechanics.jl | 68 ++++++++++++++++++------ 1 file changed, 52 insertions(+), 16 deletions(-) diff --git a/src/modeling/coupler/electromechanics.jl b/src/modeling/coupler/electromechanics.jl index 0c504ed57..d280c406e 100644 --- a/src/modeling/coupler/electromechanics.jl +++ b/src/modeling/coupler/electromechanics.jl @@ -19,19 +19,52 @@ struct ElectroMechanicalCoupledModel{SM #=<: QuasiStaticModel =#, EM, CT <: Elec coupler::CT end -struct ElectroMechanicalSynchronizer{NodalTransferOp1, NodalTransferOp2} +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 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 = 4 npoints = outer_integrator.subintegrator_tree[1][2].f.npoints - calcium_in_EP = outer_integrator.subintegrator_tree[1][2].u[npoints * (calcium_state_idx-1)+1:npoints * calcium_state_idx] + 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 = zeros(ndofs(dh_mech)) + 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) @@ -45,23 +78,24 @@ 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 = [Tensors.Vec((u_mech_flat[i:i+2]...)) for i in 1:3:length(u_mech_flat)] + 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) - dim = Ferrite.getrefdim(Ferrite.getfieldinterpolation(inner_integrator.f.dh, (1,1))) - ∇u_mech = Thunderbolt.construct_qvector(Vector{Tensor{2, dim}}, Vector{Int64}, dh.grid, qrc) - I=SymmetricTensor{2,3}((1.0,0.0,0.0,1.0,0.0,1.0)) + ∇u_mech = sync.cache.∇u_mech + ∇u_mech.data .= Ref(zero(eltype(∇u_mech.data))) integrator = Thunderbolt.BilinearDiffusionIntegrator( - ConstantCoefficient( I), - # Thunderbolt.NodalQuadratureRuleCollection(Thunderbolt.LagrangeCollection{1}()), # Allow e.g. mass lumping for explicit integrators. + 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 - F = deepcopy(∇u_mech) - F.data .= Ref(I) .+ F.data - I4_mech = Thunderbolt.construct_qvector(Vector{Float64}, Vector{Int64}, dh.grid, qrc) + ∇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) @@ -70,23 +104,25 @@ function OS.backward_sync_external!(outer_integrator::OS.OperatorSplittingIntegr 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(F, 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_vec_of_vecs = Vector{Float64}[] + 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) - push!(I4_vec_of_vecs, I4_mech_cell) + I4_vecs[cell_idx] .= I4_mech_cell + cell_idx += 1 end end close!(proj) - projected = project(proj, I4_vec_of_vecs) + 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 From a3acee612a2e2519d330611209a01b230f656097 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 22 Jul 2025 16:23:17 +0200 Subject: [PATCH 3/6] pass CI bitte --- Project.toml | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5ed3a364f..e7a532ee0 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", rev = "7ef3c95"} + [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" From dc4d90d9737c04aa91320ed3af4268a30d9b5fea Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 22 Jul 2025 16:36:12 +0200 Subject: [PATCH 4/6] maybe now the CI is alive --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e7a532ee0..cfd5ccc88 100644 --- a/Project.toml +++ b/Project.toml @@ -39,7 +39,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" [sources] -OrdinaryDiffEqOperatorSplitting = {url = "https://github.com/SciML/OrdinaryDiffEqOperatorSplitting.jl", rev = "7ef3c95"} +OrdinaryDiffEqOperatorSplitting = {url = "https://github.com/SciML/OrdinaryDiffEqOperatorSplitting.jl.git", rev = "7ef3c95"} [extensions] CuThunderboltExt = "CUDA" From 1650730de90811eb0cd462c514b99e2a8c898648 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 22 Jul 2025 16:39:16 +0200 Subject: [PATCH 5/6] once again, deps --- Project.toml | 2 +- docs/Project.toml | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index cfd5ccc88..24b641da1 100644 --- a/Project.toml +++ b/Project.toml @@ -39,7 +39,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" [sources] -OrdinaryDiffEqOperatorSplitting = {url = "https://github.com/SciML/OrdinaryDiffEqOperatorSplitting.jl.git", rev = "7ef3c95"} +OrdinaryDiffEqOperatorSplitting = {url = "https://github.com/SciML/OrdinaryDiffEqOperatorSplitting.jl.git", rev = "main"} [extensions] CuThunderboltExt = "CUDA" diff --git a/docs/Project.toml b/docs/Project.toml index eee86dd63..eeb7bbc22 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"} From 09b050c9f4de1e7383d7675e3a637b9ca58364f4 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 22 Jul 2025 16:55:45 +0200 Subject: [PATCH 6/6] get_calcium_state_idx --- src/modeling/coupler/electromechanics.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/modeling/coupler/electromechanics.jl b/src/modeling/coupler/electromechanics.jl index d280c406e..b8e100aa7 100644 --- a/src/modeling/coupler/electromechanics.jl +++ b/src/modeling/coupler/electromechanics.jl @@ -55,12 +55,13 @@ function ElectroMechanicalSynchronizer(transfer_op_EP_Mech, transfer_op_Mech_EP) ) 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 = 4 + 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