Historic variables in Gridap

[santiagobadia]

I have created a driver for isotropic damage model. It is not working yet, since we still have to decide how to implement the historic variables within Gridap. But I think that the essence of the method and design issues are clear

module DamageModel

using Gridap
using LineSearches: BackTracking

##

const E = 2.1e4
const ν = 0.3
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))

# Damage model input
const σ_u = 5.5
const r_0 = σ_u / sqrt(E)
const H = 0.1
q(r) = r_0 + H*(r-r_0)

# Here we should create FEFunctions for the historic variables r and d

@law σ(x,ε) = λ*tr(ε)*one(ε) + 2*μ*ε
@law τ(x,ε) = sqrt(inner(ε(u),σ(x,ε)))

@law σ_new(x,ε) = (1-d)*σ(x,ε)
@law dσ_new(x,dε,ε) = (1-d)*σ(x,dε)

function damage_model(x,u,r,d)

  τ_new = τ(x,ε(u))
  r_current = r(x)
  d_current = d(x)

  if (τ_new <= r_current) # elastic regime

    r_new = r_current # update r at this integration point
    d_new = d_current

    @law σ_new(x,ε) = (1-d_new)*σ(x,ε)
    @law dσ_new(x,dε,ε) = (1-d_new)*σ(x,dε)

  else # inelastic regime

    r_new = τ_new
    d_new = 1 - q(r_new)/r_new
    λ_new = (1-d_new)*λ
    μ_new = (1-d_new)*μ

    @law σ_new(x,ε) = (1-d_new)*σ(x,ε)
    @law dσ_new(x,dε,ε) = (1-d_new)*σ(x,dε) - ((q(r_new)- H*r_new)/(r_new^3))*inner(σ(x,ε),dε)*σ(x,ε)
  end

  return σ_new, dσ_new

end

res(u,v) = inner( ε(v), σ_new(ε(u)) )

jac(u,v,du) = inner( ε(v), dσ_new(ε(u)) )

# The rest shoud be quite standard

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(20,20))

labels = FaceLabels(model)
add_tag_from_tags!(labels,"diri_0",[1,3,7])
add_tag_from_tags!(labels,"diri_1",[2,4,8])

order = 2
diritags = ["diri_0", "diri_1"]

T = VectorValue{2,Float64}

fespace = FESpace(reffe=:Lagrangian, conformity=:H1, valuetype=VectorValue{2,Float64},
model=model, labels=labels, order=order, diritags=["diri_0","diri_1"])

V = TestFESpace(fespace)

# Setup integration
trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

# Setup weak form terms
t_Ω = NonLinearFETerm(res,jac,trian,quad)

# Setup non-linear solver
nls = NLSolver(
show_trace=true,
method=:newton,
linesearch=BackTracking())

solver = NonLinearFESolver(nls)

disp_max = 0.05
disp_inc = 0.05
nsteps = ceil(Int,abs(disp_max)/disp_inc)

x0 = zeros(Float64,num_free_dofs(fespace))

disp_x = disp_max

g0 = zero(T)
g1 = VectorValue(disp_x,0.0)
U = TrialFESpace(fespace,[g0,g1])

op = NonLinearFEOperator(V,U,t_Ω)

println("\n+++ Solving for disp_x $disp_x in step $step of $nsteps +++\n")

l2(u) = inner(u,u)

uh = FEFunction(U,x0)

norm = zeros(Float64,3)
for i in 1:2
  global d
  solve!(uh,solver,op)
  # norm[i] = sqrt(sum( integrate(l2(uh),trian,quad) ))
end

# rf = joinpath(pwd(),"tmp/DamageModel")

# writevtk(trian,"$rf/results_$(lpad(step,3,'0'))",cellfields=["uh"=>uh,"sigma"=>σ(∇(uh))])

##

end #module

[santiagobadia]

As far as I understand, @law takes as input (x,du,u) where u is a FEFunction (linearisation point) and du will be the TrialFESpace for the problem at hand. However, as one can see in the previous example, we have historic variables (that depend on the load process) d and r that are not part of the FE unknowns (they are a postprocess) but change and modify the problem within the nonlinear iterations.

How can we include this update of the res and jac within the nonlinear solver?

One can think about creating some kind of clausure but not sure how to trigger the update.

One could consider that the unknowns are (u,r,d) but it is a mess, I guess. It reminds me to a projection method in which one solves first for u and next for r,d... It can be done but then it is hard to hide the nonlinear solver. We should probably need more expressivity in the nonlinear solver.

By the way, it is mandatory in nonlinear solid mechanics to provide the value of sigma and dsigma at the same time, since most of the steps are shared and they are computed when advancing the historic variables after a load step. It would be quite messy not to do it this way.

We can discuss alternatives tomorrow

[fverdugo]

Could you summarize the equations to be solved? perhaps a scanned note?

[santiagobadia]

The equations to be solved are the standard nonlinear solid mechanics eq's

The historic variables are r and d and the update for computing sigma and dsigma is probably clearer in the code than the notes...

[fverdugo]

Essentially, we need to decide which are the unknowns of the problem, and define the "algebraic" non-linear operator associated with it. I.e., the operator that given a vector for unknowns (possibly having historic variables as well), it provides the corresponding residual (a vector) and jacobian (a sparse matrix). For FE computations, we were able (up to now) to hide this "algebraic" operator inside a "FEOperator". Perhaps for the case of historic variables, we need different way of constructing the "algebraic" operator.

[fverdugo]

By seeing the notes, it is perhaps better to discuss this face to face...

[santiagobadia]

Yes, one option is to define the problem in terms of the unknowns u r and d, where u is the one that will be solved in the algebraic system and d and t are nothing but functions of u that also depend on the previous step of the historic variable, e.g., r_new = f(u,r_old). It is a Gauss point (or node) evaluation.

In that case, my doubt about how to express the law vanishes and the problem is how to efficiently express the problem

u, r = x_old
jac(x_old,dx,v) = -res(x_old,v) # only solve for du,  nonlinear problem...
du, dr = dx
u = u+du
r = f(u,r)

[santiagobadia]

Ok, but the notes are not that important, the point is how to express the problem... the damage model is just an example.

We agree that it is better to talk tomorrow face-to-face.

[santiagobadia]

The situation is even more complicated, the expression of the stress tensor and the constitutive tensor do not have a closed form in terms of the historic and state variables but they are also Gauss point dependent depending on the situation (e.g., elastic or inelastic regime). I have tried to express the idea in the following pseudo-code. We can start working on this idea.

const E = 2.1e4
const ν = 0.3
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))

# Damage model input
const σ_u = 5.5
const r_0 = σ_u / sqrt(E)
const H = 0.1
q(r) = r_0 + H*(r-r_0)

# Here we should create FEFunctions for the historic variables r and d

@law σ(x,ε) = λ*tr(ε)*one(ε) + 2*μ*ε
@law τ(x,ε) = sqrt(inner(ε(u),σ(x,ε)))

function residual_and_jacobian_at_gauss_point(x_old,x_current,y,dx) # evaluated at a Gauss point

  du = x[1]
  v = y[1]

  u, r_current, d_current = x_current
  tmp, r_old, d_old = x_old
  τ_current = τ(x,ε(u))

  if (τ_current <= r_old) # elastic regime

    r_current = r_old # update r at this integration point
    d_current = d_old

    @law σ_new(x,ε) = (1-d_current)*σ(x,ε)
    @law dσ_new(x,dε,ε) = (1-d_current)*σ(x,dε)

  else # inelastic regime

    r_current = τ_current
    d_current = 1 - q(r_current)/r_current

    @law σ_new(x,ε) = (1-d_current)*σ(x,ε)
    @law dσ_new(x,dε,ε) = (1-d_current)*σ(x,dε) - ((q(r_current)- H*r_current)/(r_current^3))*inner(σ(x,ε),dε)*σ(x,ε)

  end

  x_current = tmp, r_current, d_current

  res = inner( ε(v), σ_new(ε(u)) )
  jac = inner( ε(v), dσ_new(ε(du),ε(u)) )

  return x_current, res, jac

