Visco time step [Second part]

src/ComputationalModels/PostProcessors.jl

@@ -92,15 +92,13 @@ function Cauchy(args...)
 end
 
 
-function Piola(model::Elasto, km::KinematicModel, uh, unh, state_vars, Ω, dΩ, t, Δt=0)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
+function Piola(model::Elasto, km::KinematicModel, uh, unh, state_vars, Ω, dΩ, t=0)
   σh = Piola(model,km,uh)
   interpolate_L2_tensor(σh, Ω, dΩ)
 end
 
 
-function Piola(model::ViscoElastic, km::KinematicModel, uh, unh, state_vars, Ω, dΩ, t=0, Δt=0)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
+function Piola(model::ViscoElastic, km::KinematicModel, uh, unh, state_vars, Ω, dΩ, t=0)
   σh = Piola(model, km, uh, unh, state_vars)
   interpolate_L2_tensor(σh, Ω, dΩ)
 end
@@ -113,9 +111,8 @@ function Piola(model::Elasto, km::KinematicModel, uh, vars...)
 end
 
 
-function Piola(model::ViscoElastic,  km::KinematicModel, uh, unh, states, Δt=model.Δt[])
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  _, ∂Ψu, _ = model(Δt=Δt)
+function Piola(model::ViscoElastic,  km::KinematicModel, uh, unh, states)
+  _, ∂Ψu, _ = model()
   F, _, _ = get_Kinematics(km)
   ∂Ψu ∘ (F∘∇(uh), F∘∇(unh), states...)
 end

src/Exports.jl

@@ -89,6 +89,8 @@ end
 @publish PhysicalModels Dissipation
 @publish PhysicalModels initializeStateVariables
 @publish PhysicalModels updateStateVariables!
+@publish PhysicalModels initialize_state
+@publish PhysicalModels update_time_step!
 
 @publish WeakForms residual
 @publish WeakForms jacobian

src/PhysicalModels/ElectroMechanicalModels.jl

@@ -33,8 +33,8 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
     ∂Ψφφ(F, E, N) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
-  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic})(Λ::Float64=1.0; Δt)
-    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ, Δt=Δt)
+  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic})(Λ::Float64=1.0)
+    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano()
     Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano, Λ)
     Ψ(F, E, Fn, A...) = Ψm(F, Fn, A...) + Ψem(F, E)
     ∂Ψu(F, E, Fn, A...) = ∂Ψm_u(F, Fn, A...) + ∂Ψem_u(F, E)
@@ -44,11 +44,23 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
     ∂Ψφφ(F, E, Fn, A...) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
-
 end
 
 ViscoElectricModel = ElectroMechModel{<:Electro,<:ViscoElastic}
 
+function update_time_step!(obj::ElectroMechModel, Δt::Float64)
+  update_time_step!(obj.electro, Δt)
+  update_time_step!(obj.mechano, Δt)
+end
+
+function initialize_state(obj::ElectroMechano, points::Measure)
+  initialize_state(obj.mechano, points)
+end
+
+function update_state!(obj::ElectroMechano, args...)
+  update_state!(obj.mechano, args...)
+end
+
 function initializeStateVariables(obj::ElectroMechano, points::Measure)
   initializeStateVariables(obj.mechano, points)
 end

src/PhysicalModels/PhysicalModels.jl

@@ -13,6 +13,7 @@ using ..TensorAlgebra: trAA
 using StaticArrays
 
 import Base: +
+import Gridap: update_state!
 
 export Yeoh3D
 export NeoHookean3D
@@ -74,13 +75,14 @@ export ThermoElectro
 export FlexoElectro
 export EnergyInterpolationScheme
 export SecondPiola
+export Dissipation
 
 export DerivativeStrategy
 
 export initializeStateVariables
 export updateStateVariables!
-export update_state!
-export set_time_step!
+export initialize_state
+export update_time_step!
 
 export Kinematics
 export KinematicDescription
