Reduce OS module allocations

src/Thunderbolt.jl

@@ -206,7 +206,6 @@ export
     AdaptiveForwardEulerSubstepper,
     # Integrator
     get_parent_index,
-    get_parent_value,
     # Utils
     calculate_volume_deformed_mesh,
     elementtypes,

src/modeling/core/coefficients.jl

@@ -115,13 +115,15 @@ struct CoordinateSystemCoefficient{CS}
 end
 
 function compute_nodal_values(csc::CoordinateSystemCoefficient, dh::DofHandler, field_name::Symbol)
-    nodal_values = Vector{value_type(csc.cs)}(UndefInitializer(), ndofs(dh))
+    Tv = value_type(csc.cs)
+    nodal_values = Vector{Tv}(UndefInitializer(), ndofs(dh))
+    T = eltype(Tv)
     for sdh in dh.subdofhandlers
         field_name ∈ sdh.field_names || continue
-        ip = Ferrite.getfieldinterpolation(sdh, field_name)
-        positions = Ferrite.reference_coordinates(ip)
+        ip   = Ferrite.getfieldinterpolation(sdh, field_name)
+        rdim = Ferrite.getrefdim(ip)
+        positions = Vec{rdim,T}.(Ferrite.reference_coordinates(ip))
         # This little trick uses the delta property of interpolations
-        T = eltype(first(positions))
         qr = QuadratureRule{Ferrite.getrefshape(ip)}([T(1.0) for _ in 1:length(positions)], positions)
         cc = setup_coefficient_cache(csc, qr, sdh)
         _compute_nodal_values!(nodal_values, qr, cc, sdh)

src/modeling/core/coordinate_systems.jl

@@ -3,10 +3,13 @@
 
 Standard cartesian coordinate system.
 """
-struct CartesianCoordinateSystem{sdim}
+struct CartesianCoordinateSystem{sdim,T}
+    function CartesianCoordinateSystem{sdim}() where sdim
+        return new{sdim,Float32}()
+    end
 end
 
-value_type(::CartesianCoordinateSystem{sdim}) where sdim = Vec{sdim, Float32}
+value_type(::CartesianCoordinateSystem{sdim, T}) where {sdim, T} = Vec{sdim, T}
 
 CartesianCoordinateSystem(mesh::AbstractGrid{sdim}) where sdim = CartesianCoordinateSystem{sdim}()
 
@@ -43,11 +46,13 @@ LV only part of the universal ventricular coordinate, containing
 """
 struct LVCoordinate{T}
     transmural::T
-    apicaobasal::T
+    apicobasal::T
     rotational::T
 end
 
-value_type(::LVCoordinateSystem) = LVCoordinate{Float32}
+Base.eltype(::Type{LVCoordinate{T}}) where T = T
+Base.eltype(::LVCoordinate{T}) where T = T
+value_type(::LVCoordinateSystem{T}) where T = LVCoordinate{T}
 
 
 """
@@ -68,8 +73,8 @@ Requires a mesh with facetsets
 and a nodeset
     * Apex
 """
-function compute_lv_coordinate_system(mesh::SimpleMesh{3,<:Any,T}, subdomains::Vector{String} = [""]; up = Vec((T(0.0),T(0.0),T(1.0)))) where T
-    @assert up ≈ Vec((T(0.0),T(0.0),T(1.0))) "Custom up vector not yet supported."
+function compute_lv_coordinate_system(mesh::SimpleMesh{3,<:Any,T}, subdomains::Vector{String} = [""]; up = Vec((T(0.0),T(0.0),T(-1.0)))) where T
+    @assert abs.(up) ≈ Vec((T(0.0),T(0.0),T(1.0))) "Custom up vector not yet supported."
     ip_collection = LagrangeCollection{1}()
     qr_collection = QuadratureRuleCollection(2)
     cv_collection = CellValueCollection(qr_collection, ip_collection)
