How can I solve this error 'Cannot cast ufunc 'add' output from dtype('float64') to dtype('int32') with casting rule 'same_kind'' in the following code?

[ghost]

L1 = 1
L2 = 2
y = 10
x = 20
dy = L1 / y
dx = L2 / x
mesh = Grid2D(nx=x, ny=y, dx=dx, dy=dy)
###############################################################################################

molarWeight = 0.118
ee = -0.455971
gasConstant = 8.314
temperature = 650
vbar = 1.3e-05
liquidDensity = 7354.3402662299995
vaporDensity = 82.855803327810008
epsilon = 1e-16
###############################################################################################

def f(rho):
    return ee * rho**2 / molarWeight**2 + gasConstant * temperature * rho / molarWeight * numerix.log(rho / (molarWeight - vbar * rho))
def mu(rho):
    return 2 * ee * rho / molarWeight**2 + gasConstant * temperature / molarWeight * (numerix.log(rho / (molarWeight - vbar * rho)) + molarWeight / (molarWeight - vbar * rho))
def P(rho):
    return rho * mu(rho) - f(rho)
###############################################################################################

density = CellVariable(mesh=mesh, hasOld=True, name='density', value=1)
velocity = CellVariable(mesh=mesh, hasOld=True, name='velocity', value=1)
densityPrevious = density.copy()
velocityPrevious = velocity.copy()

potentialNC = CellVariable(mesh=mesh, name='potential')

freeEnergy = (f(density) + epsilon * temperature / 2 * density.grad.mag**2).cellVolumeAverage

matrixDiagonal = CellVariable(mesh=mesh, name='matrix', value=1e+20, hasOld=True)
correctionCoeff = mesh._faceAreas * mesh._cellDistances / matrixDiagonal.faceValue
massEqn = TransientTerm(var=density) + VanLeerConvectionTerm(coeff=velocity.faceValue + correctionCoeff * (density * potentialNC.grad).faceValue, var=density) - DiffusionTerm(coeff=correctionCoeff * density.faceValue**2, var=potentialNC)

viscosity = 1e-3
ConvectionTerm = CentralDifferenceConvectionTerm
momentumEqn = TransientTerm(coeff=density, var=velocity) + ConvectionTerm(coeff=[[1]] * density.faceValue * velocity.faceValue, var=velocity)== DiffusionTerm(coeff=2 * viscosity, var=velocity) - ConvectionTerm(coeff=density.faceValue * [[1]], var=potentialNC) + ImplicitSourceTerm(coeff=density.grad[0], var=potentialNC)
###############################################################################################

#boundary condition
velocity.constrain(0, mesh.exteriorFaces)
###############################################################################################

potentialDerivative = 2 * ee / molarWeight**2 + gasConstant * temperature * molarWeight / density / (molarWeight - vbar * density)**2

potential = mu(density)

potentialNCEqn = ImplicitSourceTerm(coeff=1, var=potentialNC) == potential + ImplicitSourceTerm(coeff=potentialDerivative, var=density) - potentialDerivative * density - DiffusionTerm(coeff=epsilon * temperature, var=density)

potentialNC.faceGrad.constrain(value=[0], where=mesh.exteriorFaces)

coupledEqn = massEqn & momentumEqn & potentialNCEqn

numerix.random.seed(2011)
density[:] = (liquidDensity + vaporDensity) / 2 * (1  + 0.01 * (2 * numerix.random.random(mesh.numberOfCells) - 1))
###############################################################################################

from fipy import input
if __name__ == '__main__':
    viewers = Viewer(density), Viewer(velocity), Viewer(potentialNC)
    for viewer in viewers:
        viewer.plot()
    input('Arrange viewers, then press <return> to proceed...')
    for viewer in viewers:
        viewer.plot()

cfl = 0.1
tolerance = 1e-1
dt = 1e-14
timestep = 0
relaxation = 0.5
if __name__ == '__main__':
    totalSteps = 1e2
else:
    totalSteps = 10

for timestep in range(int(totalSteps)):
    sweep = 0
    dt *= 1.1
    residual = 1
    initialResidual = None

    if timestep % 10 == 0:
        density.updateOld()
        velocity.updateOld()
        matrixDiagonal.updateOld()

while residual > tolerance:
    densityPrevious[:] = density
    velocityPrevious[:] = velocity
    previousResidual = residual

    dt = min(dt, dx / max(abs(velocity)) * cfl)

    coupledEqn.cacheMatrix()
    residual = coupledEqn.sweep(dt=dt)

    if initialResidual is None:
        initialResidual = residual

    residual = residual / initialResidual

    if residual > previousResidual * 1.1 or sweep > 20:
        density[:] = density.old
        velocity[:] = velocity.old
        matrixDiagonal[:] = matrixDiagonal.old
        dt /= 10
        if __name__ == '__main__':
            print('Recalculate the time step')
        timestep -= 1
        break
    else:
        matrixDiagonal[:] = coupledEqn.matrix.takeDiagonal()[mesh.numberOfCells:2 * mesh.numberOfCells]
        density[:] = (relaxation * density + (1 - relaxation) * densityPrevious).astype(density.dtype)
        velocity[:] = (relaxation * velocity + (1 - relaxation) * velocityPrevious)

    sweep += 1

    if __name__ == '__main__' and timestep % 10 == 0:
        print('timestep: %e / %e, dt: %1.5e, free energy: %1.5e' % (timestep, totalSteps, dt, freeEnergy))
        for viewer in viewers:
            viewer.plot()

    timestep += 1

[guyer]

As somebody asked you on Stack Overflow: "What line is the error on? And what is the type here at the bottom of your code?" I.e., timestep += 1type here is not valid Python. What are you trying to do?

[ghost]

As somebody asked you on Stack Overflow: "What line is the error on? And what is the type here at the bottom of your code?" I.e., timestep += 1type here is not valid Python. What are you trying to do?

on this--->coupledEqn.cacheMatrix()

[guyer]

Please provide a minimal reproducible example (which would include the full error traceback and a complete code that can actually be run). There is no way to figure out what's happening based on what you've provided.

There is good guidance in http://www.catb.org/~esr/faqs/smart-questions.html

[ghost]

Dear Jonathan Guyer, This is the whole file which I face the error in.

On Thu, May 11, 2023 at 7:18 PM Jonathan Guyer @.***> wrote:

Please provide a minimal reproducible example https://stackoverflow.com/help/minimal-reproducible-example (which would include the full error traceback and a complete code that can actually be run). There is no way to figure out what's happening based on what you've provided.

There is good guidance in http://www.catb.org/~esr/faqs/smart-questions.html

— Reply to this email directly, view it on GitHub https://github.com/usnistgov/fipy/issues/926#issuecomment-1544238170, or unsubscribe https://github.com/notifications/unsubscribe-auth/A4KJBT2U7CKLBW2OUATGFJLXFUC3BANCNFSM6AAAAAAX6IEGEY . You are receiving this because you authored the thread.Message ID: @.***>

[guyer]

It cannot be. Trying to run that code would produce at least four errors before ever getting to coupledEqn.cacheMatrix().

[ghost]

On Thu, May 11, 2023 at 11:30 PM Jonathan Guyer @.***> wrote:

It cannot be. Trying to run that code would produce at least four errors before ever getting to coupledEqn.cacheMatrix().

— Reply to this email directly, view it on GitHub https://github.com/usnistgov/fipy/issues/926#issuecomment-1544597512, or unsubscribe https://github.com/notifications/unsubscribe-auth/A4KJBT6X7JHAS2HAZUKUIQ3XFVAM5ANCNFSM6AAAAAAX6IEGEY . You are receiving this because you authored the thread.Message ID: @.***>

Can you help to solve this issue? I’ve tried everything.

[guyer]

If you don't share the entire code, then no, I can't help.

What is residual? What is tolerance? What is density? What is velocity? What is dt? What is dx? What is cfl? What is coupledEqn? And on... and on... and on...

[ghost]

On Fri, May 12, 2023 at 4:48 AM Jonathan Guyer @.***> wrote:

If you don't share the entire code, then no, I can't help.

What is residual? What is tolerance? What is density? What is velocity? What is dt? What is dx? What is cfl? What is coupledEqn? And on... and on... and on...