@@ -146,14 +148,19 @@ include("PINNs.jl")
 """
 Initialize the state variables for the given constitutive model and discretization.
 """
+function initialize_state(::PhysicalModel, points::Measure)
+  return nothing
+end
 function initializeStateVariables(::PhysicalModel, points::Measure)
   return nothing
 end
 
 
 """
-Update the state variables. The state variables must be initialized using the function 'initializeStateVariables'.
+Update the state variables. The state variables must be initialized using the function 'initialize_state'.
 """
+function update_state!(::PhysicalModel, vars...)
+end
 function updateStateVariables!(::Any, ::PhysicalModel, vars...)
 end
 
@@ -177,7 +184,7 @@ end
 """
 Set the time step to be used internally by the constitutive model.
 """
-function set_time_step!(::PhysicalModel, Δt::Float64)
+function update_time_step!(::PhysicalModel, Δt::Float64)
   Δt
 end
 

src/PhysicalModels/ThermoElectroMechanicalModels.jl

@@ -153,9 +153,15 @@ struct ThermoElectroMech_Bonet{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectro
   end
 end
 
-function Dissipation(obj::ThermoElectroMech_Bonet, Δt)
+function update_time_step!(obj::ThermoElectroMech_Bonet, Δt::Float64)
+  update_time_step!(obj.thermo, Δt)
+  update_time_step!(obj.electro, Δt)
+  update_time_step!(obj.mechano, Δt)
+end
+
+function Dissipation(obj::ThermoElectroMech_Bonet)
   @unpack Cv,θr, α, κ, γv, γd = obj.thermo
-  Dvis = Dissipation(obj.mechano, Δt)
+  Dvis = Dissipation(obj.mechano)
   gd(δθ) = 1/(γd+1) * (((δθ+θr)/θr)^(γd+1) -1)
   ∂gd(δθ) = (δθ+θr)^γd / θr^(γd+1)
   D(F, E, δθ, X...) = (1 + gd(δθ)) * Dvis(F, X...)

src/PhysicalModels/ViscousModels.jl

@@ -13,9 +13,7 @@ struct ViscousIncompressible <: Visco
   function ViscousIncompressible(elasto; τ::Float64)
     new(elasto, τ, 0)
   end
-  function (obj::ViscousIncompressible)(Λ::Float64=1.0; Δt::Float64)
-    @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-    obj.Δt[] = Δt
+  function (obj::ViscousIncompressible)()
     Ψe, Se, ∂Se∂Ce   = SecondPiola(obj.elasto)
     Ψ(F, Fn, A)      = Energy(obj, Ψe, Se, ∂Se∂Ce, F, Fn, A)
     ∂Ψ∂F(F, Fn, A)   = Piola(obj, Se, ∂Se∂Ce, F, Fn, A)
@@ -24,22 +22,33 @@ struct ViscousIncompressible <: Visco
   end
 end
 
+function update_time_step!(obj::ViscousIncompressible, Δt::Float64)
+  obj.Δt[] = Δt
+end
+
+function initialize_state(::ViscousIncompressible, points::Measure)
+  v = VectorValue(1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0,0.0)
+  CellState(v, points)
+end
+
+function update_state!(obj::ViscousIncompressible, state, F, Fn)
+  _, Se, ∂Se∂Ce = SecondPiola(obj.elasto)
+  return_mapping(A, F, Fn) = ReturnMapping(obj, Se, ∂Se∂Ce, F, Fn, A)
+  update_state!(return_mapping, state, F, Fn)
+end
+
 function initializeStateVariables(::ViscousIncompressible, points::Measure)
   v = VectorValue(1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0,0.0)
   CellState(v, points)
 end
 
-function updateStateVariables!(state, obj::ViscousIncompressible, Δt, F, Fn)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  obj.Δt[] = Δt
+function updateStateVariables!(state, obj::ViscousIncompressible, F, Fn)
   _, Se, ∂Se∂Ce = SecondPiola(obj.elasto)
   return_mapping(A, F, Fn) = ReturnMapping(obj, Se, ∂Se∂Ce, F, Fn, A)
   update_state!(return_mapping, state, F, Fn)
 end
 
-function Dissipation(obj::ViscousIncompressible, Δt)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  obj.Δt[] = Δt
+function Dissipation(obj::ViscousIncompressible)
   _, Se, ∂Se∂Ce = SecondPiola(obj.elasto)
   D(F, Fn, A) = ViscousDissipation(obj, Se, ∂Se∂Ce, F, Fn, A)
 end
@@ -51,11 +60,9 @@ struct GeneralizedMaxwell <: ViscoElastic
   function GeneralizedMaxwell(longTerm::Elasto, branches::Visco...)
     new(longTerm,branches,0)
   end
-  function (obj::GeneralizedMaxwell)(Λ::Float64=1.0; Δt::Float64)
-    @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-    obj.Δt[] = Δt
+  function (obj::GeneralizedMaxwell)()
     Ψe, ∂Ψeu, ∂Ψeuu = obj.longterm()
-    DΨv = map(b -> b(Δt=Δt), obj.branches)
+    DΨv = map(b -> b(), obj.branches)
     Ψα, ∂Ψαu, ∂Ψαuu = map(i -> getindex.(DΨv, i), 1:3)
     Ψα, ∂Ψαu, ∂Ψαuu = transpose(DΨv)
     Ψ(F, Fn, A...) = mapreduce((Ψi, Ai) -> Ψi(F, Fn, Ai), +, Ψα, A; init=Ψe(F))
@@ -65,29 +72,33 @@ struct GeneralizedMaxwell <: ViscoElastic
   end
 end
 
-function set_time_step!(obj::GeneralizedMaxwell, Δt::Float64)
+function update_time_step!(obj::GeneralizedMaxwell, Δt::Float64)
+  foreach(b -> update_time_step!(b, Δt), obj.branches)
   obj.Δt[] = Δt
-  foreach(b -> b.Δt[] = Δt, obj.branches)
-  Δt
+end
+
+function initialize_state(obj::GeneralizedMaxwell, points::Measure)
+  map(b -> initialize_state(b, points), obj.branches)
+end
+
+function update_state!(obj::GeneralizedMaxwell, states, F, Fn)
+  @assert length(obj.branches) == length(states)
+  map((b, s) -> update_state!(b, s, F, Fn), obj.branches, states)
 end
 
 function initializeStateVariables(obj::GeneralizedMaxwell, points::Measure)
   map(b -> initializeStateVariables(b, points), obj.branches)
 end
 
-function updateStateVariables!(states, obj::GeneralizedMaxwell, Δt, F, Fn)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  obj.Δt[] = Δt
+function updateStateVariables!(states, obj::GeneralizedMaxwell, F, Fn)
   @assert length(obj.branches) == length(states)
   for (branch, state) in zip(obj.branches, states)
-    updateStateVariables!(state, branch, Δt, F, Fn)
+    updateStateVariables!(state, branch, F, Fn)
   end
 end
 
-function Dissipation(obj::GeneralizedMaxwell, Δt)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  obj.Δt[] = Δt
-  Dα = map(b -> Dissipation(b, Δt), obj.branches)
+function Dissipation(obj::GeneralizedMaxwell)
+  Dα = map(Dissipation, obj.branches)
   D(F, Fn, A...) = mapreduce((Di, Ai) -> Di(F, Fn, Ai), +, Dα, A)
 end
 

src/WeakForms/ElectroMechanics.jl

@@ -72,13 +72,13 @@ end
 # Monolithic strategy
 # -------------------
 
-function residual(physicalmodel::ViscoElectricModel, kine::NTuple{2,KinematicModel},  (u, φ), (v, vφ), dΩ, Λ, Δt, un, A)
-    residual(physicalmodel, Mechano, kine, (u, φ), v, dΩ, Λ, Δt, un, A) +
-    residual(physicalmodel, Electro, kine, (u, φ), vφ, dΩ, Λ, Δt, un, A)
+function residual(physicalmodel::ViscoElectricModel, kine::NTuple{2,KinematicModel},  (u, φ), (v, vφ), dΩ, Λ, un, A)
+    residual(physicalmodel, Mechano, kine, (u, φ), v, dΩ, Λ, un, A) +
+    residual(physicalmodel, Electro, kine, (u, φ), vφ, dΩ, Λ, un, A)
 end
 
-function jacobian(physicalmodel::ViscoElectricModel, kine::NTuple{2,KinematicModel}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A)
-    jacobian(physicalmodel, Mechano, kine,(u, φ), du, v, dΩ, Λ, Δt, un, A) +
-    jacobian(physicalmodel, Electro, kine,(u, φ), dφ, vφ, dΩ, Λ, Δt, un, A) +
-    jacobian(physicalmodel, ElectroMechano, kine,(u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A)
+function jacobian(physicalmodel::ViscoElectricModel, kine::NTuple{2,KinematicModel}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, un, A)
+    jacobian(physicalmodel, Mechano, kine,(u, φ), du, v, dΩ, Λ, un, A) +
+    jacobian(physicalmodel, Electro, kine,(u, φ), dφ, vφ, dΩ, Λ, un, A) +
+    jacobian(physicalmodel, ElectroMechano, kine,(u, φ), (du, dφ), (v, vφ), dΩ, Λ, un, A)
 end

src/WeakForms/WeakForms.jl

@@ -36,9 +36,8 @@ end
 
 Calculate the residual using the given constitutive model and finite element functions.
 """
-function residual(physicalmodel::Mechano, km::KinematicModel, u, v, dΩ, Λ=1.0, vars...; kwargs...)
-  haskey(kwargs, :Δt) && @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  _, ∂Ψu, _ = physicalmodel(; kwargs...)
+function residual(physicalmodel::Mechano, km::KinematicModel, u, v, dΩ, Λ=1.0, vars...)
+  _, ∂Ψu, _ = physicalmodel()
   F, _, _   = get_Kinematics(km; Λ=Λ)
   ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', vars...)))dΩ
 end
@@ -48,9 +47,8 @@ end
 
 Calculate the jacobian using the given constitutive model and finite element functions.
 """
-function jacobian(physicalmodel::Mechano, km::KinematicModel, u, du, v, dΩ, Λ=1.0, vars...; kwargs...)
-  haskey(kwargs, :Δt) && @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  _, _, ∂Ψuu = physicalmodel(; kwargs...)
+function jacobian(physicalmodel::Mechano, km::KinematicModel, u, du, v, dΩ, Λ=1.0, vars...)
+  _, _, ∂Ψuu = physicalmodel()
   F, _, _    = get_Kinematics(km; Λ=Λ)
   ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', vars...)) ⊙ ∇(du)'))dΩ
 end
