ExplicitDiffusionTerm doesn't work in coupled equations

[eugenio-pino]

With the script attached I run into a weird error caused by ExplicitDiffusionTerm. If I remove those terms from the equations, the script runs correctly.

Traceback (most recent call last):
  File "double_ld.py", line 155, in <module>
    linear_decay_gi_equation.solve(dt = dt)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/term.py", line 176, in solve
    solver = self._prepareLinearSystem(var, solver, boundaryConditions, dt)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/term.py", line 130, in _prepareLinearSystem
    buildExplicitIfOther=self._buildExplcitIfOther)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/coupledBinaryTerm.py", line 88, in _buildAndAddMatrices
    buildExplicitIfOther=buildExplicitIfOther)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/binaryTerm.py", line 34, in _buildAndAddMatrices
    buildExplicitIfOther=buildExplicitIfOther)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/binaryTerm.py", line 34, in _buildAndAddMatrices
    buildExplicitIfOther=buildExplicitIfOther)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/binaryTerm.py", line 34, in _buildAndAddMatrices
    buildExplicitIfOther=buildExplicitIfOther)
  [Previous line repeated 1 more time]
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/unaryTerm.py", line 67, in _buildAndAddMatrices
    diffusionGeomCoeff=diffusionGeomCoeff)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/terms/explicitDiffusionTerm.py", line 34, in _buildMatrix
    return (var, SparseMatrix(mesh=var.mesh), b - L * var.value)
  File "/home/eugenio/miniconda3/envs/reactive_transport/lib/python3.6/site-packages/fipy/matrices/petscMatrix.py", line 726, in __mul__
    x = other[self._localNonOverlappingColIDs]
IndexError: index 1000 is out of bounds for axis 0 with size 1000

I'm trying to solve a system of two coupled equations with a Crank-Nicolson scheme, that's why that term is present in equations. What could the error be caused by?

https://pastebin.pl/view/12a21ecd

[guyer]

This appears to be because ExplicitDiffusionTerm doesn't work with coupled equations (at least not with these coupled equations).

We'll try to figure out why this doesn't work but, in the meantime, you can substitute, e.g.,

    +(0.5 * D * u_GI.faceGrad).divergence

for

    +ExplicitDiffusionTerm(var = u_GI, coeff = 0.5 * D)

There really isn't any need for ExplicitDiffusionTerm anymore, even if it worked right.

Alternatively, you could decouple your equations and solve them iteratively. They're barely coupled as it is.

[eugenio-pino]

Alternatively, you could decouple your equations and solve them iteratively. They're barely coupled as it is. Yes these are pretty simple but I've got some other non-linear ones which instead need this term necessarily.

+(0.5 * D * u_GI.faceGrad).divergence Thank you for your time! I will write down the formulation next time.

[guyer]

There's still a bug that needs fixing