end

##

for load_stages

  u_old, r_old, d_old = x_old # Values from previous time step, i.e., initial state

  x_current = x_old # Start nonlinear iterations with previous time step values

  for nonlinear iterations

    nonlinear_solid_mechanics!(x_current,x_old) 
  # using the previous Gauss point-based expression of the forms and updating inside the
  # historic variables

  end

  x_old = x_current

end

##

[fverdugo]

Let us assume that we want to solve a mechanical problem with a body force, - div(sigma) = f.

At a generic point (e.g., a gauss point), the relation between the symmetric gradient of the displacement and the corresponding stress is given by a julia function with the following signature

σ_gp, state_gp_new = constitutive_relation_gp(ε_gp, state_gp_old)

Note that this relation is written in terms of a state variable, which is updated from state_gp_old to state_gp_new. Further evaluations of the constitutive relation need to be evaluated with the most recent version available of the state variable.

In addition, we have a linearization of this constitutive relation in a julia function with the following signature

σ_gp, dσ_gp, state_gp_new = constitutive_relation_gp(ε_gp, dε_gp, state_gp_old)

Now, using the constitutive relation at a generic point, we build a related function which takes cell-wise fe-related quantities, e.g. fe functions, cell fields, etc. The julia function will have this signature

 σ_fe = constitutive_relation_fe!(state_fe,ε_fe)

Here, state_fe is an object that represents the state variable over the entire mesh, it can be a fe function or a raw array (the user decides). The important thing is that state_fe has a mutable state that changes when constitutive_relation_fe! is called (see the ! sign).

We can also implement the corresponding linearization

 σ_fe, dσ_fe = constitutive_relation_fe!(state_fe,ε_fe,dε_fe)

With this, the corresponding res and jac functions representing the weak residual and jacobian are written as follows:

state_fe = # the user defines the inital state before functions res and jac

function res(u,v)
  ε_fe = ε(u)
  σ_fe = constitutive_relation_fe!(state_fe, ε_fe)
  inner(ε(v), σ_fe) - inner(v,f)
end

function jac(u,v,du)
  ε_fe = ε(u)
  dε_fe = ε(du)
  σ_fe, dσ_fe = constitutive_relation_fe!(state_fe, ε_fe, dε_fe)
  inner(ε(v), dσ_fe)
end

Remarks

1. The user needs to implement constitutive_relation_fe! by using the relation at a generic point constitutive_relation_gp.

2. Currently, the residual and jacobian are evaluated separately, but this can change in the future.

[santiagobadia]

Can you provide more details about the usage of constitutive_relation_gp! (I put ! to stress the fact that the update of the historic variable is within this function) within constitutive_relation_fe!, e.g. using some silly example.

[fverdugo]

Can you provide more details about the usage of constitutive_relation_gp! (I put ! to stress the fact that the update of the historic variable is within this function)

Note that constitutive_relation_gp has to be a pure function (no side effects on the inputs) since all inputs will be typically (inmutable) numbers. The function gets the old state and returns the new one.

Here an example of how to implement constitutive_relation_fe! from constitutive_relation_gp.

struct StateFE
  state_dofs
  state_space
  σ_dofs
  σ_space
end

function constitutive_relation_fe!(state_fe::StateFE, ε_fe)

  σ_space = state_fe.σ_space
  state_dofs = state_fe.state_dofs
  σ_dofs = state_fe.σ_dofs

  ε_dofs = free_values(interpolate(ε_fe, σ_space))

  _fill_sigma_and_state!(state_dofs, σ_dofs, ε_dofs)

  σ_fe = FEFunction(σ_space,σ_dofs)

  return σ_fe
end

