FEM derivative matrix construction

Question

Assuming $A$ is a $3 \times 3$ matrix (and a function of $x$, and $z$) and $K_i, \beta_i$, and $\alpha(T)$ are known parameters, I need to solve the following equation, (with an implied sum over $\mu, i$, and $j$),

$$K_1 \partial^2_j A_{\mu i} + K_{23} \partial_i (\partial_j A_{\mu j}) = 2\beta_1 Tr(A A^T) A^*_{\mu i} + 2\beta_2 Tr(A A^\dagger) A_{\mu i} + 2 \beta_3 [A A^T A^*]_{\mu i} + 2 \beta_4 [A A^\dagger A]_{\mu i} + 2 \beta_5 [A^* A^T A]_{\mu i} + \alpha(T) Tr(A A^\dagger)$$

This system of equations can be simplified with a known form of $A$. For example,

$$A = \begin{pmatrix} A_{xx} & 0 & 0 \\ 0 & A_{yy} & 0 \\ 0 & 0 & A_{zz} \end{pmatrix}$$

the equations reduce to (when normalized),

$$a_{||}'' = -\frac{1}{5} a_{||} ( 2 a_{||}^2 + a_\perp^2 ) + \frac{2}{5} a_{||} ( 2 |a_{||}|^2 + |a_\perp|^2 ) + \frac{2}{5} |a_{||}|^2 a_{||} - a_{||}$$

and

$$3 a_\perp'' = -\frac{1}{5} a_\perp^* ( 2 a_{||}^2 + a_\perp^2 ) + \frac{2}{5} a_\perp ( 2 |a_{||}|^2 + |a_\perp|^2 ) + \frac{2}{5} |a_\perp|^2 a_\perp - a_\perp$$

These appear rather simple, and I figured that Mathematica could solve them with NDSolve, but it hasn't worked. With more complicated forms of $A$ and more complicated regions, I'm looking to use Mathematica's FEM solver, but I just don't understand how to build the derivative matrix. I've looked at Wolfram's PDE solving guide and their FEM guide, but I only see coefficients for gradients, divergences, etc., while I will also need mixed derivatives ($\partial_i (\partial_j A_{\mu j})$). How do I build the mixed derivative matrix in InitializePDECoefficients? I have some working code in C++ that using the finite difference method, and I build the matrices there, but shouldn't there be a better way than inserting every element? Here it seems that the matrices were large, but Mathematica documentation here, they use a $2 \times 2$ identity matrix, but for the coefficient? I'm a little confused there.

Or, how can I get NDSolve or NDSolveValue to work? The boundary conditions seem to be causing problems in these functions. The BC's for my equations are, (for the first BC, I just assume that zMax (in code below) is close enough to infinity.)

$$\partial_z A_{\alpha j}|_{z=0} = \frac{1}{b}A_{\alpha j}|_{z=0}, \text{ for } b \in [0,\infty)$$

and

$$\lim_{z \to \infty} A_{\alpha j} = 1.$$

My code so far (using NDSolve first),

b = 1;
{xMin, xMax} = {-5, 5};
{zMin, zMax} = {0, 20};
sol = NDSolve[
  {D[p[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[p[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) p[x, z] + 2/5 Abs[p[x, z]]^2 p[x, z] - p[x, z],
  3 D[s[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[s[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) s[x, z] + 2/5 Abs[s[x, z]]^2 s[x, z] - s[x, z],
  (D[s[x, z], z] /. {z -> zMin}) == 1/b s[x, zMin],
  s[x, zMax] == 1,
  p[x, zMax] == 1,
  (D[p[x, z], z] /. {z -> zMin}) == 0 },
  {p, s}, {x, xMin, xMax}, {z, zMin, zMax}
]

This code gives the same error twice (I thought I had other errors previously, but I can't seem to recreate them).

NDSolve: The expression s^(0,1)[x, 0] == s given as a spatial boundary condition for the possibly automatically chosen finite element method should not have explicit derivatives of the dependent variables. NeumannValue should be used to specify spatial derivatives at the boundary.

And my other attempt,

Needs["NDSolve`FEM`"]
b = 1;
{xMin, xMax} = {-5, 5};
{zMin, zMax} = {0, 20};
nRegion = ToNumericalRegion[Rectangle[{xMin, zMin}, {xMax, zMax}]];

Γpi = NeumannValue[0, z == 0];
Γpf = DirichletCondition[p[x, z] == 1, z == zMax];
Γsi = NeumannValue[s[x, z]/b, z == 0];
Γsf = DirichletCondition[s[x, z] == 1, z == zMax];

eqP = -D[p[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[p[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) p[x, z] + 2/5 Abs[p[x, z]]^2 p[x, z] - p[x, z];

eqS = -3 D[s[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[s[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) s[x, z] + 2/5 Abs[s[x, z]]^2 s[x, z] - s[x, z];

sol = NDSolveValue[{{eqP == Γpi, Γpf}, {eqS == Γsi, Γsf}}, {p, s}, {x, z} ∈ nRegion]

gives these errors, twice:

DiscretizePDE: The FEMStiffnessElements operator failed.

FindRoot: The minimal damping factor of 1/10000 has been reached.

FindRoot: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

NDSolveValue: PDESolve could not find a solution.

Answer 1

I do not fully understand what you want to do but try this:

Needs["NDSolve`FEM`"]
b = 1;
{xMin, xMax} = {-5, 5};
{zMin, zMax} = {0, 20};
nRegion = ToNumericalRegion[Rectangle[{xMin, zMin}, {xMax, zMax}]];

\[CapitalGamma]pi = NeumannValue[0, z == 0];
\[CapitalGamma]pf = DirichletCondition[p[x, z] == 1, z == zMax];
\[CapitalGamma]si = NeumannValue[s[x, z]/b, z == 0];
\[CapitalGamma]sf = DirichletCondition[s[x, z] == 1, z == zMax];

eqP = -Inactive[Div][{{0, 0}, {0, 1}} . 
      Inactive[Grad][p[x, z], {x, z}], {x, 
      z}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[p[x, z]] + 
    2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) p[x, z] + 
    2/5 Abs[p[x, z]]^2 p[x, z] - p[x, z];

eqS = -3 Inactive[Div][{{0, 0}, {0, 1}} . 
      Inactive[Grad][s[x, z], {x, z}], {x, 
      z}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[s[x, z]] + 
    2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) s[x, z] + 
    2/5 Abs[s[x, z]]^2 s[x, z] - s[x, z];

sol = NDSolveValue[{{eqP == \[CapitalGamma]pi, \[CapitalGamma]pf}, \
{eqS == \[CapitalGamma]si, \[CapitalGamma]sf}}, {p, 
   s}, {x, z} \[Element] nRegion]

This returns two interpolating functions, not sure if correct.