Hi,

I’m in the process of solving a Landau theory model concerning a crystal where strain and electric polarization are coupled. I’m trying to minimize the total free energy given by:

```
F = Fl(P, α) + Fc(ε, C) + Fq(P, ε, q) + Fg(∇P, G)
function Fl(P::JuAFEM.Vec{3}, α)
P² = P.^2
P⁴ = P.^4
return α[1, 1, 1] * sum(P²)
+ α[1, 1, 1] * sum(P⁴)
+ α[1, 1, 2] * sum(P.^6)
+ α[1, 2, 1] * (P²[1]*P²[2] + P²[2]*P²[3] + P²[1]*P²[3])
+ α[1, 1, 2] * (P⁴[1]*(P²[2] + P²[3]) + P⁴[2]*(P²[1] + P²[3]) + P⁴[3]*(P²[1] + P²[2]))
+ α[1, 2, 3] * prod(P²)
end
Fc(ε::Tensor{2, 3}, C::Tensor{4, 3}) = 0.5ε ⊡ C ⊡ ɛ
Fq(P::Vec{3}, ε::Tensor{2, 3}, q::Tensor{4, 3}) = -P ⋅ (ε ⊡ q) ⋅ P
Fg(∇P::Tensor{2, 3}, G::Tensor{4, 3}) = 0.5∇P ⊡ G ⊡ ∇P
```

`ε`

is the gradient of the displacement field in all directions, `P`

is the polarization, `C`

gives the penalty for creating strain, `q`

the coupling between strain and polarization `P`

, `G`

the Ginzburg coefficients for spatial variation of `P`

and `α`

the different Landau components for free energy from `P`

.

Since we are dealing with 6 variables `u_1,2,3, P_1,2,3`

in 2D, on a mesh of size `> 2000x600`

, throwing it into `Optim.jl`

and brute force optimize it is pretty much impossible. Next option is to analytically minimize the free energy with regards to the two fields (`u`

, `P`

) and solve the equations on a FEM model. This results in a system of nonlinear equations.

Now we come to the questions. Looking at the examples of `JuAFEM`

I can almost understand how to solve a nonlinear model, but since I’m not really an expert on FEM I don’t fully understand how to connect `JuAFEM`

to a nonlinear solver like `NLsolve`

, is this possible? And if so, how does one go about it? The second question is more about how to use `ForwardDiff`

, can I exploit this package to not have to analytically compute the derivatives of `F`

, which result in the nonlinear equations? This would be a great aid, because the derivation is pretty cumbersome and prone to errors. The main thing I don’t fully understand is if it’s possible to derive the `Fq`

equation, which couples the two fields, through the gradient `ε`

. From what I read, `ForwardDiff`

only works when a function can be called with one `Vector`

input or basically 1 parameter, can I use a tuple to derive this? Will that impact performance?

Is there another suite of packages I should have a look at, or that encompass everything I want to do?

Thanks in advance for any comments/ help!

Cheers