# Finding the weak form of a PDE with a tensor argument

Let $$A$$ be a $$3×3$$ matrix, a (possibly) complex function of $$x, z \in \mathbb{R}$$ representing the Order Parameter in the Ginzburg Landau equations, $$A = \begin{pmatrix} A_{uu} & A_{uv} & A_{uw} \\ A_{vu} & A_{vv} & A_{vw} \\ A_{wu} & A_{wv} & A_{ww} \end{pmatrix}.$$ $$K_i$$,$$\beta_i$$, and $$\alpha (T)$$ are known parameters. I need to solve the following equation, (with the implied sum over $$j$$, for each component identified by $$\mu$$, $$i$$), $$K_1 \partial^2_j A_{\mu i} + K_{23} \partial_i (\partial_j A_{\mu j}) = 2 \beta_1 Tr(AA^T) A_{\mu i} + 2 \beta_2 Tr(AA^\dagger)A_{\mu i} + 2 \beta_3[AA^TA]_{\mu i} + 2 \beta_4[AA^\dagger A]_{\mu i} + 2 \beta_5[AA^TA]_{\mu i} + \alpha(T) Tr(AA^\dagger) = (\text{rhs})_{\mu i}$$

I have a code in C++ that implements the FDM with relaxation, but we have found that our mixed derivative approximations have large enough error that the code doesn't always converge to a solution. We are looking at the option of using others' FEM solvers (like FreeFEM or MOOSE), but I'm having a hard time getting the weak form of our set of equations. I am following the description here, and obtained $$0 = K_1 \oint_{\Gamma} {\psi \vec{\nabla}A_{\mu i} \cdot \hat{n}} - K_1 \int_{\Omega} {\vec{\nabla} \psi \cdot \vec{\nabla}A_{\mu i}} - \int_{\Omega} {\psi (\text{rhs})_{\mu i}} + K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})},$$ but I don't know how to get the last term. I saw this post and thought that maybe I'd have to explicitly write out all 9 equations and use 9 test functions?

Edit: I played around with the equations a little more and figured I could write them generally as, $$(\text{rhs})_{\mu i} = K_1 \vec{\nabla}^2 A_{\mu i} + K_{23} \partial_i \Big[ \vec{\nabla} \cdot \vec{A}_{\mu;r} \Big],$$ where $$\vec{A}_{\mu;r}$$ is the vector formed from the $$\mu$$-th row of $$A$$ (is there a better way to represent that?). Thus, $$0 = K_1 \oint_{\Gamma} {\psi \vec{\nabla}A_{\mu i} \cdot \hat{n}} - K_1 \int_{\Omega} {\vec{\nabla} \psi \cdot \vec{\nabla}A_{\mu i}} - \int_{\Omega} {\psi (\text{rhs})_{\mu i}} + K_{23}\int_{\Omega} {\psi \partial_i \Big[ \vec{\nabla} \cdot \vec{A}_{\mu;r} \Big]},$$ which might be easier (or more obvious how) to integrate by parts?

So, would I still have to use 9 test functions? We expect the solution function for each element to be different, so I might still need to have all 9 separately...

Edit 2: Taking the last term, $$K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})}$$, explicitly writing out the sum over $$j$$, and using the product rule, we can say, $$K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})} \to K_{23} \sum_{j=u,v,w}{ \Bigg[ \int_{\Omega}{\partial_i (\psi \partial_j A_{\mu j})} - \int_{\Omega}{(\partial_i \psi)(\partial_j A_{\mu j})} \Bigg] }$$

The first term here ($$\int_{\Omega}{\partial_i (\psi \partial_j A_{\mu j})}$$) looks like the divergence term that we previously rewrote as a boundary integral...except that there are no vectors. Is this now in the right form so that I can use it in a FEM solver?

• You would be better off asking this on the scientific computing stackexchange. Jun 6, 2022 at 22:41
• You need test functions also to be matrices, i.e. $\psi_{\mu i}$ and need summation by $i, \mu$ Jun 7, 2022 at 12:29

