Data-Driven Simulation

[LoveFrootLoops]

Hi,
After addressing some initial challenges, I must say that I love working with it. To get started, I explored the "Linear Elasticity" example available at https://gridap.github.io/Tutorials/dev/pages/t003_elasticity/.

Now, I am working on implementing a fascinating concept called "Data Driven Computational Mechanics." Instead of relying on Hookes Law, I define a dataset $\mathcal{D}$ comprising pairs of strain-stress $(\hat{\epsilon}_i, \hat{\sigma}_i)$. My objective is to find the material state $(\epsilon, \sigma)$ that best matches the dataset while satisfying the equilibrium and kinematics equations.

The optimization problem can be formulated as follows:

$\min\limits_{(\hat{\epsilon}_i, \hat{\sigma}_i) \in \mathcal{D}} \frac{1}{2}\int (\hat{\epsilon}_i - \epsilon)^2 + (\hat{\sigma}_i - \sigma)^2 ,d\Omega$

subject to: $\mathrm{div} \sigma = f$ and $\epsilon = \nabla^{\text{sym}} u$.

While defining the equations is not an issue, my challenge lies in finding the data point in the dataset that yields the minimum distance from the material state during the simulation. I have tried using the following code, but it seems I cannot perform calculations directly on a CellField. Can somebody maybe help me regarding this?

using Gridap, PDENLPModels, Random, Distributions

### Material and Dataset
# Material parameters
const E = 70.0e9
const ν = 0.3
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))
σmat(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

# Helper functions
function ToVoigt_sig(A) 
    return Gridap.VectorValue(A.data[1], A.data[2], A.data[3])
end

function ToVoigt_eps(A) 
    return Gridap.VectorValue(A.data[1], A.data[2], 2*A.data[3])
end

# Dataset
ndata = 100
eps = [TensorValue(rand(Uniform(-0.01,0.01), 2, 2)) for i in 1:ndata]
eps_vec = [ToVoigt_eps(epsi) for epsi in eps]
sig_vec = [ToVoigt_sig(σmat(epsi)) for epsi in eps]


### Domain and Spaces
# Domain
domain = (0,1,0,1)
partition = (10, 10)
model = CartesianDiscreteModel(domain,partition)
writevtk(model,"model")
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri_0",[1,3,7])
add_tag_from_tags!(labels,"diri_1",[2,4,8])
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)


# Solution space u
order = 1
reffeᵤ = ReferenceFE(lagrangian, VectorValue{2,Float64}, order)
Vᵤ = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags = ["diri_0", "diri_1"])
disp_x = 0.01
g1 = VectorValue(0.0,0.0)
g2 = VectorValue(disp_x, 0.0)
U = TrialFESpace(Vᵤ,[g1,g2])

# Solution space σ
reffeₛ = ReferenceFE(lagrangian, VectorValue{3, Float64}, order)
Vₛ = TestFESpace(model, reffeₛ; conformity=:H1)
Σ = TrialFESpace(Vₛ)

# Multifield Space
X = MultiFieldFESpace([Vᵤ, Vₛ])
Y = MultiFieldFESpace([U, Σ])

### Simulation
# Objective Function (distance minimazation)
function f(u, σ)
    ## Here should be a minimizer to find the optimal data point in data set D = (eps_vec, sig_vec)
    1/2*∫(*(eps_vec - (ToVoigt_eps∘ε(u)))⋅(eps_vec  - (ToVoigt_eps∘ε(u))) + (sig_vec  - σ)⋅(sig_vec  - σ))*dΩ 
  end

# Constraints (div σ = f)
function res(u, σ, v)
    ∫(σ ⋅ (ToVoigt_eps∘ε(v)))*dΩ 
end

# Solution simulation
op = FEOperator(res, Y, Vᵤ)
ndofs = Gridap.FESpaces.num_free_dofs(Y)
xin = zeros(ndofs)
nlp = GridapPDENLPModel(xin, f, Ω, U, Σ, Vᵤ, Vₛ, op, name = "Test")

[JordiManyer]

@LoveFrootLoops The main issue I see is that sig_vec - σ and eps_vec - (ToVoigt_eps∘ε(u)) make no real sense. Let me explain:

Within the integral, u and σ are CellFields, i.e objects which contain the fields to evaluate for each cell of the mesh. These can be evaluated on CellPoints, i.e objects that contain the quadrature points for each cell of the mesh. The ∫(...)*dΩ syntax will return a DomainContribution, which will contain the sub-integrals in each cell of the mesh.

So there are two things you should probably do:

1. I imagine you want to evaluate the integral, which means you need to add all cell contributions to obtain a scalar quantity. This can be done by using sum(∫(...)*dΩ). This will yield a scalar for a single data point.

2. What I imagine you want to do is to evaluate this integral for each scalar contained in (eps_vec,sig_vec), but that is something Gridap is not prepared to do automatically. The simplest way to do this is to make a for loop that evaluates the integral for each individual data point, i.e

i = 1
for (eps,sig) in zip(eps_vec,sig_vec) 
  res[i] = sum(∫(...)*dΩ)
  i = i+1
end

This is inefficient, since it evaluates u and σ for each data point. But it's by far the easiest.

[fverdugo]

Hi! Very interesting application btw. If you need to evaluate your fe stress and strain at arbitrary points, take a look into tutorial 15 in our tutorial page.

[LoveFrootLoops]

Thanks for the fast and detailed reply.

First of all @fverdugo I cannot find tutorial 15. Can you maybe share a link?

Secondly, @JordiManyer, thanks for the detailed answer and explanations. It helps a lot. But there is one issue which maybe caused some misunderstanding. Let's say I operate on Gauss point level. For each Gauss point I have a corresponding strain and stress pair. For this strain and stress pair I want to find the data point in the data set which is closest to this pair. Then I want to sum the distance over all gauss points.

In your code I would just sum over all data points. But this is not what I try to do. So in pseudo code I would yes to have something like

For each gp in GP
# Find data point in data set closest to strain and stress at the gp
# save the distance (which is defined by the integral)
End
res= Sum of all distances

Is this somehow achievable?

[JordiManyer]

@LoveFrootLoops Manipulating directly at gauss-point level need to be done through composition, just like you do for ToVoigt_sig:

You should define your distance functions as:

function d_eps(ε)
  return minimum(map(eps -> (eps - ToVoigt_eps(ε)) ⋅ (eps - ToVoigt_eps(ε)) , eps_vec))
end

and then compose it to with CellField by calling d_eps∘(ε(u)). Similarly for σ.

Note that within the function d_eps, ε is the already evaluated CellField at a single gauss-point (and therefore is a Number subtype depending on the FE Space). This is why there is no longer need to compose with ToVoigt_eps and you can simply call the function.

But just so that we are clear: This will potentially choose a different data point each for gauss-point (which makes no sense in my mind). If the choice of the data-point needs to be global to the mesh, my first solution is what you need...