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 ]