To start with I strongly suggest not to mix vector and index notation. Let's start with indexed version of the integration by parts formula: $$\int_{\Omega} (\partial_i u_{\alpha,\beta,\dots}) v_{\lambda,\mu,\dots} d\Omega = -\int_{\Omega} u_{\alpha,\beta,\dots} (\partial_i v_{\lambda,\mu,\dots}) d\Omega + \int_{\Gamma} n_i u_{\alpha,\beta,\dots} v_{\lambda,\mu,\dots} d\Gamma. \tag{*}$$ Indices $$\alpha, \beta, \dots$$ and $$\lambda, \mu, \dots$$ are arbitrary and may include $$i$$. Note that in the last term $$\int_\Omega \partial_i$$ became $$\int_\Gamma n_i$$.
The strong form of the equation is: $$K_1 \partial_j^2 A_{\mu i} + K_{23} \partial_i (\partial_j A_{\mu j}) = ({\rm rhs})_{\mu i}.$$ Weak form is obtained by multiplying with $$\psi_{\mu i}$$, summation by $$\mu, i$$ and integration over $$\Omega$$: $$K_1 \int_{\Omega} \psi_{\mu i} \partial_j^2 A_{\mu i} d\Omega + K_{23} \int_{\Omega} \psi_{\mu i} \partial_i (\partial_j A_{\mu j}) d\Omega = \int_{\Omega} \psi_{\mu i} ({\rm rhs})_{\mu i} d\Omega.$$ This form is not good enough since it contains second order derivatives. Let's eliminate them using (*) $$-K_1 \int_{\Omega} (\partial_j \psi_{\mu i}) (\partial_j A_{\mu i}) d\Omega -K_{23} \int_{\Omega} (\partial_i \psi_{\mu i}) (\partial_j A_{\mu j}) d\Omega+\\ +K_1 \int_{\Gamma} \psi_{\mu i} n_j (\partial_j A_{\mu i}) d\Gamma +K_{23} \int_{\Gamma} \psi_{\mu i} n_i (\partial_j A_{\mu j}) d\Gamma = \int_{\Omega} \psi_{\mu i} ({\rm rhs})_{\mu i} d\Omega.$$
Note that $$\int_{\Omega}$$ terms are symmetric by $$\psi \leftrightarrow A$$, which means that their discretization also would be a symmetric matrix.
You did not specify your boundary conditions, they are required to process further the $$\int_\Gamma$$ terms. In case of Dirichlet boundary conditions they simply vanish (due to $$\psi\big|_\Gamma = 0$$). In Neumann case they would be $$K_1 n_j (\partial_j A_{\mu i}) + K_{23} n_i (\partial_j A_{\mu j}) = g_{\mu i}$$ and the integrals migrate to the right hand side as $$\int_\Gamma \psi_{\mu i} g_{\mu i} d\Gamma$$ term.
• This is very helpful! Just to clarify a few things: the $K_1$ integral over $\Gamma$ should have $d\Gamma$? The $\int_{\Omega}$ terms aren't quite symmetric since one $A$ has an $i$ index while the other has $j$? For mixed boundary conditions, the integrals can be trivially broken up to apply different kinds on different domain boundaries? And the $n_i$ represents what, in this case with matrices? Jun 7, 2022 at 17:03
• That was a copy-paste typo, coorected it, thank you. Each of the terms $\int_\Omega$ is symmetric with respect to interchanging of $A$ and $\psi$. It might not be so obvious for the second, but note that $i$ and $j$ are dumb indices. Mixed boundary conditions are similar to Neumman, you eliminate derivatives using boundary conditions. But in this case there remain terms like $\int_\Gamma \psi A d\Gamma$ in the left hand side. The vector $n_i$ is the unit outward normal to $\Omega$, it depends only on the domain and not on the trial or test function. Jun 7, 2022 at 17:42