Hi Andrew,
Can you attach './HalfCirc.msh' so I can run the script.
Cheers
On Jan 22, 2008 4:43 PM, asreeve <[EMAIL PROTECTED]> wrote:
>
> Daniel,
>
> vel is a 3 by number of faces array... x,y,z components of velocity at
> each face. In the array I'm using, vel[0] and vel[1] are all nearly zero,
> and vel[2] is + or - 2.5e-4. The magnitude of these numbers is
> correct.
>
> passing 'vel' to the conductance term doesn't work. I get an error:
>
>
> <type 'exceptions.TypeError'>: The coefficient must be a vector value.
> WARNING: Failure executing file: <PeatCylSimple.py>
>
>
> Both of the other options don't result in an error, just a wrong answer.
>
> Should the coeff in the convection term be a vector normal to the face?
> As an aside, the equation I'm solving is the standard Advection-diffusion
> equation:
>
> dC/dx=d/dX(D dC/dX) - d(v*C)/dX
>
> v=velocity
> C=concentration
> D=diffusion coef
>
> I get a sensible answer only when I switch the sign on the diffusion term
> in my transport equation.
>
>
>
>
>
> Here's the script I'm using (running anisotropy Branch of FiPy). The last
> couple of lines calculate the analytic solution.:
>
> from scipy.special import erfc
> def addLayer(verts):
> addArray = verts.copy()
> addArray[:] = 0
> addArray[2] = length/slices
> return verts + addArray
>
> ###create mesh, hemi-cylinder of peat, 40 cm long, 2.54 cm radius
> length=40. ##cm
> slices=40
> mesh=GmshImporter2D('./HalfCirc.msh').extrude(layers=slices,
> extrudeFunc=addLayer)
>
> ###calculate hydraulic head in peat
> head=CellVariable(name='hydraulic head',
> mesh=mesh,value=100.,hasOld=1) ## in cm
> Dist=random.uniform(1.,1.,(mesh.getNumberOfCells()))
> Kxx=1.e-6*100*Dist ##hyd. cond. in cm/sec
> Ss=1.e-2/100. ##Sp.Storage in 1/cm
> kond=CellVariable(name='hydraulic conductivity',mesh=mesh,rank=0,value=Kxx)
> spStore=CellVariable(name='spec. storage',mesh=mesh,rank=0,value=Ss)
> effPor=CellVariable(name='effective porosity',mesh=mesh,rank=0,value=0.2)
>
> gwFlowEq=TransientTerm(spStore)==DiffusionTerm(kond.getHarmonicFaceValue())
>
> ##hydraulic gradients across peat are 0.5 or 0.8, units in cm, 40 cm long tube
> X,Y,Z=mesh.getFaceCenters()
>
> ZMax=max(Z)-0.05
> bndry=(FixedValue(mesh.getFaces().where((Z<0.05)),120.),
> FixedValue(mesh.getFaces().where((Z>ZMax)),100.))
> delTime=1.e12
> timeSteps=3
> for step in range(timeSteps):
> head.updateOld()
> resid=1.
> while resid>1.e-12:
> resid=gwFlowEq.sweep(var=head,
> solver=LinearPCGSolver(iterations=500,
> tolerance=1.e-13),
> boundaryConditions=bndry,
> dt=delTime)
> print 'hyd.head',resid
>
> discharge=dot(head.getFaceGrad(),kond.getHarmonicFaceValue())
> vel=discharge/effPor.getArithmeticFaceValue()
> conc=CellVariable(name='solute conc',mesh=mesh,
> value=1.,hasOld=1)
> D=.0002
> Disp=FaceVariable(name='dispersion tensor',
> mesh=mesh,rank=0)
> Disp.setValue(D)
>
> X,Y,Z=mesh.getFaceCenters()
> ZMax=max(Z)-0.05
> bndry2=(FixedValue(mesh.getFaces().where((Z<0.05)),101.),
> FixedValue(mesh.getFaces().where((Z>ZMax)),1.))
>
> transportEq2=(TransientTerm()==
> -PowerLawConvectionTerm(coeff=vel())
> -DiffusionTerm(coeff=Disp))
>
> t=0
> delTime=1.e3
> timeSteps=40
> for step in range(timeSteps):
> t+=delTime
> conc.updateOld()
> resid=1.
> while resid>1.e-12:
> #for i in range(10):
> resid=transportEq2.sweep(var=conc,
> solver=LinearGMRESSolver(iterations=500,
> tolerance=1.e-15),
> boundaryConditions=bndry2,
> dt=delTime)
> print 'conc',resid
>
> X,Y,Z=mesh.getCellCenters()
> v=2.5e-4
> Dis=D
> OB=[1+100*.5*(erfc((x-v*t)/(2*sqrt(Dis*t)))+exp(v*x/Dis)*erfc((x+v*t)/(2*sqrt(Dis*t))))
> for x in Z]
> ----------------------
> Andrew Reeve
> Dept. of Earth Science
> University of Maine
> Orono, ME 04469
>
> > On Jan 22, 2008 1:41 PM, asreeve <[EMAIL PROTECTED]> wrote:
> >>
> >> Hello,
> >
> > Hi, hope all is well.
> >
> >> I've having some difficulty figuring out how to use the convection terms
> >> in a 3-D problem.
> >
> > I haven't solved a lot of 3D convection problems so I am not sure how well
> > tested everything is. Should work though.
> >
> >> I've calculated a velocity term for a 3-D simulations:
> >
> >> vel=dot(head.getFaceGrad(),kond.getHarmonicFaceValue())
> >
> > I'm confused about which axis is getting dotted. What are the shapes of
> > "kond"
> > and "vel"?
> >
> >> and would like to pass these velocities to a convection term.
> >
> > I assume that "vel" is rank 1. Is it?
> >
> >> I've simplified my problem to the point that there are non-zero velocities
> >> only in the z direction (but I need to solve a problem where vel varies in
> >> space and direction).
> >
> > Shouldn't be a problem.
> >
> >>
> >> I've tried using :
> >>
> >> PowerLawConvectionTerm(coeff=vel())
> >>
> >> and
> >>
> >> PowerLawConvectionTerm(coeff=(vel()[0],vel()[1],vel()[2]))
> >>
> >> in the equation term, but both of these yield an incorrect result.
> >>
> >> However, when I assign the following (values taken from vel variable),
> >>
> >> PowerLawConvectionTerm(coeff=(0,0,2.5e-4))
> >>
> >> the output is correct.
> >>
> >> What's the proper syntax to assign the calculated 'vel' to a convection
> >> term?
> >
> > Have you tried?
> >
> > PowerLawConvectionTerm(coeff=vel)
> >
> > Cheers
> >
> > --
> > Daniel Wheeler
> >
>
>
--
Daniel Wheeler
