Dion - I'm glad this was helpful to you.
For future reference, help(DiffusionTerm) gives the mathematical form that the term represents. http://www.ctcms.nist.gov/fipy/documentation/numerical/equation.html and the following pages describe each of the terms and how they are discretized in more detail. The change I made certainly seems to have resolved the stability problems and looks qualitatively consistent with the paper you cite, but you'll want to work through the results carefully and check their error against whatever metrics you have. I know Daniel spent a lot of time working through the SIMPLE algorithm and ultimately didn't feel even it was the final word, but it's the best we can do for flow problems right now. - Jon On Nov 1, 2015, at 11:09 AM, Dion Häfner <[email protected]> wrote: > Hey Jon, > > that looks amazing, thank you very much! It seems like all problems are > resolved for now with the new treatment of the advection term. I was in > fact not aware how exactly (if at all) Fipy was treating partial > derivatives of the convection coefficients, but your answer makes > perfect sense. I will run some tests and implement the SIMPLE scheme if > I should run into any more trouble. > > Best regards, > Dion > > On 31.10.2015 19:03, Guyer, Jonathan E. Dr. wrote: >> Dion - >> >> Thanks for your example and explanation. I've converted to an IPython >> notebook to make it easier to show the math derivation. See >> https://github.com/guyer/haefner/blob/master/killworth.ipynb (math rendering >> is kind of random on the github site, but should be OK if you download it >> and have MathJax installed). >> >>> * Is it really necessary to have the face velocity as a separate variable? >>> Could I not somehow use the cell velocities and achieve better results >>> (solve for the convection coefficient implicitly)? >> I think it should be possible to declare velocity as a rank-1 cell variable, >> write a single velocity equation, and just declare fVelocity = >> velocity.faceValue, but I can't get this to work right now. There may be a >> reason that Daniel didn't write any of his flow examples to work that way. >> >>> * Did I translate the equations correctly? I was not sure if the convection >>> terms really represent the form from the paper. >> There's an additional term missing related to changing "u ux + v uy" from >> advective form to FiPy's conservative convective form. See the IPython >> notebook for the specifics. >> >> I also think there's a sign error in coriolisTermY, but that probably just >> flips you into the southern hemisphere. >> >>> * Which type of convection term should I use? I do not really care much >>> about speed, I just want an accurate solution. >> Daniel uses a CentralDifferenceConvectionTerm in his reactive wetting >> examples. I don't know why, although I'm sure he explained it to me at one >> point. I don't do much with flow in my work. >> >>> * Is the height.grad term correct? Could I somehow write this as an >>> (implicit) source, or is this done automatically? >> I tried doing it as a ConvectionTerm, but it introduced a lot of ringing >> artifacts. Compare >> >> https://github.com/guyer/haefner/blob/ffd31134748555afd5eff415356e90cc8a10d240/output/killworth/Column%20height-10.00days.png >> >> to >> >> https://github.com/guyer/haefner/blob/711fdd6154ae6c4beedd0cbe64385f258b51764f/output/killworth/Column%20height-10.00days.png >> >> >> >>> * Should I implement a SIMPLE algorithm (as in the cavity flow example)? >> Yes, I think this is the right way to couple the height and velocity >> variables. >> >> - Jon >> >> >> On Oct 29, 2015, at 12:32 PM, Dion Häfner <[email protected]> wrote: >> >>> Dear Fipy developers, >>> >>> >>> >>> I am currently trying to implement a 2D shallow water model in Fipy. This >>> has worked out fine so far, but now I am facing some instability, and don't >>> really know how to resolve that. The model is taken from >>> >>> >>> http://journals.ametsoc.org/doi/pdf/10.1175/1520-0485%281991%29021%3C1581%3ACEGA%3E2.0.CO%3B2 >>> (eq. 6.1 - 6.3) >>> >>> >>> and I am trying to reproduce figure 17 (same parameters as in the paper). I >>> implemented the PDE system as follows: >>> >>> >>> >>> # VARIABLE SETUP >>> height = CellVariable(mesh=mesh, name='Column height', hasOld=True) >>> xVelocity = CellVariable(mesh=mesh, name='U', value=0., hasOld=True) >>> yVelocity = CellVariable(mesh=mesh, name='V', value=0., hasOld=True) >>> fVelocity = FaceVariable(mesh=mesh, name="Face velocity", rank=1, value=0.) >>> >>> [ set initial conditions... ] >>> >>> # EQUATION SETUP >>> coriolisTermX = ImplicitSourceTerm(var=yVelocity,coeff=.5*mesh.y) >>> coriolisTermY = ImplicitSourceTerm(var=xVelocity,coeff=.5*mesh.y) >>> if par["Ah"] != 0: >>> xVelocityEq = TransientTerm(var=xVelocity) + \ >>> ConvectionTerm(coeff=fVelocity,var=xVelocity) - \ >>> coriolisTermX - \ >>> DiffusionTerm(coeff=par["Ah"],var=xVelocity) + \ >>> height.grad.dot([1.,0.]) \ >>> == 0 >>> yVelocityEq = TransientTerm(var=yVelocity) + \ >>> ConvectionTerm(coeff=fVelocity,var=yVelocity) - \ >>> coriolisTermY - \ >>> DiffusionTerm(coeff=par["Ah"],var=yVelocity) + \ >>> height.grad.dot([0.,1.]) \ >>> == 0 >>> heightEq = TransientTerm(var=height) + \ >>> ConvectionTerm(coeff=fVelocity,var=height) \ >>> == 0 >>> eqSystem = xVelocityEq & yVelocityEq & heightEq # couple equations >>> >>> # BOUNDARY CONDITIONS >>> xVelocity.constrain(0., mesh.exteriorFaces) >>> yVelocity.constrain(0., mesh.exteriorFaces) >>> height.faceGrad.constrain(0., mesh.exteriorFaces) >>> >>> >>> and eventually solve it by calling >>> >>> t = 0. >>> while t < (par["t1"] * DAY_TO_T): >>> xVelocity.updateOld() >>> yVelocity.updateOld() >>> height.updateOld() >>> residual = 1 >>> while residual > par["absTol"]: >>> fVelocity[0] = xVelocity.arithmeticFaceValue >>> fVelocity[1] = yVelocity.arithmeticFaceValue >>> fVelocity[..., mesh.exteriorFaces.value] = 0. >>> residual = eqSystem.sweep(dt = dt) >>> t += dt >>> >>> >>> >>> As you see, I am treating each velocity component as its own variable and >>> then couple all equations together. The system is then solved by sweeping >>> and setting the face velocity in every iteration, since I need it for the >>> convection terms. But when I run this with the parameters from the paper, >>> the model becomes unstable, and all water eventually piles up in one corner >>> of the domain. Thus, a couple of questions: >>> >>> >>> * Is it really necessary to have the face velocity as a separate variable? >>> Could I not somehow use the cell velocities and achieve better results >>> (solve for the convection coefficient implicitly)? >>> >>> * Did I translate the equations correctly? I was not sure if the convection >>> terms really represent the form from the paper. >>> >>> * Which type of convection term should I use? I do not really care much >>> about speed, I just want an accurate solution. >>> >>> * Is the height.grad term correct? Could I somehow write this as an >>> (implicit) source, or is this done automatically? >>> >>> * Should I implement a SIMPLE algorithm (as in the cavity flow example)? >>> >>> >>> I have attached the full version of the script and the parameter file I >>> used. This should be runnable without further packages. >>> >>> >>> Thanks a bunch! >>> >>> Dion Häfner >>> <adjustment.py><parameters.ini>_______________________________________________ >>> fipy mailing list >>> [email protected] >>> http://www.ctcms.nist.gov/fipy >>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >> >> _______________________________________________ >> fipy mailing list >> [email protected] >> http://www.ctcms.nist.gov/fipy >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
