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 ]
