usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
504 stars 148 forks source link

How can I solve this error 'all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)' in the following code? #934

Closed ghost closed 1 year ago

ghost commented 1 year ago
def Flow():
    #Import fipy
    from fipy import CellVariable, FaceVariable, Grid2D, DiffusionTerm, Viewer, TransientTerm, ImplicitSourceTerm
    ############################################################################################################################

    # Define the initial & physical parameters
    viscosity = 0.1
    diffusivity = 1.0
    interfacialAreaCoeff = 0.1
    u = 1.0
    v = 0.0
    k = 1.0
    ############################################################################################################################

    # Set up the domain
    L1 = 5.0    #cm
    L2 = 20    #cm
    x = 10.0
    y = 2.0
    dy = L1 / y
    dx = L2 / x
    mesh = Grid2D(nx=x, ny=y, dx=dx, dy=dy)
    x, y = mesh.cellCenters
    ############################################################################################################################

    # Define the fluid velocity
    velocity = FaceVariable(name="Velocity", mesh=mesh, value=(u,v))

    #Define the pressure
    pressure = CellVariable(name="Pressure", mesh=mesh)

    # Define the species concentration
    C = CellVariable(name="Concentration", mesh=mesh, hasOld=True)

    # Define the volume fraction
    phi = CellVariable(name="Volume Fraction", mesh=mesh, hasOld=True)

    # Define the interfacial area
    A = CellVariable(name="Interfacial Area", mesh=mesh, hasOld=True)
    ############################################################################################################################

    # Define the initial and boundary conditions
    C.setValue(0.0)
    C.constrain(1.0, mesh.facesBottom)
    phi.setValue(0.5)
    phi.constrain(0.5, mesh.facesRight)
    A.setValue(0.0)
    velocity.constrain(0.0, mesh.exteriorFaces)
    velocity.constrain(1.0, mesh.facesLeft)
    ############################################################################################################################

    # Define the chemical reaction
    reaction = -k * C

    # Define the transport coefficients
    fluidCoeff = FaceVariable(mesh=mesh, value=viscosity)
    speciesCoeff = FaceVariable(mesh=mesh, value=diffusivity)
    interfacialAreaCoeff = FaceVariable(mesh=mesh, value=interfacialAreaCoeff)
    ############################################################################################################################

    # Define the Navier-Stokes equations
    equ1 = TransientTerm(var=velocity) == DiffusionTerm(coeff=fluidCoeff, var=velocity) - pressure.grad
    equ2 = TransientTerm(var=C) == DiffusionTerm(coeff=speciesCoeff, var=C) + reaction
    equ3 = TransientTerm(var=phi) == -velocity.divergence
    equ4 = TransientTerm(var=A) == DiffusionTerm(coeff=interfacialAreaCoeff, var=phi) - A / phi * velocity.divergence
    equ5 = ImplicitSourceTerm(coeff=1./fluidCoeff.divergence, var=pressure) == velocity.divergence
    ############################################################################################################################

    # Solve the equations
    eqns = equ1 & equ2 & equ3 & equ4 & equ5
    timeStepDuration = 0.1
    steps = 10
    for step in range(steps):
        eqns.solve(dt=timeStepDuration)
        #DiffusionTerm().solve(eqns)
    ############################################################################################################################

    # Visualize the results
    viewer = Viewer(vars=(C, velocity), datamin=0.0, datamax=1.0)
    viewer.plot()
Flow() 

Acid_Model_9.zip

guyer commented 1 year ago

The error means that the variable velocity of eqn1 is rank-1 (vector) and the variable C of eqn2 is rank-0 (scalar). FiPy doesn't support this.

Further issues:

ghost commented 1 year ago

The error means that the variable velocity of eqn1 is rank-1 (vector) and the variable C of eqn2 is rank-0 (scalar). FiPy doesn't support this.

Further issues:

  • velocity is a FaceVariable and the solution variable of an equation must be a CellVariable
  • fluidCoeff is a scalar field, so fluidCoeff.divergence doesn't make sense

Would you mind correcting the code?

guyer commented 1 year ago

I don't know what you're trying to represent with ImplicitSourceTerm(coeff=1./fluidCoeff.divergence, var=pressure), so I couldn't guess what this term should be.

I have previously referred you to https://github.com/usnistgov/fipy/blob/master/examples/reactiveWetting/liquidVapor2D.py. That example treats velocity as it needs to be in 2D in FiPy. I strongly recommend you use that example as the foundation for your code.

ghost commented 1 year ago

I don't know what you're trying to represent with ImplicitSourceTerm(coeff=1./fluidCoeff.divergence, var=pressure), so I couldn't guess what this term should be.

I have previously referred you to https://github.com/usnistgov/fipy/blob/master/examples/reactiveWetting/liquidVapor2D.py. That example treats velocity as it needs to be in 2D in FiPy. I strongly recommend you use that example as the foundation for your code.

I've tried to write the continuity equrstion in Navier-Stokes, is it correct?