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

I am trying to solve for the order parameter ($$A$$) in the Ginzburg Landau equations. I had asked on the math SE site but was recommended to ask here. We are trying to solve the following equation, (with the implied sum over $$j$$, for each component identified by $$\mu$$, $$i$$), $$K_i$$, $$\beta_i$$, and $$\alpha (T)$$ are known parameters (constants wrt position), $$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},$$ assuming, $$A = \begin{pmatrix} A_{uu} & A_{uv} & A_{uw} \\ A_{vu} & A_{vv} & A_{vw} \\ A_{wu} & A_{wv} & A_{ww} \end{pmatrix}.$$

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?

My main question is: can this system of equations be reasonably written for a FEM solver? If so, how?

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?). (A note here: the rows and columns of the order parameter do not necessarily correspond to spatial variables, and for our ultimate goal will be in cylindrical coordinates; so it may not be a good idea to represent this as a vector.) 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...

Trying a slightly different approach: 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?

After applying the product rule, I found that I could select a different order of derivatives (because these derivative operators commute, right?) and instead obtain, $$\to K_{23} \Bigg[ \int_{\Omega}{ \vec{\nabla} \cdot (\psi \partial_i \vec{A}_{\mu;r})} - \sum_{j=u,v,w}{ \int_{\Omega}{(\partial_j \psi)(\partial_i A_{\mu j})} \Bigg] }$$ Now the first term could be rewritten using the divergence theorem? (Again, assuming that $$\vec{A}_{\mu;r}$$ can be treated as a "normal" vector.)