function _fill_sigma_and_state!(state_dofs, σ_dofs, ε_dofs)
  for gp in 1:length(σ_dofs)

    ε_gp = ε_dofs[gp]
    state_gp_old = state_dofs[gp]

    σ_gp, state_gp_new = constitutive_relation_gp(ε_gp, state_gp_old)

    σ_dofs[gp] = σ_gp
    state_dofs[gp] = state_gp_new

  end
end

For the linearized version, it would be analogous.

[fverdugo]

For the linearized version, it would be analogous.

The difference is that we need to return a cell basis instead of a cell field.

[santiagobadia]

What is the analogous of σ_fe = FEFunction(σ_space,σ_dofs) for dσ_fe that will work in inner(ε(v), dσ_fe)?

You want a cell array such that at every cell at every Gauss points provides you with a method that applies the constitutive tensor dσ over all the test function, i.e., nshape second order tensors. And the resulting cell array of ngauss x nshape second order tensors should be such that it can be used within inner(ε(v), dσ_fe). I am not saying that it does not work, I just want to make sure how it is done.

[fverdugo]

Yes, we need to build a CellBasis with tensorial value that represents dσ_fe. The problem is that building this would require to go quite deep into the code...far beyond a normal user would master.

To think how to solve it.

[santiagobadia]

OK, at least we have now a better understanding of what we need. I will mark the issue as new functionality.

[fverdugo]

@santiagobadia I have implemented a way of using constitutive laws with state variables.

The main buildings block are two new FETerm types called AffineFETermFromCellMatVec and FETermFromCellJacRes that allows the user to fully control the (simultaneous) integration of elemental matrices and vectors.

E.g., for a Poisson problem on could create the following FE term:

using Gridap
using Gridap.FESpaces # needed to use the advanced function

function cellmatvec!(mat,vec,∇v,∇u,v,j,w)
  Q = length(w)
  M,N = size(mat)
  for q in 1:Q

    dV = det(j[q])*w[q]

    for n in 1:N
      for m in 1:M
        mat[m,n] += ∇v[q,m]*∇u[q,n]*dV
      end
    end

    for m in 1:M
      vec[m] += v[q,m]*dV
    end

  end
end

trian = ...
quad = ...

q = get_coordinates(quad)
w_q = get_weights(quad)
ϕ = get_cell_map(trian)
jac = ∇(ϕ)
jac_q = evaluate(jac,q)
x_q = evaluate(ϕ,q)

function cellmatvec(v,u)
  v_q = evaluate(v,q)
  ∇v_q = evaluate(∇(v),q)
  apply_cellmatvec(cellmatvec!, ∇v_q, ∇v_q, v_q, jac_q, w_q, x_q)
end

t_Ω = AffineFETermFromCellMatVec(cellmatvec,trian)

This will be useful (of course not in the previous case) in situations where the user needs full control (e.g., with state variables).

E.g., a p-laplacian driver with a dummy state variable:

using Test
using LinearAlgebra

using Gridap
import Gridap: ∇

# Needed to use some advanced functions
using Gridap.FESpaces

domain = (1,2,1,2)
partition =(3,3)
model = CartesianDiscreteModel(domain,partition)

u(x) = x[1] + x[2]
∇u(x) = VectorValue(1,1)
∇(::typeof(u)) = ∇u

const p = 3
flux(∇u) = norm(∇u)^(p-2) * ∇u
dflux(∇du,∇u) = (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2)*∇du

order = 1
T = Float64

V = TestFESpace(
 model=model,
 order=order,
 reffe=:Lagrangian,
 conformity=:H1,
 valuetype=T,
 dirichlet_tags="boundary")

U = TrialFESpace(V,u)

degree = 2*order
trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)

# Prepare some quantities at gp level
q = get_coordinates(quad)
w_q = get_weights(quad)
ϕ = get_cell_map(trian)
j_q = evaluate(∇(ϕ),q)