@@ -87,16 +85,14 @@ function residual(physicalmodel::ThermoElectroMechano, ::Type{Electro}, kine::NT
   return -1.0*∫(∇(vφ)' ⋅ (∂Ψφ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)))dΩ
 end
 
-function residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0, Δt=0.0, vars...)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
+function residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0, vars...)
   κ = physicalmodel.thermo.κ
-  D, ∂D = Dissipation(Δt)
+  D, ∂D = Dissipation()
   return ∫(κ * ∇(θ) ⋅ ∇(vθ) -D(u, φ, θ, vars...))dΩ
 end
 
-function transient_residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), (un, φn, θn), vθ, dΩ, Λ=1.0, Δt=0.0, vars...)
-  @warn "The argument 'Δt' will be removed shortly. Just kept to avoid breaking benchmarks..."
-  DΨ = physicalmodel(Δt=Δt)
+function transient_residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), (un, φn, θn), vθ, dΩ, Λ, Δt, vars...)
+  DΨ = physicalmodel()
   η = -DΨ[4]
   return ∫((1/Δt)*(θ*η(F,E,θ,vars...) - θn*η(Fn,En,θn,vars...) + η(F,E,θ)*(θ - θn)))dΩ
 end

test/TestConstitutiveModels/ElectroMechanicalTests.jl

