Bug introduced in 12.0.0 - Fixed in Version: 12.2.0
For example, consider the following piar of PDEs with two unknown functions $u(x,y)$ and $v(x,y)$: $$ \partial_x^2 u+\partial_y\left(v\partial_y u\right)=0\,,\quad\partial_x^2 v+\partial_y\left(u\partial_y v\right)=0\,, $$ as well as some suitable boundary conditions. We could try the following code
Eq1 = Inactive[Div][{{1, 0}, {0, v[x, y]}}.Inactive[Grad][u[x, y], {x, y}], {x, y}]==0;
Eq2 = Inactive[Div][{{1, 0}, {0, u[x, y]}}.Inactive[Grad][v[x, y], {x, y}], {x, y}]==0;
NDSolve[{Eq1, Eq2, DirichletCondition[u[x, y] == Cos[x + y], True], DirichletCondition[v[x, y] == Cos[x + y], True]}, {u[x, y], v[x, y]}, {x, y} \[Element] Disk[]]
Then we would have an error
NDSolve::femper: PDE parsing error of Div[{{1,0},{0,u}}.Grad$146922]. Inconsistent equation dimensions.
In contrast, if we try to solve another similar pair of PDEs $$ \partial_x^2 u+\partial_y\left(u\partial_y u\right)=0\,,\quad\partial_x^2 v+\partial_y\left(v\partial_y v\right)=0\,, $$ with similar codes, no error is reported and the codes run well.
More generally, as long as the coefficients of the second derivatives of one unknown function depend on other unknown functions, this error comes out.
What should I do to avoid this error?
InitializePDECoefficientsand the low level fem code. $\endgroup$