FiPyexamples.reactiveWetting.liquidVapor2D¶
A 2D version of the 1D example.
>>> molarWeight = 0.118
>>> ee = -0.455971
>>> gasConstant = 8.314
>>> temperature = 650.
>>> vbar = 1.3e-05
>>> liquidDensity = 7354.3402662299995
>>> vaporDensity = 82.855803327810008
>>> from fipy import CellVariable, Grid2D, TransientTerm, VanLeerConvectionTerm, DiffusionTerm, ImplicitSourceTerm, ConvectionTerm, CentralDifferenceConvectionTerm, Viewer
>>> from fipy.tools import numerix
>>> 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))
>>> Lx = 1e-6
>>> nx = 100
>>> dx = Lx / nx
>>> mesh = Grid2D(nx=nx, ny=nx, dx=dx, dy=dx)
>>> density = CellVariable(mesh=mesh, hasOld=True, name=r'$\rho$')
>>> velocityX = CellVariable(mesh=mesh, hasOld=True, name=r'$u_x$')
>>> velocityY = CellVariable(mesh=mesh, hasOld=True, name=r'$u_y$')
>>> velocityVector = CellVariable(mesh=mesh, name=r'$\vec{u}$', rank=1)
>>> densityPrevious = density.copy()
>>> velocityXPrevious = velocityX.copy()
>>> velocityYPrevious = velocityY.copy()
>>> potentialNC = CellVariable(mesh=mesh, name=r'$\mu^{NC}$')
>>> epsilon = 1e-16
>>> freeEnergy = (f(density) + epsilon * temperature / 2 * density.grad.mag**2).cellVolumeAverage
>>> matrixDiagonal = CellVariable(mesh=mesh, name=r'$a_f$', value=1e+20, hasOld=True)
>>> correctionCoeff = mesh._faceAreas * mesh._cellDistances / matrixDiagonal.faceValue
>>> massEqn = TransientTerm(var=density) \
... + VanLeerConvectionTerm(coeff=velocityVector.faceValue + correctionCoeff \
... * (density * potentialNC.grad).faceValue, \
... var=density) \
... - DiffusionTerm(coeff=correctionCoeff * density.faceValue**2, var=potentialNC)
>>> viscosity = 1e-3
>>> ConvectionTerm = CentralDifferenceConvectionTerm
>>> ##matXX = numerix.array([[[2], [0]], [[0], [1]]])
>>> ##matYY = numerix.array([[[1], [0]], [[0], [2]]])
>>> ##matXY = numerix.array([[[0], [0.5]], [[0.5], [0]]])
>>> ##matYX = matXY
>>> matXX = 1
>>> matYY = 1
>>> matXY = 0
>>> matYX = 0
>>> momentumXEqn = TransientTerm(coeff=density, var=velocityX) \
... + ConvectionTerm(coeff=density.faceValue * velocityVector.faceValue, var=velocityX) \
... == DiffusionTerm(coeff=(viscosity * matXX,), var=velocityX) \
... + DiffusionTerm(coeff=(viscosity * matXY,), var=velocityY) \
... - ConvectionTerm(coeff=density.faceValue * [[1], [0]], var=potentialNC) \
... + ImplicitSourceTerm(coeff=density.grad[0], var=potentialNC)
>>> momentumYEqn = TransientTerm(coeff=density, var=velocityY) \
... + ConvectionTerm(coeff=density.faceValue * velocityVector.faceValue, var=velocityY) \
... == DiffusionTerm(coeff=(viscosity * matYY,), var=velocityY) \
... + DiffusionTerm(coeff=(viscosity * matYX,), var=velocityX) \
... - ConvectionTerm(coeff=density.faceValue * [[0], [1]], var=potentialNC) \
... + ImplicitSourceTerm(coeff=density.grad[1], var=potentialNC)
>>> velocityX.constrain(0, mesh.facesLeft & mesh.facesRight)
>>> velocityY.constrain(0, mesh.facesTop & mesh.facesBottom)
>>> 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], [0]], where=mesh.exteriorFaces)
>>> coupledEqn = massEqn & momentumXEqn & momentumYEqn & potentialNCEqn
>>> numerix.random.seed(2012)
>>> 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(velocityVector), Viewer(potentialNC)
... for viewer in viewers:
... viewer.plot()
... input('arrange viewers')
... for viewer in viewers:
... viewer.plot()
>>> cfl = 0.1
>>> tolerance = 1e-1
>>> dt = 1e-14
>>> timestep = 0
>>> relaxation = 0.5
>>> sweeps = 0
>>> if __name__ == '__main__':
... totalSteps = 1e+10
... totalSweeps = 1e+10
... else:
... totalSteps = 1
... totalSweeps = 1
Beginning with FiPy :error:`4`, solver Convergence tolerance is normalized by the magnitude of the right-hand-side (RHS) vector. For this particular problem, the initial residual is much smaller than the RHS and so the solver gets “stuck”. Changing the normalization to use the initial residual at the beginning of each sweep allows the solution to progress.
>>> from fipy.solvers import solver_suite
>>> if solver_suite == "petsc":
... # PETSc's default ILU preconditioner does not behave well
... # for this problem
... from fipy import SSORPreconditioner
... precon = SSORPreconditioner()
... else:
... precon = "default"
>>> solver = coupledEqn.getDefaultSolver(criterion="initial", precon=precon)
Note
PETSc intrinsically wants to use “preconditioned” normalization, which prevents the solver from getting “stuck” because the preconditioner effectively scales the RHS to be similar in magnitude to the residual. Unfortunately, this normalization is not available for the other solver suites, so we don’t use it as the default.
>>> while timestep < totalSteps:
...
... sweep = 0
... dt *= 1.1
... residual = 1.
... initialResidual = None
...
... density.updateOld()
... velocityX.updateOld()
... velocityY.updateOld()
... matrixDiagonal.updateOld()
...
... while residual > tolerance and sweeps < totalSweeps:
... sweeps += 1
... densityPrevious[:] = density
... velocityXPrevious[:] = velocityX
... velocityYPrevious[:] = velocityY
... previousResidual = residual
... velocityVector[0] = velocityX
... velocityVector[1] = velocityY
...
... dt = min(dt, dx / max(abs(velocityVector.mag)) * cfl)
...
... coupledEqn.cacheMatrix()
... residual = coupledEqn.sweep(dt=dt, solver=solver)
...
... if initialResidual is None:
... initialResidual = residual
...
... residual = residual / initialResidual
...
... if residual > previousResidual * 1.1 or sweep > 20:
... density[:] = density.old
... velocityX[:] = velocityX.old
... velocityY[:] = velocityY.old
... matrixDiagonal[:] = matrixDiagonal.old
... dt = 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
... velocityX[:] = relaxation * velocityX + (1 - relaxation) * velocityXPrevious
... velocityY[:] = relaxation * velocityY + (1 - relaxation) * velocityYPrevious
...
... sweep += 1
...
... if __name__ == '__main__' and timestep % 1 == 0:
... print('timestep: %e / %e, dt: %1.5e, free energy: %1.5e' % (timestep, totalSteps, dt, freeEnergy))
... for viewer in viewers:
... viewer.plot()
...
... timestep += 1
>>> from fipy import input
>>> if __name__ == '__main__':
... input('finished')
Last updated on Jun 26, 2024.
Created using Sphinx 7.1.2.