— Reply to this email directly, view it on GitHub https://github.com/usnistgov/fipy/issues/926#issuecomment-1544967660, or unsubscribe https://github.com/notifications/unsubscribe-auth/A4KJBT3BSPVTFKXCYXPN233XFWFVPANCNFSM6AAAAAAX6IEGEY . You are receiving this because you authored the thread.Message ID: @.***>

I shared the entire file yesterday but I think you haven’t received it yet. I’ll email it to you directly. Would you mind sending your email address?

[guyer]

I see you've posted the full code into your Stack Overflow question (please only ask in one place). I've copied it here:

L1 = 1
L2 = 2
y = 10
x = 20
dy = L1 / y
dx = L2 / x
mesh = Grid2D(nx=x, ny=y, dx=dx, dy=dy)
############################################################################

molarWeight = 0.118
ee = -0.455971
gasConstant = 8.314
temperature = 650
vbar = 1.3e-05
liquidDensity = 7354.3402662299995
vaporDensity = 82.855803327810008
epsilon = 1e-16
############################################################################

def f(rho):
    return ee * rho**2 / molarWeight**2 + gasConstant * temperature * rho / molarWeight * numerix.log(rho / (molarWeight - vbar * rho))
def mu(rho):
    return 2 * ee * rho / molarWeight**2 + gasConstant * temperature / molarWeight * (numerix.log(rho / (molarWeight - vbar * rho)) + molarWeight / (molarWeight - vbar * rho))
def P(rho):
    return rho * mu(rho) - f(rho)
############################################################################

density = CellVariable(mesh=mesh, hasOld=True, name='density', value=1)
velocity = CellVariable(mesh=mesh, hasOld=True, name='velocity', value=1)
densityPrevious = density.copy()
velocityPrevious = velocity.copy()

potentialNC = CellVariable(mesh=mesh, name='potential')

freeEnergy = (f(density) + epsilon * temperature / 2 * density.grad.mag**2).cellVolumeAverage

matrixDiagonal = CellVariable(mesh=mesh, name='matrix', value=1e+20, hasOld=True)
correctionCoeff = mesh._faceAreas * mesh._cellDistances / matrixDiagonal.faceValue
massEqn = TransientTerm(var=density) + VanLeerConvectionTerm(coeff=velocity.faceValue + correctionCoeff * (density * potentialNC.grad).faceValue, var=density) - DiffusionTerm(coeff=correctionCoeff * density.faceValue**2, var=potentialNC)

viscosity = 1e-3
ConvectionTerm = CentralDifferenceConvectionTerm
momentumEqn = TransientTerm(coeff=density, var=velocity) + ConvectionTerm(coeff=[[1]] * density.faceValue * velocity.faceValue, var=velocity)== DiffusionTerm(coeff=2 * viscosity, var=velocity) - ConvectionTerm(coeff=density.faceValue * [[1]], var=potentialNC) + ImplicitSourceTerm(coeff=density.grad[0], var=potentialNC)
############################################################################

#boundary condition
velocity.constrain(0, mesh.exteriorFaces)
############################################################################

potentialDerivative = 2 * ee / molarWeight**2 + gasConstant * temperature * molarWeight / density / (molarWeight - vbar * density)**2

potential = mu(density)

potentialNCEqn = ImplicitSourceTerm(coeff=1, var=potentialNC) == potential + ImplicitSourceTerm(coeff=potentialDerivative, var=density) - potentialDerivative * density - DiffusionTerm(coeff=epsilon * temperature, var=density)

potentialNC.faceGrad.constrain(value=[0], where=mesh.exteriorFaces)

coupledEqn = massEqn & momentumEqn & potentialNCEqn

numerix.random.seed(2011)
density[:] = (liquidDensity + vaporDensity) / 2 * (1  + 0.01 * (2 * numerix.random.random(mesh.numberOfCells) - 1))
############################################################################

from fipy import input
if __name__ == '__main__':
    viewers = Viewer(density), Viewer(velocity), Viewer(potentialNC)
    for viewer in viewers:
        viewer.plot()
    input('Arrange viewers, then press <return> to proceed...')
    for viewer in viewers:
        viewer.plot()

cfl = 0.1
tolerance = 1e-1
dt = 1e-14
timestep = 0
relaxation = 0.5
if __name__ == '__main__':
    totalSteps = 1e2
else:
    totalSteps = 10

