How to uses skfem.Basis in coupled problem

[grupp-tim]

Hello @kinnala,

I am trying to implement the Folgar-Tucker-Model (in 2D) which is a phenomenological model describing the evolution of fiber orientation distributions.

$$\frac{\text{D}\mathbf{A}}{\text{D}T} = \mathbf{W} \mathbf{A} - \mathbf{A} \mathbf{W}+ \xi \bigl(\mathbf{D} \mathbf{A} + \mathbf{A} \mathbf{D} - 2\mathbb{A} [\mathbf{D}] \bigr) + 2 C_{\text{I}} \dot{\gamma}(\mathbf{I} - 2 \mathbf{A})$$

With the rate-of-deformation tensor $\mathbf{D}$ and the vorticity tensor $\mathbf{W}$, defined as:

$$D_{ij} = \frac{1}{2} ( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i})$$

$$W_{ij} = \frac{1}{2} ( \frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i})$$

In 2 dimensions, the symmetric fiber orientation tensor (a second order tensor) with trace equal to 1

$$ \begin{pmatrix} a_{11} & a_{12} \\ a_{12} & 1 - a_{11} \end{pmatrix} $$

has only two independent entries, $a_{11}$ and $a_{12}$, which I want to solve for.

The fourth order tensor is approximated by the so called linear closure:
$$A_{ijkl}^{\text{L}} = -\frac{1}{35} (\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}) +\delta_{il}\delta_{jk}) +\frac{1}{7} (A_{ij}\delta_{kl} +A_{ik}\delta_{jl}+A_{il}\delta_{jk} +A_{jl}\delta_{ik} +A_{jk}\delta_{il} +A_{kl}\delta_{ij}$$

I derived the coupled system for $a_{11}$ and $a_{12}$:

$$\frac{\text{D}a_{11}}{\text{D}T} = \frac{\text{d} a_{11}}{\text{d} t} + v_{1} \frac{\partial a_{11} }{\partial x_{1}} + v_{2} \frac{\partial a_{11}}{\partial x_{2}} = a_{12} \left( \frac{\partial v_{1}}{\partial x_{2}} - \frac{\partial v_{2}}{\partial x_{1}} \right)+ \frac{1}{4} \frac{\partial v_{1}}{\partial x_{1}} - \frac{1}{4} \frac{\partial v_{2}}{\partial x_{2}} + 2 C_{\text{I}} \dot{\gamma} - 4 C_{\text{I}} \dot{\gamma} a_{11}$$

$$\frac{\text{D}a_{12}}{\text{D}T} = \frac{\text{d} a_{12}}{\text{d} t} + v_{1} \frac{\partial a_{12} }{\partial x_{1}} + v_{2} \frac{\partial a_{12}}{\partial x_{2}} = a_{11} \left( \frac{\partial v_{2}}{\partial x_{1}} - \frac{\partial v_{1}}{\partial x_{2}} \right)+ \frac{3}{4} \frac{\partial v_{1}}{\partial x_{2}} - \frac{1}{4} \frac{\partial v_{2}}{\partial x_{1}} - 4 C_{\text{I}} \dot{\gamma} a_{12}$$

Which discretized becomes:

$$ \begin{bmatrix} C_{ij} & 0 \\ 0 & C_{ij} \\ \end{bmatrix} \begin{pmatrix} \dot{a_{11j}} \\ \dot{a_{12j}} \end{pmatrix} + \begin{bmatrix} K_{Iij} & K_{IIij} \\ -K_{IIij} & K_{Iij} \end{bmatrix} \begin{pmatrix} a_{11j} \\ a_{12j} \end{pmatrix} = \begin{pmatrix} R_i \\ Q_i \end{pmatrix} $$

The coupling Term $K_{IIij}$ is:

$$\int_{\Omega} N_i N_j \left( \frac{\partial v_{2}}{\partial x_{1}} - \frac{\partial v_{1}}{\partial x_{2}} \right) \text{d} \Omega $$

with $N_i$ and $N_j$ coming from the discretization of the weight function and $a_{11}$, $a_{12}$.

I've created a simple rectangular mesh with height $h$ and length $L$.

m = skfem.MeshTri.init_tensor(
    np.linspace(0, L, int(L) + 1),
    np.linspace(0, h, int(h) + 1),
)

The velocity is defined as:

$$v_1 = V_0 \frac{x_2}{h}$$

$$v_2 = 0$$

which is a simple shear flow with the upper plate moving with velocity $V_0$.

def v1(x):
    v1 = V0 * x[1] / h
    return v1

I found several other discussion here, e.g. #616 or #1028 where for each unknown an own basis is used. Therefore, I created two basis, one for each unknown, with a TriP1 Element:

e = skfem.ElementTriP1()

basis = {
    "a11": skfem.Basis(m, e),
    "a12": skfem.Basis(m, e),
}

The velocities are projected to the mesh and then interpolated at the quadrature points:

v1_proj = basis["a11"].project(v1)
v1_intpol = basis["a11"].interpolate(v1_proj)

v2_proj = basis["a11"].project(v2)
v2_intpol = basis["a11"].interpolate(v2_proj)

The bilinear form is defined as:

@skfem.BilinearForm
def k12(u, v, w):
    return v * u * (w.v2.grad[0] - w.v1.grad[1])

which gets assembled to:

K12 = skfem.asm(
    k12,
    basis["a12"],
    basis["a11"],
    v1=v1_intpol,
    v2=v2_intpol,
)

K21 = skfem.asm(
    k12,
    basis["a11"],
    basis["a12"],
    v1=v1_intpol,
    v2=v2_intpol,
)
...
K = bmat([[K11, K12], [-K21, K22]], "csr")

1. My question now is, how do I use the basis properly for the two unknowns I am solving for? How do I know which order to choose when calling the two bases in the asm command? I have to admit, I havn't fully understood the concept of the basis yet.

2. And which basis do I use to project the velocity to the mesh? For now I just choose the basis["a11"]

3. I want to extend the model by using a so called quadratic closure $A_{ijkl}^{\text{Q}} = A_{ij} A_{kl}$, making the system non-linear. For that I guess I would need to change the implementation to something shown here: #1028 and #1040?

I hope my question becomes clear to you. If not I am happy to extend this question.
Thank you!
Tim

PS: I just started with this project, so I might come along with some more questions :)

