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) > <jonathan.gu...@nist.gov> 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("<Enter to continue>...") > # 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("<Enter to finsh>… ") > > >> On Aug 14, 2018, at 11:22 AM, fgendr01 <fabien.gend...@univ-lr.fr> 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) : >> <PastedGraphic-1.png> >> 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("<Enter to continue>...") >> # 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("<Enter to finsh>… ») >> >> My error : all the input arrays must have same number of dimensions >> >> I’d like to obtain something like this : >> <PastedGraphic-3.png> >> >> If someone can help me, I will be very happy. >> >> Thanks a lot. >> >> Fabien G. >> _______________________________________________ >> fipy mailing list >> fipy@nist.gov >> http://www.ctcms.nist.gov/fipy >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]