for timestep in range(int(totalSteps)):
    sweep = 0
    dt *= 1.1
    residual = 1
    initialResidual = None

    if timestep % 10 == 0:
        density.updateOld()
        velocity.updateOld()
        matrixDiagonal.updateOld()

while residual > tolerance:
    densityPrevious[:] = density
    velocityPrevious[:] = velocity
    previousResidual = residual

    dt = min(dt, dx / max(abs(velocity)) * cfl)

    coupledEqn.cacheMatrix()
    residual = coupledEqn.sweep(dt=dt)

    if initialResidual is None:
        initialResidual = residual

    residual = residual / initialResidual

    if residual > previousResidual * 1.1 or sweep > 20:
        density[:] = density.old
        velocity[:] = velocity.old
        matrixDiagonal[:] = matrixDiagonal.old
        dt /= 10
        if __name__ == '__main__':
            print('Recalculate the time step')
        timestep -= 1
        break
    else:
        matrixDiagonal[:] = coupledEqn.matrix.takeDiagonal()[mesh.numberOfCells:2 * mesh.numberOfCells]
        density[:] = (relaxation * density + (1 - relaxation) * densityPrevious).astype(density.dtype)
        velocity[:] = (relaxation * velocity + (1 - relaxation) * velocityPrevious)

    sweep += 1

    if __name__ == '__main__' and timestep % 10 == 0:
        print('timestep: %e / %e, dt: %1.5e, free energy: %1.5e' % (timestep, totalSteps, dt, freeEnergy))
        for viewer in viewers:
            viewer.plot()

    timestep += 1

There are a couple of issues:

• Solution variables must be floats, not integers. Change

  - density = CellVariable(mesh=mesh, hasOld=True, name='density', value=1)
  - velocity = CellVariable(mesh=mesh, hasOld=True, name='velocity', value=1)
  + density = CellVariable(mesh=mesh, hasOld=True, name='density', value=1.)
  + velocity = CellVariable(mesh=mesh, hasOld=True, name='velocity', value=1.)

• It doesn't make sense to try to cast the type of a FiPy expression (and it's not in the example you appear to have gotten this from). Change

  - density[:] = (relaxation * density + (1 - relaxation) * densityPrevious).astype(density.dtype)
  + density[:] = (relaxation * density + (1 - relaxation) * densityPrevious)

[ghost]

I see you've posted the full code into your Stack Overflow question (please only ask in one place). I've copied it here:

L1 = 1

L2 = 2

y = 10

x = 20

dy = L1 / y

dx = L2 / x

mesh = Grid2D(nx=x, ny=y, dx=dx, dy=dy)

############################################################################

molarWeight = 0.118

ee = -0.455971

gasConstant = 8.314

temperature = 650

vbar = 1.3e-05

liquidDensity = 7354.3402662299995

vaporDensity = 82.855803327810008

epsilon = 1e-16

############################################################################

def f(rho):

    return ee * rho**2 / molarWeight**2 + gasConstant * temperature * rho / molarWeight * numerix.log(rho / (molarWeight - vbar * rho))

def mu(rho):

    return 2 * ee * rho / molarWeight**2 + gasConstant * temperature / molarWeight * (numerix.log(rho / (molarWeight - vbar * rho)) + molarWeight / (molarWeight - vbar * rho))

def P(rho):

    return rho * mu(rho) - f(rho)

############################################################################

density = CellVariable(mesh=mesh, hasOld=True, name='density', value=1)

velocity = CellVariable(mesh=mesh, hasOld=True, name='velocity', value=1)

densityPrevious = density.copy()

velocityPrevious = velocity.copy()

potentialNC = CellVariable(mesh=mesh, name='potential')

freeEnergy = (f(density) + epsilon * temperature / 2 * density.grad.mag**2).cellVolumeAverage

matrixDiagonal = CellVariable(mesh=mesh, name='matrix', value=1e+20, hasOld=True)

correctionCoeff = mesh._faceAreas * mesh._cellDistances / matrixDiagonal.faceValue

massEqn = TransientTerm(var=density) + VanLeerConvectionTerm(coeff=velocity.faceValue + correctionCoeff * (density * potentialNC.grad).faceValue, var=density) - DiffusionTerm(coeff=correctionCoeff * density.faceValue**2, var=potentialNC)

viscosity = 1e-3

