Re: [Opengrads-devel] Fishy finite differencing

[Michael F. <Mic...@no...>]

hi arlindo,

i'll take this on, for the lat-lon grids and if you can send me the refs 
and the algorithm, i'll give it a go.  can you point me to the precise 
fish udx code i need to work?  thanks, also it seem more consistent to 
leave the gauss-spec version as a separte udx. 

best /r mike

Arlindo da Silva wrote:
> Mike,
>
>   Or should I have said "finite differencing for fish()"?. (This is 
> intended primary for Mike; I'm cc'ing the list in case someone wants 
> to comment on the  algorithm.)
>
>    A little background. Fish() is a poisson solver that computes 
> streamfunction and velocity potential given vorticity and divervence:
>
>    psi = fish(vort)
>    chi = fish(divg)
>
> Unlike the classic UDF psichi,  I  have  separated  the computation of 
> vort/div from the possion solver proper. Like the classic psichi,  
> fish() uses fishpak.f, and old dusty-deck f77 eliptic solver that 
> stood the test of time. (Fishpak is not really f77 but rather f66, 
> with tons of goto's and arithmetic ifs, but it is very fast and the 
> numerics are solid.) Fishpak is freely available from netlib.
>
> Now, in principle vorticity and divergence could be computed with 
> grads intrinsics such as:
>
>   vort = hcurl(u,v)
>
> The problem is, hcurl() does not handles the boundaries very well, 
> placing UNDEFS in the first and last points (even if this is a global 
> periodic grid). Altough fish replaces undefs with 0. (like the classic 
> psichi), this is not the right thing to do. So, we need to provide a 
> replacement in the fish() package for properly handling the 
> boundaries, something that we could contribute later to COLA for 
> improving the current intrinsics.
>
> Here is my proposal for handling finite differencing at the boundaries:
>
> a) Zonal direction. If the grid is global, exploit the periodicidy, 
> and use centered differences for all points. If the grid (read: 
> dimension environment) is regional, use uncentered differences at the 
> boundaries.
>
> b) Meridional direction. If the grid (dimension environment) is 
> regional, use uncentered differences at the boundaries. Now, if the 
> grid is global, the handling of the extrema depends on whether the 
> pole is a grid point or not:
>
> i) Grid is off the pole by delta-lat/2. For example, on a 1 degree 
> grid the first and last gridpoints are -89.5 and 89.5. In this case 
> use uncentered differences at these points.
>
> ii) The poles are grid-points. This is the grid used by the 
> finite-volume dynamical core of S.J.-Lin used at GMAO, NCAR and GFDL.  
> The trick here is to consider these polar points as little polar caps 
> and use vector calculus to compute vorticity/divergence. Using Stokes' 
> theorem one can relate the mean vorticity in this polar cap to the 
> "circulation" at the boundaries of the polar cap, basically the sum of 
> "u"  at the boundary of the polar cap for some 2 pi factor thrown in 
> for good measure; for 1 deg grids, this boundary is at -89.5 and 89.5, 
> so "u" needs to be interpolated as the average of the first and second 
> meridion grid points (similarly for last and 1 from last gridpoints).  
> A similar algorithm is used for divergence, but now one relies on the 
> Gauss (divergence) theorem to relate the mean divergence at the polar 
> cap to the mass flux into the polar cap, basically the sum of "v" at 
> the boundary with some 2pi normalization factor.
>
> iii) Gaussian grids. The best solution would be to do a spectral 
> transform and compute vort/div that way.
>
> Poles are always messy, and usually a big source of trouble. So it 
> pays to be a bit careful.
>
> Mike: I added the stubs already to the fish() code base, perhaps more 
> than we need. I guess all that we need at this point is vort() and 
> div(). Could you add your code for the zonal direction? You can 
> implement the polar polar, and if you wish I can do it or write down 
> the algorithm more precisely (or send you the references.)
>
> I need to dig some spectral transform code for the gaussian grid case, 
> or else do a simple uncentered calculation  at the pole. If one is 
> going to bring spectral transforms in, then it is better to do the 
> calculation of psi/chi spectrally and not use fish() which is intended 
> for lat/lon grids. This may be another extension altogether.
>
> I'm cc'ing Brian in case he has any suggestions for us, as he may want 
> to incorporate some of this into hcur/divg later.
>
>     As always, any suggestions greatly appreciated.
>
>       Cheers!
>
>          Arlindo
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> -- 
> Arlindo da Silva
> da...@al... <mailto:da...@al...>
> ------------------------------------------------------------------------
>
> -------------------------------------------------------------------------
> This SF.net email is sponsored by: Microsoft
> Defy all challenges. Microsoft(R) Visual Studio 2008.
> http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/
> ------------------------------------------------------------------------
>
> _______________________________________________
> Opengrads-devel mailing list
> Ope...@li...
> https://lists.sourceforge.net/lists/listinfo/opengrads-devel
>