Thanks a lot Jon.
I appreciate your help.
I will try to apply your remarks but I will certainly still need your help
again ;-) to develop the code.
Thanks again.
Fabien
> Le 15 août 2018 à 19:46, Guyer, Jonathan E. Dr. (Fed)
> a écrit :
>
> Fabien -
>
> I think the code below should get you going.
>
> The changes I made were:
>
> - `xVelocity` and `zVelocity` changed to rank-0 CellVariables. FiPy *always*
> solves at cell centers.
> - Created a rank-1 FaceVariable to hold the velocity. The cell components
> must be interpolated and manually inserted into this field at each sweep.
> - Changed the sign of the DiffusionTerm in the heat equation; it's
> unconditionally unstable otherwise
> - Slowed down the time steps to better see the evolution
>
> The result is not identical to the image you showed, but I think that's down
> to boundary conditions, sign of gravity, etc.
>
> - Jon
>
> # -*- coding: utf-8 -*-
> from fipy import *
> import matplotlib.pyplot as plt
> # Parameter
> L = 1.
> N = 50.
> dL = L/N
> alpha = 0.0002# Thermical dilatation
> landa = 0.6 # Thermical conductivity
> ro0 = 1023. # Average volumic Mass
> g = 10. # Gravity
> # Mesh
> mesh = Grid2D(nx=N, ny=N, dx=dL, dy=dL)
> # Variables
> dT = CellVariable(name = 'dT', mesh = mesh, value = 0.)
> xVelocity = CellVariable(mesh=mesh, name='Xvelocity', value = 0.)
> zVelocity = CellVariable(mesh=mesh, name='Zvelocity', value = 0.)
> velocity = FaceVariable(mesh=mesh, name='velocity', rank=1)
> # Init Condition
> x = mesh.cellCenters[0]
> dT.setValue(1.)
> dT.setValue(-1., where = x > L/2)
> # Viewer
> viewer = None
> if __name__ == '__main__':
> viewer = Viewer(vars=dT, datamin=-1., datamax=1.)
> viewer.plotMesh()
> raw_input("...")
> # Boussinesq equations
> D = landa/ro0
> eqX = (TransientTerm(var=xVelocity) + ConvectionTerm (var=xVelocity,
> coeff=velocity) == 0)
> eqZ = (TransientTerm(var=zVelocity) + ConvectionTerm (var=zVelocity,
> coeff=velocity) + alpha*g*dT == 0)
> eqT = (TransientTerm(var=dT) + ConvectionTerm (var=dT, coeff=velocity) ==
> DiffusionTerm(var=dT, coeff=D))
> eq = eqX & eqZ & eqT
> # Solving Boussinesq equations
> timeStepDuration = 1 * dL**2 / (2 * D)
> steps = 50
> sweeps = 5
> for step in range(steps):
>for sweep in range(sweeps):
>eq.sweep(dt=timeStepDuration)
>velocity[0] = xVelocity.arithmeticFaceValue
>velocity[1] = zVelocity.arithmeticFaceValue
>velocity[..., mesh.exteriorFaces.value] = 0.
> if viewer is not None:
> viewer.plot()
> plt.pause(0.1)
> raw_input("… ")
>
>
>> On Aug 14, 2018, at 11:22 AM, fgendr01 wrote:
>>
>> Hello,
>>
>> I’m a PhD Student in La Rochelle University (France) and I’m working on
>> Boussinesq Approximation.
>> I’m starting with fipy and I’d like to make a little algorithm to solve
>> these equations in 2D (x,z) :
>>
>> In my problem T = T0 + dT.
>> (ux, uz) is the velocity.
>>
>> So, here is my algorithm :
>> # -*- coding: utf-8 -*-
>> from fipy import *
>> import matplotlib.pyplot as plt
>> # Parameter
>> L = 1.
>> N = 50.
>> dL = L/N
>> alpha = 0.0002 # Thermical dilatation
>> landa = 0.6 # Thermical conductivity
>> ro0 = 1023. # Average volumic Mass
>> g = 10. # Gravity
>> # Mesh
>> mesh = Grid2D(nx=N, ny=N, dx=dL, dy=dL)
>> # Variables
>> dT = CellVariable(name = 'dT', mesh = mesh, value = 0.)
>> xVelocity = FaceVariable(mesh=mesh, name='Xvelocity', value = 0., rank=1)
>> zVelocity = FaceVariable(mesh=mesh, name='Zvelocity', value = 0., rank=1)
>> # Init Condition
>> x = mesh.cellCenters[0]
>> dT.setValue(1.)
>> dT.setValue(-1., where = x > L/2)
>> # Viewer
>> viewer = None
>> if __name__ == '__main__':
>> viewer = Viewer(vars=dT, datamin=-1., datamax=1.)
>> viewer.plotMesh()
>> raw_input("...")
>> # Boussinesq equations
>> D = landa/ro0
>> velocity = ((xVelocity), (zVelocity))
>> eqX = (TransientTerm(var=xVelocity) + ConvectionTerm (var=xVelocity,
>> coeff=velocity) == 0)
>> eqZ = (TransientTerm(var=zVelocity) + ConvectionTerm (var=zVelocity,
>> coeff=velocity) + alpha*g*dT == 0)
>> eqT = (TransientTerm(var=dT) + ConvectionTerm (var=dT, coeff=velocity) ==
>> -DiffusionTerm(var=dT, coeff=D))
>> eq = eqX & eqZ & eqT
>> # Solving Boussinesq equations
>> timeStepDuration = 10 * dL**2 / (2 * D)
>> steps = 50
>> for step in range(steps):
>> eq.solve(dt=timeStepDuration)
>> if viewer is not None:
>> viewer.plot()
>> plt.pause(0.1)
>> raw_input("… »)
>>
>> My error : all the input arrays must have same number of dimensions
>>
>> I’d like to obtain something like this :
>>
>>
>> If someone can help me, I will be very happy.
>>
>> Thanks a lot.
>>
>> Fabien G.
>> ___
>> fipy mailing list
>> fipy@nist.gov
>>
