Make mixed second-derivative operators more correct

[johnomotani]

Previously, D2DYDZ was using an explicit loop, doing 2nd order centred differences without taking into account whether yup/ydown should be used. Instead it should be OK to just use DDZ(DDY(f)) since the z-grid spacing is constant and the y-grid spacing does not depend on z, so the order of derivatives is not important.

[d7919]

Given that we're working towards (or at least thinking about) treating z the same as other directions (i.e. could be non-uniform etc.) is this introducing a "problem" we will have to fix? Could we switch the order of the derivatives at the expense of a communicate?

[johnomotani]

I think we definitely have a "problem", but it's not being introduced here as such: I think in principle all the mixed derivatives should have 2 forms, e.g. D2DXDY and D2DYDX to support variable grid spacing. At the moment we get away with this because usually g12==0 so the D2DXDY terms don't contribute, and the mixed z-derivatives are OK because dz is constant and dx, dy don't depend on z.

The problem with D2DXDY isn't simple to fix, because communication isn't enough. To calculate D2DXDY(f) you really need to:

1. calculate dfdy = DDY(f)
2. communicate dfdy
3. apply a boundary condition to dfdy consistent with the boundary condition on f
4. return DDX(dfdy)

I don't know how generally we can apply a sensible boundary condition. If f's x-boundary condition is zero-value Dirichlet or Neumann, I think you could apply the same one to dfdy, but anything more complicated would be tricky. For example if you had a Dirichlet boundary condition that varied in the y-direction, you'd have to take a y-derivative of the boundary condition, and then apply that to dfdy...

If/when we'd make z the same as other directions, the mixed z-derivatives would have the same problem as D2DXDY. That's why I think it's OK to rely on special-ness of the z-direction for now, because if we fix the problem above, we could then apply it to all mixed derivatives at the same time. If this is a good argument, it could be worth changing D2DXDZ to take the DDX first (like this implementation of D2DYDZ with Y->X), as it would simplify the function a bit?

[bendudson]

Mathematically D2DXDY should be equal to D2DYDX unless the coordinates are non-commutative (which the shifted metric coordinates are). It shouldn't matter if the grid is non-uniform. One way might be to calculate the derivatives in both orders (X then Y, and Y then X) and take the average of the result.

[johnomotani]

I'd made a mistake about the non-uniform grids (see bottom). I think DDX(DDY(f)) and DDY(DDX(f)) should both converge to the correct answer. They're not identical numerically though:

DDX(DDY(f))
= DDX( (f[i,j+1] - f[i,j-1]) / (2*dy[i,j]) )
= (f[i+1,j+1] - f[i+1,j-1]) / (4*dx[i,j]*dy[i+1,j]) - (f[i-1,j+1] - f[i-1,j-1]) / (4*dx[i,j]*dy[i-1,j])

DDY(DDX(f))
= DDY( (f[i+1,j] - f[i-1,j]) / (2*dx[i,j]) )
= (f[i+1,j+1] - f[i-1,j+1]) / (4*dx[i,j+1]*dy[i,j]) - (f[i+1,j-1] - f[i-1,j-1]) / (4*dx[i,j-1]*dy[i,j])

My mistake about the non-uniform grids was because I was doing some testing previously (the stuff in #1315) that I think I set up wrong. I was trying to test non-uniform grids, but was initializing fields to test using expressions, and put in a scaling factor for varying dy/dx. I was taking symbolic derivatives of that scaling factor, which made me think it the order of derivatives mattered, when I think I just shouldn't have been: I should actually have been doing something like 'real Y' #1131.

[ZedThree]

Was there a consensus on the correct thing to do here?

[johnomotani]

Was talking to Ben about this yesterday. We think:

• In a fixed, globally defined, coordinate system derivatives commute, so analytically ddx(ddy(f))=ddy(ddx(f))
• Complications in shifted-metric or FCI are because we are in a field-aligned coordinate system, but only locally - the coordinates are defined around one y-value but at the next or previous grid point in y, the local coordinate system is different
• In the local coordinate system, derivatives still commute - and f.yup()/f.ydown() are the values of f in the local field-aligned coordinate system
  • if we define dfdx_up = DDX(f.yup()) and dfdx_down = DDX(f.ydown()) then the mixed derivative (dfdx_up[jy+1] - dfdx_down[jy-1]) / (2*dy[jy]) is equivalent to DDX(DDY(f))[jy]*
  • similarly if dfdz_up = DDZ(f.yup()) and dfdz_down = DDZ(f.ydown()) then (dfdz_yup[jy+1] - dfdz_ydown[jy-1]) / (2*dy[jy]) is equivalent to DDZ(DDY(f))[jy]
    • for ShiftedMetric f.yup() and f.ydown() are only shifted by an angle in z, so I think DDY(DDZ(f)) should also be the same, but this would not be true for FCI

The conclusion is that D2DX2(f) = DDX(DDY(f)) and D2DYDZ(f) = DDZ(DDY(f)) should be correct, and we should use those forms as they are simpler to calculate than the alternative with yup/ydown.

I'll polish up this PR soon...

--
* or at least they both converge to the same thing at 2nd order accuracy. The points where dx and dy are evaluated are slightly different for the two orders, but the differences should be higher order corrections, especially as dx is constant in y and dy is constant in x, so we don't get the non-uniform grid corrections, like the term in D2DX2 due to DDX(dx).

[johnomotani]

I've added to this PR:

• D2DXDZ does DDZ(DDX(f)) - this means we don't have to communicate or set boundary conditions
• D2DXDY takes the y-derivative first, then the x-derivative, as discussed above.
  • applies a boundary condition to the intermediate dfdy before taking DDX; by default uses free_o3, but this can be changed with an argument to D2DXDY.
• Section added to the manual to record this discussion.

[ZedThree]

Thanks @johnomotani, I love a PR that includes documentation updates! 😄