# Init the state variables at gp (a float per gp)
s_q = [ zeros(Float64,size(qi)) for qi in q ]

# Definition of the integration of the jacobian and residual at cell level
function celljacres!(jac,res,∇v,∇du,∇uh,j,w,s)

  Q = length(w)
  M,N = size(jac)
  for q in 1:Q

    dV = abs(det(j[q]))*w[q]
    flux_q = flux(∇uh[q])

    # update the state variable with a dummy expression
    s[q] = max(s[q],norm(flux_q))

    for m in 1:M
      for n in 1:N
        jac[m,n] += ∇v[q,m]*dflux(∇du[q,n],∇uh[q])*dV
      end
      res[m] += ∇v[q,m]*flux_q*dV
    end

  end
end

# Build a FETerm from the definition of the cell jacobian and residual

function celljacres(uh,v,du)
  # Prepare remaining quantities at gp level
  ∇uh_q = evaluate(∇(uh),q)
  ∇du_q = evaluate(∇(du),q) 
  ∇v_q = evaluate(∇(v),q)
  # This builds an array of tupes (jac,res) representing the contribution of each cell
  apply_cellmatvec(celljacres!,∇v_q,∇du_q,∇uh_q,j_q,w_q,s_q)
end

t_Ω = FETermFromCellJacRes(celljacres,trian)

op = FEOperator(V,U,t_Ω)

nls = NLSolver(show_trace=true, method=:newton)
solver = FESolver(nls)

x = rand(T,num_free_dofs(U))
uh0 = FEFunction(U,x)
uh, = solve!(uh0,solver,op)

# Print state variables and solution
x_q = evaluate(ϕ,q)
writevtk(x_q,"x",nodaldata=["s"=>s_q])
writevtk(trian,"trian",nsubcells=order,cellfields=["uh"=>uh])

e = u - uh

l2(u) = inner(u,u)
sh1(u) = inner(∇(u),∇(u))
h1(u) = sh1(u) + l2(u)

el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))
ul2 = sqrt(sum( integrate(l2(uh),trian,quad) ))
uh1 = sqrt(sum( integrate(h1(uh),trian,quad) ))

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

The same can also be done in multi field computations. Eg.

function cellmatvec!(A,B,y,x,j,w)

  A_vu = A[1,1]
  A_vp = A[1,2]
  A_qp = A[2,2]
  B_v = B[1]
  B_q = B[2]
  u = x[1]
  p = x[2]
  v = y[1]
  q = y[2]

  N_v, N_u = size(A_vu)
  N_q, N_p = size(A_qp)
  S = length(w)

  for s in 1:S
    dV = det(j[s])*w[s]
    for n_v in 1:N_v
      for n_u in 1:N_u
        A_vu[n_v,n_u] += v[s,n_v]*u[s,n_u]*dV
      end
      for n_p in 1:N_p
        A_vp[n_v,n_p] += v[s,n_v]*p[s,n_p]*dV
      end
    end
    for n_q in 1:N_q
      for n_p in 1:N_p
        A_qp[n_q,n_p] -=  q[s,n_q]*p[s,n_p]*dV
      end
    end
    for n_v in 1:N_v
      B_v[n_v] += v[s,n_v]*4*dV
    end
    for n_q in 1:N_q
      B_q[n_q] += q[s,n_q]*dV
    end
  end

end

TODO: For the moment, not for integrals on the skeleton

(Work not yet in master)

[santiagobadia]

This development is good, because it provides flexibility. In any case, I will still think we should develop an integration point definition of the residual and Jacobian in the spirit of the discussion above. It would be more in line with the standard usage of the library. Don't you think so?

[fverdugo]

Yes, I have done this in order to have an "advanced" mode of defining cell matrices and vectors which provides full flexibility and also improved performance in some cases since it does not introduce complex operation trees.

In any case, I am thinking about implementing a macro @statelaw which will be the counter part of @law with state variables. This will be the one users are supposed to use.