@@ -117,7 +117,8 @@ end
   Uvn *= det(Uvn)^(-1/3)
   λvn = 1e-3
   Avn = VectorValue(Uvn..., λvn)
-  Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = model(Δt=0.01)
+  update_time_step!(model, 0.01)
+  Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = model()
   @test norm(∂Ψu(F(∇u), E(∇φ), F(∇un), Avn)) ≈ 25.049301121178615
   @test norm(∂Ψuu(F(∇u), E(∇φ), F(∇un), Avn)) ≈ 3110.7607787445168
 end
@@ -130,7 +131,7 @@ end
   viscous_branch2 = ViscousIncompressible(short_term, τ=60.)
   visco_elastic = GeneralizedMaxwell(hyper_elastic, viscous_branch1, viscous_branch2)
   dielectric = IdealDielectric(ε=1.0)
-  model =dielectric+visco_elastic
+  model = dielectric+visco_elastic
   Ke=Kinematics(Electro,Solid)
   E = get_Kinematics(Ke)
   K=Kinematics(Mechano,Solid)
@@ -139,7 +140,8 @@ end
   Uvn *= det(Uvn)^(-1/3)
   λvn = 1e-3
   Avn = VectorValue(Uvn..., λvn)
-  Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = model(Δt=0.01)
+  update_time_step!(model, 0.01)
+  Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = model()
   @test norm(∂Ψu(F(∇u), E(∇φ), F(∇un), Avn, Avn)) ≈ 25.102080194257017
   @test norm(∂Ψuu(F(∇u), E(∇φ), F(∇un), Avn, Avn)) ≈ 3110.9722775475557
 end

