Dear all,
I am trying to solve with fipy the contaminant transport problem which
looks like:
Latex:
\frac{partial phi}{partial t}=\nabla(S\nabla phi)-\nabla(phi.v)
where S is the identity matrix times delta ie S=delta x I , where delta is
a given parameter.
and v=(v1,v2) a given vector
the analytical solution is:
phi(x,y,t)=(1/F)exp((-50*(x-x0-v1*t)**2+(y-y0-v2*t)**2) /F)
where F=(200*nu*t+1) and nu is a given parameter and x=(x0,y0) a given
point.
this analytical solution is Gaussian peak starting at the point x=(x0,y0),
being transported by the convective field v=(v1,v2), and diffusing.
we consider this domain: (0,3)x(0,3) and T=2
I write a program and I got some plots, but these plots don't give what I am
expecting. here is my program:
I get some problem when the diffusion coefficient is a matrix and not a
scalar and the analytical solution is a function of x,y and t
**************************************************************************************************************************************
gamma=0.1
nu=.1
convCoeff=v=(0.8, 0.4)
x=[0.5,1.35]
L=3
nx=30
ny=nx
dx=3./nx
dy=dx
from fipy.meshes.grid2D import Grid2D
from fipy import numerix
from scipy import *
from scipy.linalg import *
diffcoeff=S=gamma*numerix.array([[1, 0],[0 ,1]]) #my coefficient here is a
matrix, not a scalar as in all example in the documentation I don't know if
this a good way to set this matrix in fipy
mesh=Grid2D(dx=dx,dy=dy,nx=nx,ny=ny)
from fipy.variables.cellVariable import CellVariable
phi=CellVariable(name="solution variable", mesh=mesh) #initializing the
solution variable to zero
#here I am building my equation in fipy, is this correct?
from fipy.terms.transientTerm import TransientTerm
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.exponentialConvectionTerm import ExponentialConvectionTerm
diffTerm= ImplicitDiffusionTerm(coeff=diffcoeff)
eqX=TransientTerm()==diffTerm-ExponentialConvectionTerm(coeff=convCoeff,
diffusionTerm=diffTerm)
xx=mesh.getCellCenters()
D=norm(S)# because S is a matrix I take it matrix norm and affect it to D, I
don't know if this is correct, to take the norm of S, or maybe there another
way to deal with this matrix coefficient
xcoord=xx[...,0]
ycoord=xx[...,1]
phi.setValue(numerix.exp(-50*(xcoord-x[0])**2+(ycoord-x[1])**2))# here I am
setting the initial condition at t=0
# is this time step duration correct? where D is the matrix norm of S,
should I use the matrix norm of S?
timeStepDuration=0.9*dx**2/(2*D)
steps=100
t=timeStepDuration*steps
phiAnalytical=CellVariable(name="analytical value", mesh=mesh)
try:
phiAnalytical.setValue((1/(200*nu*t+1))*numerix.exp((-50*((xcoord-x[0]-v[0]*t)**2+(ycoord-x[1]-v[1]*t)**2))/(200*nu*t+1)))
except importError:
print "The Scipy library is not available to test the solution to the
transient diffusion equation"
if __name__ =='__main__':
from fipy import viewers
viewer=viewers.make(vars=(phi, phiAnalytical), limits={'datamin':0,
'datamax':1})
viewer.plot()
for step in range(steps):
eqX.solve(var=phi,dt=timeStepDuration)
if __name__ == '__main__':
viewer.plot()
print phi.allclose(phiAnalytical, atol = 7e-4)
Moreover, I need to use the final time T=2 in my code,
Thanks for your help in giving me some hint to deal with this problem
Regards
Franck
--
***********************************************
***********************************************
Franck Kalala Mutombo
[EMAIL PROTECTED]
African Institut for Mathematical Sciences, Muizenberg Cape Town, South
Africa
[EMAIL PROTECTED]
************************************************
************************************************