[kinnala]

Here you have the same mesh and the same element for both Basis objects so their contents will be exactly the same and it does not matter at all which order you use them or which one you use for projections. In fact, if you want to save memory, define only one basis and use it everywhere. However, were you decide to use different element for basis["a12"] and basis["a11"] then it becomes important.

[kinnala]

Nonlinear problems have to be linearized before they can be put into scikit-fem which is for linear problems. There are many ways to linearize a given PDE problem. They go imto two general categories:

1. Newton(-type) linearization where a (Frechet) derivative is used to find a perturbation towards the correct solution given the current approximation
2. Fixed point linearization where some suitable parts of the nonlinear term are taken as fixed values from the previous iteration

The first approach requires calculating the derivative of the nonlinear terms. For this purpose we have the beta quality skfem.experimental.autodiff which uses automatic differentiation to calculate the derivative. The derivative can be also calculated by hand.

The second approach is often easier to code and understand but depending on the problem, the nonlinear iteration might converge more slowly. Example in docs: https://scikit-fem.readthedocs.io/en/latest/howto.html#discrete-functions-in-forms

[grupp-tim]

Thank you very much for your help.
I changed it to just one basis as you suggested and it works fine, using a simple Euler backwards scheme to solve the system.

basis = skfem.Basis(m, e)
...
v1_proj = basis.project(v1)
v1_intpol = basis.interpolate(v1_proj)
...
K11 = skfem.asm(
    k11,
    basis,
    v1=v1_intpol,
    v2=v2_intpol,
)

I checked the assembly process by looking at the matrices, and it seems good. You can see the four quadrants (divided at $x,y = 55$) corresponding to $a_{11}$, $a_{12}$ and the coupling quadrants.

plt.imshow(C.toarray(), interpolation="none")