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?
Best Answer
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.