@@ -83,7 +88,6 @@ function compute_lv_coordinate_system(mesh::SimpleMesh{3,<:Any,T}, subdomains::V
     # Assemble Laplacian
     # TODO use bilinear operator for performance
     K = allocate_matrix(dh)
-
     assembler = start_assemble(K)
     for sdh in dh.subdofhandlers
         cellvalues = getcellvalues(cv_collection, getcells(mesh, first(sdh.cellset)))
@@ -119,38 +123,19 @@ function compute_lv_coordinate_system(mesh::SimpleMesh{3,<:Any,T}, subdomains::V
     close!(ch)
     update!(ch, 0.0);
 
-    K_transmural = copy(K)
+    K_transmural = K
     f = zeros(ndofs(dh))
 
     apply!(K_transmural, f, ch)
-    transmural = K_transmural \ f;
+    sol = solve(LinearSolve.LinearProblem(K_transmural, f), LinearSolve.KrylovJL_CG())
+    transmural = sol.u
 
     # Apicobasal coordinate
-    #TODO refactor check for node set existence
-    if !haskey(mesh.grid.nodesets, "Apex") #TODO this is just a hotfix, assuming that z points towards the apex
-        apex_node_index = 1
-        nodes = getnodes(mesh)
-        for (i,node) ∈ enumerate(nodes)
-            if nodes[i].x[3] > nodes[apex_node_index].x[3]
-                apex_node_index = i
-            end
-        end
-        addnodeset!(mesh, "Apex", OrderedSet{Int}((apex_node_index)))
-    end
-
-    ch = ConstraintHandler(dh);
-    dbc = Dirichlet(:coordinates, getfacetset(mesh, "Base"), (x, t) -> 0)
-    Ferrite.add!(ch, dbc);
-    dbc = Dirichlet(:coordinates, getnodeset(mesh, "Apex"), (x, t) -> 1)
-    Ferrite.add!(ch, dbc);
-    close!(ch)
-    update!(ch, 0.0);
-
-    K_apicobasal = copy(K)
-    f = zeros(ndofs(dh))
-
-    apply!(K_apicobasal, f, ch)
-    apicobasal = K_apicobasal \ f;
+    apicobasal = zeros(ndofs(dh))
+    apply_analytical!(apicobasal, dh, :coordinates, x->x ⋅ up)
+    apicobasal .-= minimum(apicobasal)
+    apicobasal = abs.(apicobasal)
+    apicobasal ./= maximum(apicobasal)
 
     rotational = zeros(ndofs(dh))
     rotational .= NaN
@@ -175,7 +160,7 @@ function compute_lv_coordinate_system(mesh::SimpleMesh{3,<:Any,T}, subdomains::V
                     rotational[dofs[qp.i]] = 0.0
                 else
                     x = x_planar / xlen
-                    rotational[dofs[qp.i]] = (π + atan(x[1], x[2]))/2 # TODO tilted coordinate system
+                    rotational[dofs[qp.i]] = 1/2 + atan(x[1], x[2])/(2π) # TODO tilted coordinate system
                 end
             end
         end
@@ -194,7 +179,7 @@ Requires a mesh with facetsets
     * Myocardium
 """
 function compute_midmyocardial_section_coordinate_system(mesh::SimpleMesh{3,<:Any,T}, subdomains::Vector{String} = [""]; up = Vec((T(0.0),T(0.0),T(1.0)))) where T
-    @assert up ≈ Vec((T(0.0),T(0.0),T(1.0))) "Custom up vector not yet supported."
+    @assert abs.(up) ≈ Vec((T(0.0),T(0.0),T(1.0))) "Custom up vector not yet supported."
     ip_collection = LagrangeCollection{1}()
     qr_collection = QuadratureRuleCollection(2)
     cv_collection = CellValueCollection(qr_collection, ip_collection)
@@ -244,26 +229,22 @@ function compute_midmyocardial_section_coordinate_system(mesh::SimpleMesh{3,<:An
     close!(ch)
     update!(ch, 0.0);
 
-    K_transmural = copy(K)
+    K_transmural = K
     f = zeros(ndofs(dh))
 
     apply!(K_transmural, f, ch)
-    transmural = K_transmural \ f;
-
-    ch = ConstraintHandler(dh);
-    dbc = Dirichlet(:coordinates, getfacetset(mesh, "Base"), (x, t) -> 0)
-    Ferrite.add!(ch, dbc);
-    dbc = Dirichlet(:coordinates, getfacetset(mesh, "Myocardium"), (x, t) -> 0.15)
-    Ferrite.add!(ch, dbc);
-    close!(ch)
-    update!(ch, 0.0);
-
-    K_apicobasal = copy(K)
-    f = zeros(ndofs(dh))
-
-    apply!(K_apicobasal, f, ch)
-    apicobasal = K_apicobasal \ f;
+    sol = solve(LinearSolve.LinearProblem(K_transmural, f), LinearSolve.KrylovJL_CG())
+    transmural = sol.u
 
+    # Apicobasal coordinate
+    apicobasal = zeros(ndofs(dh))
+    apply_analytical!(apicobasal, dh, :coordinates, x->x ⋅ up)
+    apicobasal .-= minimum(apicobasal)
+    apicobasal = abs.(apicobasal)
+    apicobasal ./= maximum(apicobasal)
+    apicobasal .*= 0.15
+
+    # Rotational coordinate
     rotational = zeros(ndofs(dh))
     rotational .= NaN
 
@@ -283,7 +264,7 @@ function compute_midmyocardial_section_coordinate_system(mesh::SimpleMesh{3,<:An
                 x_planar = x_dof - (x_dof ⋅ up) * up # Project into plane
                 x = x_planar / norm(x_planar)
 
-                rotational[dofs[qp.i]] = (π + atan(x[1], x[2]))/2 # TODO tilted coordinate system
+                rotational[dofs[qp.i]] = 1/2 + atan(x[1], x[2])/(2π) # TODO tilted coordinate system
             end
         end
     end
@@ -328,11 +309,13 @@ Biventricular universal coordinate, containing
 """
 struct BiVCoordinate{T}
     transmural::T
-    apicaobasal::T
+    apicobasal::T
     rotational::T
     transventricular::T
 end
 
+Base.eltype(::Type{BiVCoordinate{T}}) where T = T
+Base.eltype(::BiVCoordinate{T}) where T = T
 value_type(::BiVCoordinateSystem) = BiVCoordinate
 
 getcoordinateinterpolation(cs::BiVCoordinateSystem, cell::Ferrite.AbstractCell) = Ferrite.getfieldinterpolation(cs.dh, (1,1))

src/solver/operator_splitting/integrator.jl

@@ -72,18 +72,15 @@ function DiffEqBase.__init(
 
     callback = DiffEqBase.CallbackSet(callback)
 
-    cache = init_cache(prob, alg; dt, kwargs...)
-
-    u = cache.u
-    uprev = cache.uprev
+    cache = init_cache(prob, alg; u0, t0, dt, kwargs...)
 
-    subintegrators = build_subintegrators_recursive(prob.f, prob.f.synchronizers, p, cache, u, uprev, t0, dt, 1:length(u0), u, tstops, _tstops, saveat, _saveat)
+    subintegrators = build_subintegrators_recursive(prob.f, prob.f.synchronizers, p, cache, t0, dt, 1:length(u0), cache.u, tstops, _tstops, saveat, _saveat)
 
     integrator = OperatorSplittingIntegrator(
         prob.f,
         alg,
-        u,
-        uprev,
+        cache.u,
+        cache.uprev,
         p,
         t0,
         copy(t0),
@@ -285,8 +282,8 @@ function __step!(integrator)
     synchronize_subintegrators!(integrator)
     tnext = integrator.t + integrator.dt
 
-     # Solve inner problems
-    advance_solution_to!(integrator, tnext)
+    # Solve inner problems
+    advance_solution_to!(integrator, tnext; uparent=integrator.u)
     stepsize_controller!(integrator)
 
     # Update integrator
@@ -307,8 +304,8 @@ function __step!(integrator)
 end
 
 # solvers need to define this interface
-function advance_solution_to!(integrator, tnext)
-    advance_solution_to!(integrator, integrator.cache, tnext)
+function advance_solution_to!(integrator, tnext; uparent)
+    advance_solution_to!(integrator, integrator.cache, tnext; uparent)
 end
 
 DiffEqBase.get_dt(integrator::OperatorSplittingIntegrator) = integrator._dt
@@ -344,21 +341,19 @@ end
     end
 end
 
-advance_solution_to!(integrator::OperatorSplittingIntegrator, cache::AbstractOperatorSplittingCache, tnext::Number) = advance_solution_to!(integrator.subintegrators, cache, tnext)
+function advance_solution_to!(integrator::OperatorSplittingIntegrator, cache::AbstractOperatorSplittingCache, tnext::Number; uparent)
+    advance_solution_to!(integrator.subintegrators, cache, tnext; uparent)
+end
 
 # Dispatch for tree node construction
-function build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::Tuple, p::Tuple, cache::AbstractOperatorSplittingCache, u::AbstractArray, uprev::AbstractArray, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat)
+function build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::Tuple, p::Tuple, cache::AbstractOperatorSplittingCache, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat)
     return ntuple(i ->
         build_subintegrators_recursive(
             get_operator(f, i),
             synchronizers[i],
             p[i],
             cache.inner_caches[i],
             # TODO recover this
-            # cache.inner_caches[i].u,
-            # cache.inner_caches[i].uprev,
-            similar(u, length(f.dof_ranges[i])),
-            similar(uprev, length(f.dof_ranges[i])),
             t, dt, f.dof_ranges[i],
             # We pass the full solution, because some parameters might require
             # access to solution variables which are not part of the local solution range
@@ -367,18 +362,14 @@ function build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::
         ), length(f.functions)
     )
 end
-function build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::NoExternalSynchronization, p::Tuple, cache::AbstractOperatorSplittingCache, u::AbstractArray, uprev::AbstractArray, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat)
+function build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::NoExternalSynchronization, p::Tuple, cache::AbstractOperatorSplittingCache, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat)
     return ntuple(i ->
         build_subintegrators_recursive(
             get_operator(f, i),
             synchronizers,
             p[i],
             cache.inner_caches[i],
             # TODO recover this
-            # cache.inner_caches[i].u,
-            # cache.inner_caches[i].uprev,
-            similar(u, length(f.dof_ranges[i])),
-            similar(uprev, length(f.dof_ranges[i])),
             t, dt, f.dof_ranges[i],
             # We pass the full solution, because some parameters might require
             # access to solution variables which are not part of the local solution range
@@ -388,14 +379,14 @@ function build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::
     )
 end
 
-@unroll function prepare_local_step!(subintegrators::Tuple)
+@unroll function prepare_local_step!(uparent, subintegrators::Tuple)
     @unroll for subintegrator in subintegrators
-        prepare_local_step!(subintegrator)
+        prepare_local_step!(uparent, subintegrator)
     end
 end
 
-@unroll function finalize_local_step!(subintegrators::Tuple)
+@unroll function finalize_local_step!(uparent, subintegrators::Tuple)
     @unroll for subintegrator in subintegrators
-        finalize_local_step!(subintegrator)
+        finalize_local_step!(uparent, subintegrator)
     end
 end

src/solver/operator_splitting/solver.jl

@@ -13,36 +13,33 @@ end
 struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplittingCache
     u::uType
     uprev::uType # True previous solution
-    uprev2::tmpType # Previous solution used during time marching
     tmp::tmpType # Scratch
     inner_caches::iiType
 end
 
 # Dispatch for outer construction
-function init_cache(prob::OperatorSplittingProblem, alg::LieTrotterGodunov; dt, kwargs...) # TODO
+function init_cache(prob::OperatorSplittingProblem, alg::LieTrotterGodunov; u0, kwargs...) # TODO
     @unpack f = prob
     @assert f isa GenericSplitFunction
 
-    u          = copy(prob.u0)
-    uprev      = copy(prob.u0)
-
     # Build inner integrator
-    return construct_inner_cache(f, alg, u, uprev)
+    return construct_inner_cache(f, alg; uparent=u0, u0, kwargs...)
 end
 
 # Dispatch for recursive construction
-function construct_inner_cache(f::AbstractOperatorSplitFunction, alg::LieTrotterGodunov, u::AbstractArray, uprev::AbstractArray)
+function construct_inner_cache(f::AbstractOperatorSplitFunction, alg::LieTrotterGodunov; uparent, u0, kwargs...)
     dof_ranges = f.dof_ranges
 
-    uprev2     = similar(uprev)
+    u          = copy(u0)
+    uprev      = copy(u0)
     tmp        = similar(u)
-    inner_caches = ntuple(i->construct_inner_cache(get_operator(f, i), alg.inner_algs[i], similar(u, length(dof_ranges[i])), similar(u, length(dof_ranges[i]))), length(f.functions))
-    LieTrotterGodunovCache(u, uprev, uprev2, tmp, inner_caches)
+    inner_caches = ntuple(i->construct_inner_cache(get_operator(f, i), alg.inner_algs[i]; uparent, u0=view(uparent,dof_ranges[i]), kwargs...), length(f.functions))
+    LieTrotterGodunovCache(u, uprev, tmp, inner_caches)
 end
 
-@inline @unroll function advance_solution_to!(subintegrators::Tuple, cache::LieTrotterGodunovCache, tnext)
+@inline @unroll function advance_solution_to!(subintegrators::Tuple, cache::LieTrotterGodunovCache, tnext; uparent)
     # We assume that the integrators are already synced
-    @unpack u, uprev2, uprev, inner_caches = cache
+    @unpack u, uprev, inner_caches = cache
 
     # Store current solution
     uprev .= u
@@ -51,8 +48,8 @@ end
     i = 0
     @unroll for subinteg in subintegrators
         i += 1
-        prepare_local_step!(subinteg)
-        advance_solution_to!(subinteg, inner_caches[i], tnext)
-        finalize_local_step!(subinteg)
+        prepare_local_step!(uparent, subinteg)
+        advance_solution_to!(subinteg, inner_caches[i], tnext; uparent)
+        finalize_local_step!(uparent, subinteg)
     end
 end 

src/solver/time/euler.jl

@@ -86,7 +86,7 @@ function setup_solver_cache(f::TransientDiffusionFunction, solver::BackwardEuler
     uprev = create_system_vector(solver.solution_vector_type, f)
     tmp   = create_system_vector(solver.solution_vector_type, f)
 
-    T = eltype(A)
+    T = eltype(u0)
 
     qr = create_quadrature_rule(f, solver, field_name)