By the way, I have realized that trying to compute the residual and Jacobin simultaneously is actually counterproductive. The nlsolve function that we use in the non-linear solver computes the residual and Jacobians separately (it only builds them simultaneously once at the first iteration...).

[fverdugo]

I have implemented the @statelaw macro. I think we can close the issue.

Example:

using Test
using LinearAlgebra
using Gridap
import Gridap: ∇

domain = (1,2,1,2)
partition =(3,3)
model = CartesianDiscreteModel(domain,partition)

u(x) = x[1] + x[2]
∇u(x) = VectorValue(1,1)
∇(::typeof(u)) = ∇u

const p = 3

# The function annotated with @statelaw should return also the updated state
# which will be set internally into the corresponding array
@statelaw function flux(∇u,s)
  f = norm(∇u)^(p-2) * ∇u
  new_s = max(s,norm(f))
  f, new_s
end

@law dflux(∇du,∇u) = (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2)*∇du

order = 3
T = Float64

V = TestFESpace(
 model=model,
 order=order,
 reffe=:Lagrangian,
 conformity=:H1,
 valuetype=T,
 dirichlet_tags="boundary")

U = TrialFESpace(V,u)

degree = 2*order
trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)

# Initialize the state variable
s = QPointCellField(0.0,trian,quad)

# Note the call flux(∇(u),s)
res(u,v) = inner( ∇(v), flux(∇(u),s) )
jac(u,v,du) = inner(  ∇(v) , dflux(∇(du),∇(u)) )

t_Ω = FETerm(res,jac,trian,quad)

op = FEOperator(V,U,t_Ω)

nls = NLSolver(show_trace=false, method=:newton)
solver = FESolver(nls)

x = rand(T,num_free_dofs(U))
uh0 = FEFunction(U,x)
uh, = solve!(uh0,solver,op)

# Visualize state variable ad gps
q = get_coordinates(quad)
x = get_physical_coordinate(trian)
x_q = evaluate(x,q)
s_q = evaluate(s,q)
writevtk(x_q,"x",nodaldata=["state"=>s_q])

# Visualize state variable as a field using a L2 projection into the solution space
t_l2 = AffineFETerm((v,u)->v*u, (v) -> v*s, trian, quad)
op_l2 = AffineFEOperator(V,U,t_l2)
sh = solve(op_l2)
writevtk(trian,"trian",nsubcells=order,cellfields=["uh"=>uh,"state"=>sh])

e = u - uh

l2(u) = inner(u,u)
sh1(u) = inner(∇(u),∇(u))
h1(u) = sh1(u) + l2(u)

el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

@test el2 < 1.e-8
@test eh1 < 1.e-7

[fverdugo]

Using the @statelaw, I have started to implement a tutorial on a damage law:

https://github.com/gridap/Tutorials/blob/damage_model_tutorial/src/damage_model.jl

[santiagobadia]

In our last meeting, we have decided that the update of the state variables, which is enforced by using a @statelaw instead of a @law, should be done for the Jacobian. On the other hand, the residual must be a @law. The reason for that is that the residual provides the same result irrespectively of whether the state variable has been updated or not, but this is not the case for the Jacobian. This way, the approach works for whatever implementation of the nonlinear solver.

On the other hand, the damage mode in https://github.com/gridap/Tutorials/blob/damage_model_tutorial/src/damage_model.jl should be updated and completed (there are missing terms at the Jacobian).

[fverdugo]

First damage model working in https://github.com/gridap/Tutorials/blob/damage_model_tutorial/src/damage_model.jl

Thanks to @macaicedos was quite straight forward to set-up the problem

I need to add some test in Gridap that tests this functionality and then I'll close the issue.

[fverdugo]

On the other hand, the damage mode in https://github.com/gridap/Tutorials/blob/damage_model_tutorial/src/damage_model.jl should be updated and completed (there are missing terms at the Jacobian).

Finally, the update is not at the Jacobian level within the non-linear solver. The state variables are updated at the load step level, not within the non-linear solver