Merge pull request #240 from kaipartmann/newton-raphson-serial

Experimental: Newton Raphson solver (only serial)

Author: kaipartmann

Project.toml

@@ -1,13 +1,14 @@
 name = "Peridynamics"
 uuid = "4dc47793-80f3-4232-b30e-ca78ca9d621b"
-authors = ["Kai Partmann"]
 version = "0.5.0-DEV"
+authors = ["Kai Partmann"]
 
 [deps]
 AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LibGit2 = "76f85450-5226-5b5a-8eaa-529ad045b433"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -26,6 +27,7 @@ AbaqusReader = "0.2"
 Base64 = "1.8"
 CodecZlib = "0.7"
 Dates = "1.8"
+IterativeSolvers = "0.9.4"
 LibGit2 = "1.8"
 LightXML = "0.9"
 LinearAlgebra = "1.8"

docs/src/private_api_reference.md

@@ -8,26 +8,8 @@
 Pages = ["private_api_reference.md"]
 ```
 
+## Types
 ```@docs
-Peridynamics.failure_permit!
-Peridynamics.get_frac_params
-Peridynamics.set_failure_permissions!
-Peridynamics.has_fracture
-Peridynamics.check_pos_and_vol
-Peridynamics.pre_submission_check
-Peridynamics.get_paramsetup
-Peridynamics.get_params
-Peridynamics.check_if_sets_intersect
-Peridynamics.check_if_set_is_defined
-Peridynamics.find_points
-Peridynamics.apply_precracks!
-Peridynamics.apply_precrack!
-Peridynamics.point_sets_intersect
-Peridynamics.displacement_bc!
-Peridynamics.velocity_databc!
-Peridynamics.forcedensity_databc!
-Peridynamics.invreg
-
 Peridynamics.InterfaceError
 Peridynamics.HaloExchange
 Peridynamics.JobOptions
@@ -56,7 +38,37 @@ Peridynamics.BACPointParameters
 Peridynamics.SingleDimIC
 Peridynamics.PosDepSingleDimIC
 Peridynamics.ShortRangeForceContact
+```
 
+## Functions
+```@docs
+Peridynamics.failure_permit!
+Peridynamics.get_frac_params
+Peridynamics.set_failure_permissions!
+Peridynamics.has_fracture
+Peridynamics.check_pos_and_vol
+Peridynamics.pre_submission_check
+Peridynamics.get_paramsetup
+Peridynamics.get_params
+Peridynamics.check_if_sets_intersect
+Peridynamics.check_if_set_is_defined
+Peridynamics.find_points
+Peridynamics.apply_precracks!
+Peridynamics.apply_precrack!
+Peridynamics.point_sets_intersect
+Peridynamics.invreg
+```
+
+## Macros
+```@docs
 Peridynamics.@storage
 Peridynamics.@autoinfiltrate