test/TestConstitutiveModels/ViscousModelsTests.jl

@@ -70,7 +70,8 @@ end
 
 
 function test_viscous_derivatives_numerical(model; rtolP=1e-12, rtolH=1e-12)
-  Ψ, ∂Ψu, ∂Ψuu = model(Δt = 1e-2)
+  update_time_step!(model, 1e-2)
+  Ψ, ∂Ψu, ∂Ψuu = model()
   F = TensorValue(1.:9...) * 1e-3 + I3
   Fn = TensorValue(1.:9...) * 5e-4 + I3
   Uvn = isochoric_F(TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3)
@@ -94,7 +95,7 @@ end
 
 
 struct EmptyElastic <: Elasto
-  function (::EmptyElastic)(Λ=0.0)
+  function (::EmptyElastic)()
     Ψ(F) = 0.0
     ∂Ψu(F) = TensorValue(zeros(9)...)
     ∂Ψuu(F) = TensorValue(zeros(81)...)
@@ -130,7 +131,8 @@ end
 
 @testset "ViscousIncompressible2" begin
   visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=1.0), τ=10.0)
-  Ψ, ∂Ψu, ∂Ψuu = visco(Δt = 0.1)
+  update_time_step!(visco, 0.1)
+  Ψ, ∂Ψu, ∂Ψuu = visco()
   F    =  1e-2*TensorValue(1,2,3,4,5,6,7,8,9) + I3
   Fn   =  5e-3*TensorValue(1,2,3,4,5,6,7,8,9) + I3
   Fvn  =  2e-2*TensorValue(1.0,2.0,3.0,4.0,5.0,8.7,6.5,4.3,6.5) + I3