ConvectionTerm = CentralDifferenceConvectionTerm

momentumEqn = TransientTerm(coeff=density, var=velocity) + ConvectionTerm(coeff=[[1]] * density.faceValue * velocity.faceValue, var=velocity)== DiffusionTerm(coeff=2 * viscosity, var=velocity) - ConvectionTerm(coeff=density.faceValue * [[1]], var=potentialNC) + ImplicitSourceTerm(coeff=density.grad[0], var=potentialNC)

############################################################################

#boundary condition

velocity.constrain(0, mesh.exteriorFaces)

############################################################################

potentialDerivative = 2 * ee / molarWeight**2 + gasConstant * temperature * molarWeight / density / (molarWeight - vbar * density)**2

potential = mu(density)

potentialNCEqn = ImplicitSourceTerm(coeff=1, var=potentialNC) == potential + ImplicitSourceTerm(coeff=potentialDerivative, var=density) - potentialDerivative * density - DiffusionTerm(coeff=epsilon * temperature, var=density)

potentialNC.faceGrad.constrain(value=[0], where=mesh.exteriorFaces)

coupledEqn = massEqn & momentumEqn & potentialNCEqn

numerix.random.seed(2011)

density[:] = (liquidDensity + vaporDensity) / 2 * (1  + 0.01 * (2 * numerix.random.random(mesh.numberOfCells) - 1))

############################################################################

from fipy import input

if __name__ == '__main__':

    viewers = Viewer(density), Viewer(velocity), Viewer(potentialNC)

    for viewer in viewers:

        viewer.plot()

    input('Arrange viewers, then press <return> to proceed...')

    for viewer in viewers:

        viewer.plot()

cfl = 0.1

tolerance = 1e-1

dt = 1e-14

timestep = 0

relaxation = 0.5

if __name__ == '__main__':

    totalSteps = 1e2

else:

    totalSteps = 10

for timestep in range(int(totalSteps)):

    sweep = 0

    dt *= 1.1

    residual = 1

    initialResidual = None

    if timestep % 10 == 0:

        density.updateOld()

        velocity.updateOld()

        matrixDiagonal.updateOld()

while residual > tolerance:

    densityPrevious[:] = density

    velocityPrevious[:] = velocity

    previousResidual = residual

    dt = min(dt, dx / max(abs(velocity)) * cfl)

    coupledEqn.cacheMatrix()

    residual = coupledEqn.sweep(dt=dt)

    if initialResidual is None:

        initialResidual = residual

    residual = residual / initialResidual

    if residual > previousResidual * 1.1 or sweep > 20:

        density[:] = density.old

        velocity[:] = velocity.old

        matrixDiagonal[:] = matrixDiagonal.old

        dt /= 10

        if __name__ == '__main__':

            print('Recalculate the time step')

        timestep -= 1

        break

    else:

        matrixDiagonal[:] = coupledEqn.matrix.takeDiagonal()[mesh.numberOfCells:2 * mesh.numberOfCells]

        density[:] = (relaxation * density + (1 - relaxation) * densityPrevious).astype(density.dtype)

        velocity[:] = (relaxation * velocity + (1 - relaxation) * velocityPrevious)

    sweep += 1

    if __name__ == '__main__' and timestep % 10 == 0:

        print('timestep: %e / %e, dt: %1.5e, free energy: %1.5e' % (timestep, totalSteps, dt, freeEnergy))

        for viewer in viewers:

            viewer.plot()

    timestep += 1

There are a couple of issues:

• Solution variables must be floats, not integers. Change

  
  - density = CellVariable(mesh=mesh, hasOld=True, name='density', value=1)
  
  - velocity = CellVariable(mesh=mesh, hasOld=True, name='velocity', value=1)
  
  + density = CellVariable(mesh=mesh, hasOld=True, name='density', value=1.)
  
  + velocity = CellVariable(mesh=mesh, hasOld=True, name='velocity', value=1.)
  

• It doesn't make sense to try to cast the type of a FiPy expression (and it's not in the example you appear to have gotten this from). Change

  
  - density[:] = (relaxation * density + (1 - relaxation) * densityPrevious).astype(density.dtype)
  
  + density[:] = (relaxation * density + (1 - relaxation) * densityPrevious)
  

Thank you🙏🏼