+```
+
+## Experimental Features
+```@docs
+Peridynamics.displacement_bc!
+Peridynamics.velocity_databc!
+Peridynamics.forcedensity_databc!
+Peridynamics.NewtonRaphson
 ```

src/Peridynamics.jl

@@ -1,7 +1,7 @@
 module Peridynamics
 
 using Base.Threads, Printf, LinearAlgebra, StaticArrays, PointNeighbors, ProgressMeter,
-      WriteVTK, TimerOutputs, MPI, PrecompileTools
+      WriteVTK, TimerOutputs, MPI, PrecompileTools, IterativeSolvers
 @static if Sys.islinux()
     using ThreadPinning
 end
@@ -135,6 +135,7 @@ include("core/mpi_multibody_data_handler.jl")
 
 include("time_solvers/velocity_verlet.jl")
 include("time_solvers/dynamic_relaxation.jl")
+include("time_solvers/newton_raphson.jl")
 
 include("physics/bond_based.jl")
 include("physics/dh_bond_based.jl")

src/auxiliary/docs.jl

@@ -11,3 +11,14 @@ function internal_api_warning()
     end
     return msg
 end
+
+function experimental_api_warning()
+    msg = """
+    !!! danger "Experimental feature"
+        Please note that this is an experimental feature. It is *not*
+        part of the public API of Peridynamics.jl, and thus can be altered (or removed)
+        at any time without it being considered a breaking change. Also, the feature may be
+        incomplete and/or contain bugs. Please use with caution.
+    """
+    return msg
+end

src/core/time_solvers.jl

@@ -29,30 +29,6 @@ function required_fields_timesolvers()
     return Tuple(unique(all_fields))
 end
 
-# function required_point_data_fields_timesolvers()
-#     all_fields = Vector{Symbol}()
-#     for solver in registered_solvers()
-#         push!(all_fields, req_point_data_fields_timesolver(solver)...)
-#     end
-#     return Tuple(unique(all_fields))
-# end
-
-# function required_bond_data_fields_timesolvers()
-#     all_fields = Vector{Symbol}()
-#     for solver in registered_solvers()
-#         push!(all_fields, req_bond_data_fields_timesolver(solver)...)
-#     end
-#     return Tuple(unique(all_fields))
-# end
-
-# function required_data_fields_timesolvers()
-#     all_fields = Vector{Symbol}()
-#     for solver in registered_solvers()
-#         push!(all_fields, req_data_fields_timesolver(solver)...)
-#     end
-#     return Tuple(unique(all_fields))
-# end
-
 function req_point_data_fields_timesolver(::Type{TS}) where {TS}
     return throw(InterfaceError(TS, "req_point_data_fields_timesolver"))
 end

src/discretization/bond_system.jl

@@ -261,6 +261,22 @@ function calc_force_density!(chunk::AbstractBodyChunk{<:AbstractBondSystem}, t,
     return nothing
 end
 
+function calc_force_density!(storage::AbstractStorage, system::AbstractBondSystem,
+                             mat::AbstractBondSystemMaterial,
+                             paramsetup::AbstractParameterSetup, idxs::AbstractVector{Int},
+                             t, Δt)
+    (; dmgmodel) = mat
+    @inbounds storage.b_int[:, idxs] .= 0.0
+    @inbounds storage.n_active_bonds[idxs] .= 0
+    for i in idxs
+        calc_failure!(storage, system, mat, dmgmodel, paramsetup, i)
+        calc_damage!(storage, system, mat, dmgmodel, paramsetup, i)
+        force_density_point!(storage, system, mat, paramsetup, t, Δt, i)
+    end
+    nancheck(storage, t, Δt)
+    return nothing
+end
+
 function calc_force_density!(storage::AbstractStorage, system::AbstractBondSystem,
                              mat::AbstractBondSystemMaterial,
                              paramsetup::AbstractParameterSetup, t, Δt)
@@ -322,6 +338,13 @@ function calc_damage!(storage::AbstractStorage, system::AbstractBondSystem,
     return nothing
 end
 
+function get_affected_points(system::AbstractBondSystem, i)
+    points = [bond.neighbor for bond in system.bonds[each_bond_idx(system, i)]]
+    pushfirst!(points, i)
+    sort!(points)
+    return points
+end
+
 function log_system(::Type{System}, options::AbstractJobOptions,
                     dh::AbstractDataHandler) where {System<:AbstractBondSystem}
     n_bonds = calc_n_bonds(dh)

src/discretization/interaction_system.jl

@@ -428,6 +428,22 @@ function calc_force_density!(chunk::AbstractBodyChunk{<:InteractionSystem}, t, 
     return nothing
 end
 
+function calc_force_density!(storage::AbstractStorage, system::InteractionSystem,
+                             mat::AbstractInteractionSystemMaterial,
+                             paramsetup::AbstractParameterSetup, idxs::AbstractVector{Int},
+                             t, Δt)
+    (; dmgmodel) = mat
+    @inbounds storage.b_int[:, idxs] .= 0
+    @inbounds storage.n_active_one_nis[idxs] .= 0
+    for i in idxs
+        calc_failure!(storage, system, mat, dmgmodel, paramsetup, i)
+        calc_damage!(storage, system, mat, dmgmodel, paramsetup, i)
+        force_density_point!(storage, system, mat, paramsetup, t, Δt, i)
+    end
+    nancheck(storage, t, Δt)
+    return nothing
+end
+
 function calc_failure!(storage::AbstractStorage, system::InteractionSystem,
                        mat::AbstractInteractionSystemMaterial, dmgmodel::CriticalStretch,
                        paramsetup::AbstractParameterSetup, i)
@@ -467,6 +483,13 @@ end
     return storage.one_ni_active[one_ni_id]
 end
 
+function get_affected_points(system::InteractionSystem, i)
+    points = [one_ni.neighbor for one_ni in system.one_nis[each_one_ni_idx(system, i)]]
+    pushfirst!(points, i)
+    sort!(points)
+    return points
+end
+
 function log_msg_interaction_system(n_one_nis::Int, n_two_nis::Int, n_three_nis::Int)
     msg = msg_qty("number of one-neighbor-interactions", n_one_nis)
     msg *= msg_qty("number of two-neighbor-interactions", n_two_nis)

src/physics/ba_correspondence.jl

@@ -168,10 +168,18 @@ end
     @pointfield b_ext::Matrix{Float64}
     @pointfield density_matrix::Matrix{Float64}
     @pointfield damage::Vector{Float64}
-    bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
     @pointfield stress::Matrix{Float64}
     @pointfield von_mises_stress::Vector{Float64}
+    bond_active::Vector{Bool}
+    residual::Vector{Float64}
+    jacobian::Matrix{Float64}
+    displacement_copy::Matrix{Float64}
+    b_int_copy::Matrix{Float64}
+    temp_force_a::Vector{Float64}
+    temp_force_b::Vector{Float64}
+    Δu::Vector{Float64}
+    affected_points::Vector{Vector{Int}}
 end
 
 function init_field(::BACMaterial, ::AbstractTimeSolver, system::BondAssociatedSystem,

src/physics/bond_based.jl

@@ -126,10 +126,18 @@ end
     @pointfield b_ext::Matrix{Float64}
     @pointfield density_matrix::Matrix{Float64}
     @pointfield damage::Vector{Float64}
-    bond_length::Vector{Float64}
-    bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
     @pointfield strain_energy_density::Vector{Float64}
+    bond_length::Vector{Float64}
+    bond_active::Vector{Bool}
+    residual::Vector{Float64}
+    jacobian::Matrix{Float64}
+    displacement_copy::Matrix{Float64}
+    b_int_copy::Matrix{Float64}
+    temp_force_a::Vector{Float64}
+    temp_force_b::Vector{Float64}
+    Δu::Vector{Float64}
+    affected_points::Vector{Vector{Int}}
 end
 
 function init_field(::AbstractBondBasedMaterial, ::AbstractTimeSolver, system::BondSystem,

src/physics/boundary_conditions.jl

@@ -438,6 +438,7 @@ end
     velocity_databc!(body, data, set_name, dims)
 
 $(internal_api_warning())
+$(experimental_api_warning())
 
 Specifies velocity boundary conditions for points of the set `set_name` in `body`.
 The value of the boundary condition is assigned by reading the corresponding positions
@@ -577,6 +578,7 @@ end
     forcedensity_databc!(body, data, set_name, dims)
 
 $(internal_api_warning())
+$(experimental_api_warning())
 
 Specifies forcedensity boundary conditions for points of the set `set_name` in `body`.
 The value of the boundary condition is assigned by reading the corresponding positions
@@ -633,12 +635,12 @@ end
     displacement_bc!(f::Function, body::Body, points::Vector{Int}, dim::Int)
 
 $(internal_api_warning())
+$(experimental_api_warning())
 
-Apply a displacement boundary condition to specified points in a given dimension.
-
-!!! warn "Compatibility limited"
-    This boundary condition type only works with static solvers for now.
+!!! warning "Compatibility limited"
+    This boundary condition type only works with the NewtonRaphson solver for now.
 
+Apply a displacement boundary condition to specified points in a given dimension.
 A factor `β` is calculated as `n / n_steps`, where `n` is the current step number and
 `n_steps` is the total number of steps. Then the return value of the function `f` is
 multiplicated with `β` at each time step. This means that the displacement will be applied