Hi, I am trying to get numerical solution of multicomponent diffusion in one dimension in mineral. I have system of transient PDEs with transient non-linear coefficient, boundary and intial conditions (see image with general equation and information below). $$\frac{\partial}{\partial t}X_{i}=\frac{\partial}{\partial x}\sum_{j=1}^{n-1}D_{ij}\frac{\partial}{\partial x}X_{j}$$ where $D_{ij}$ is $$D_{ij}=D_{i}\delta_{ij}-\frac{D_{i}X_{i}Z_{j}Z_{i}(D_{j}-D_{n})}{\sum_{l=1}^nX_{l}Z_{l}^2D_{l}}$$ where Z is charge of element
There are four components with equal charge (XMg+XFe+XCa+XMn=1, n=4, and we consider Mn as dependent component), interval of diffusion (in meter, x) and time interval (in sec, t), intial concentration (at t=0, Piecewise function), assumption that concentrations of components at the ends of interval do not change (boundary condition) and self-diffusion coefficient ( $D_{i}$ or $D_{j}$).
So, XMg, XFe, and XCa depend on t and x, and $D_{ij}$ depends on XMg[t,x], XFe[t,x], and XCa[t,x]. Depend variables occur as in denominator as in numerator in $D_{ij}$, so coefficients of system of 3 transient PDE is transient and non-linear. I read that Wolfram Mathematica v 12.1 can solve such equations, so I tried to implement code for getting of numerical solutions.
In general, I got not bad results with using Method of Lines and TensorProductGrid (NDSolve chose it automatically). I checked it by "ElementMesh" (=None) and {state} = NDSolveProcessEquations (={NDSolve
StateData["<" 0. ">"]}). But I want to get result with Method of Lines and FEM or with purely FEM. Here I have a problem.
When I re-write equations in the Coefficient Form (with Div, Grad, Inactive), condition in DirichletCondition form and try to evalute NDSolve, Wolfram Mathematica send me error message (see in attached notebook) "NDSolveValue::femper"; "PDE parsing error of Div .... Inconsistent equation dimensions."
What is going wrong? How should I correct my code? If Method of Line with TensorProductGrid works, at least Method of Line with FiniteElement will work too.
P.S. I tried to prediscribe fuctions for two of three components (for exmaple, XCa[t,x] and XFe[t,x]) and NDSolve can get result for third component (hence, equations are correct), but it did not work with two or more dependent variables.
Attachments: