Add unit test for LaplaceXZPetsc (3D metrics)

[ZedThree]

Doesn't currently work, but I've probably messed up either the forward operator or the boundary conditions. Copied and stripped down from the 3D AMG unit test.

[ZedThree]

Just rebased on top of the latest 3D metrics branch. @bshanahan please could you take a look at the forward operator? It's different from the version in the test-laplacexz integrated test, but that test doesn't work in 3D (I'm not sure why it doesn't!)

[bshanahan]

Looks good to me. Thanks @ZedThree!

[dschwoerer]

Did you forget to add the CMakeLists.txt?

 CMake Error at tests/integrated/CMakeLists.txt:26 (add_subdirectory):
  The source directory

    /home/runner/work/BOUT-dev/BOUT-dev/tests/integrated/test-laplacexz

  does not contain a CMakeLists.txt file.

[ZedThree]

Oops, good catch!

Fiddling with the integrated test, I noticed it fails when using a more complicated geometry taken from another test:

# circular, large-aspect-ratio geometry
eps = .1
r = eps*(.9 + .1*x)
theta = y-pi
Rxy = 1. + r*cos(theta)
Bt = 1./Rxy
Bp = 1/r
Bxy = sqrt(Bt^2 + Bp^2)
hthe = r
J = hthe/Bp
nu = Bt*hthe/(Bp*Rxy)

dx = .1/nx
dy = 2.*pi/ny
dz = 2.*pi/nz

g11 = (Rxy*Bp)^2
g22 = 1./hthe^2
g33 = (Bxy/(Rxy*Bp))^2
g12 = 0.
g13 = 0.0
g23 = -nu/hthe^2
g_11 = 1./(Rxy*Bp)^2
g_22 = (Bxy*hthe/Bp)^2
g_33 = Rxy^2
g_12 = 0.
g_13 = 0.
g_23 = Bt*hthe*Rxy/Bp

which I think is because the forward operator used is:

Field3D g = coord->g11 * D2DX2(f) + coord->g13 * D2DXDZ(f) + coord->g33 * D2DZ2(f);

and I think this assumes A = 1, B = 0, J = 1, and there's no spatial variation in any of the metric components.

Using this form in the unit test, with A = 1, B = 0 does get much closer, but still not correct. And the version of solve that starts with an initial guess of zero is much worse.

So I'm still not clear why the unit test is failing. I think the forward operator used:

    const auto AJ = A * coords->J;
    const auto xx_coef = AJ * coords->g11;
    const auto zz_coef = AJ * coords->g33;
    const auto xz_coef = AJ * coords->g13;
    const auto ddx_f = DDX(f);
    const auto ddz_f = DDZ(f);

    const auto xx = (DDX(xx_coef) * ddx_f) + (xx_coef * D2DX2(f));
    const auto zz = (DDZ(zz_coef) * ddz_f) + (zz_coef * D2DZ2(f));
    const auto xz = (DDX(xz_coef) * ddz_f) + (xz_coef * D2DXDZ(f));
    const auto zx = (DDZ(xz_coef) * ddx_f) + (xz_coef * D2DXDZ(f));

    auto result = ((xx + zz + xz + zx) / coords->J) + (B * f);

reduces to the integrated test form (factor 2 in the cross-term? doesn't seem to make much difference removing zx though), but even setting A = 1, B = 0 doesn't make it pass. Could just be not enough points? Being a unit test, it's a very small grid.

[johnomotani]

which I think is because the forward operator used is:

Field3D g = coord->g11 * D2DX2(f) + coord->g13 * D2DXDZ(f) + coord->g33 * D2DZ2(f);

If I remember right, LaplaceXZ uses a finite-volume discretisation, which is different from the normal finite-difference one, so D2DX2, etc. won't be quite the right forward operators. Operators from fv_ops.hxx might be better (@bendudson would know)?

On the unit test - if the forward operator is exactly correct, the grid size shouldn't matter. If it's a slightly different discretisation, then the coarse grid will make the errors large.

[ZedThree]

Ah, that could very well be it. Is Div_a_Laplace_perp the correct FV operator?

[bendudson]

Ah, that could very well be it. Is Div_a_Laplace_perp the correct FV operator?

@ZedThree Yes I think so

[ZedThree]

Naively setting the unit test forward operator to return FV::Div_a_Laplace_perp(A, f) + (B * f); causes the solver to not converge: LaplaceXZ failed to converge. Reason -9 which is KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)

[ZedThree]

I think there might actually be a bug in LaplaceXZPetsc itself. In the integrated test, I compare the following versions of the forward operator:

  // Note additional B*f term compared to existing implementation in integrated test
  Field3D g_original_forward = coord->g11 * D2DX2(f) + coord->g13 * D2DXDZ(f) + coord->g33 * D2DZ2(f) + (B * f);
  Field3D g_fv_forward = FV::Div_a_Laplace_perp(A, f) + (B * f);
  Field3D g_full_forward = ((xx + zz + xz + zx) / coord->J) + (B * f);
  Field3D g_cutdown_forward = ((xx_coef * D2DX2(f) + zz_coef * D2DZ2(f) + xz_coef * D2DXDZ(f)) / coord->J) + (B * f);

Just removing coord->g13 = 1.8 (and using the identity metric) from the integrated test is enough to make them all pass! The original forward operator then fails when changing the coefficients from A = 1, B = 0, whereas the others all pass.

Now here's the results of using A = 1, B = 0, setting each metric component to a non-identity value (i.e. start with g11 = g22 = g33 = 1, g12 = g13 = g23 = 0):

g11 = (Rxy*Bp)^2:
      Checking tolerance of f_original_forward... Fail, maximum difference = 0.06169788056809944
      Checking tolerance of f_fv_forward... Fail, maximum difference = 5.462163618297211e-05
      Checking tolerance of f_full_forward... Fail, maximum difference = 0.00011438781469941262
      Checking tolerance of f_cutdown_forward... Fail, maximum difference = 0.06169788056809655
g22 = 1./hthe^2
      Checking tolerance of f_original_forward... Fail, maximum difference = 0.0563677865727088
      Checking tolerance of f_fv_forward... Pass
      Checking tolerance of f_full_forward... Pass
      Checking tolerance of f_cutdown_forward... Fail, maximum difference = 0.056367786572708134
g33 = (Bxy/(Rxy*Bp))^2
      Checking tolerance of f_original_forward... Fail, maximum difference = 0.0563677865727088
      Checking tolerance of f_fv_forward... Pass
      Checking tolerance of f_full_forward... Pass
      Checking tolerance of f_cutdown_forward... Fail, maximum difference = 0.056367786572708134
g12 = -nu/hthe^2
      Checking tolerance of f_original_forward... Fail, maximum difference = 0.018805880619605864
      Checking tolerance of f_fv_forward... Pass
      Checking tolerance of f_full_forward... Fail, maximum difference = 1.644526634159149e-05
      Checking tolerance of f_cutdown_forward... Fail, maximum difference = 0.018805880619597426
g13 = -nu/hthe^2
      Checking tolerance of f_original_forward... Fail, maximum difference = 0.026260950470192057
      Checking tolerance of f_fv_forward... Fail, maximum difference = 0.9135161012218411
      Checking tolerance of f_full_forward... Fail, maximum difference = 0.9138612111108522
      Checking tolerance of f_cutdown_forward... Fail, maximum difference = 0.026260950470189837
g23 = -nu/hthe^2
      Checking tolerance of f_original_forward... Fail, maximum difference = 0.021845718373473932
      Checking tolerance of f_fv_forward... Fail, maximum difference = 0.003322584546862295
      Checking tolerance of f_full_forward... Fail, maximum difference = 1.7919812282940484e-05
      Checking tolerance of f_cutdown_forward... Fail, maximum difference = 0.02184571837347593

What I think this all means is that, yes, the finite volume operator is the correct forward operator, but there is at least one bug in either FV::Div_a_Laplace_perp or in LaplaceXZpetsc -- possibly both -- when handling g11, g13, g23.

Setting the diagonal metric elements to scalars is not enough to make it fail, it looks like they need to spatially varying. Setting the off-diagonal components to arbitrary scalars in the input file is tricky, because the covariant terms can be negative.

Looking at the code: FV::Div_a_Laplace_perp only has one-sided differences in positive x and z, whereas LaplaceXZPetsc has the full 9-point stencil. Also FV::Div_a_Laplace_perp has terms for the flux in Y, whereas they obviously don't appear in LaplaceXZPetsc.

That's probably about as far as I'm going to get, I think it needs someone who knows the operators better than me to investigate further.

[ZedThree]

Reimplementing the integrated test in next and using the finite volume operator as the forward operation, it fails if g11 != 1., so more evidence of a bug in either the finite volume operator or LaplaceXZPetsc.

I think we might need a test with an analytic answer to work out where the issue is.

Bugs persist

[bshanahan]

It seems like the INVERT_SET and INVERT_RHS boundary conditions have identical code. Even if this isn't the bug, it means laplacexz-petsc needs closer review.

Curiouser and curiouser -- changing the definition of row (for instance here) in each of the X-boundary calls solves the boundary issue that plagued the old integrated test, but doesn't do much otherwise.

[dschwoerer]

Does it make sense to close this / don't wait on this, and merge #2359 once it is working?

Having #2025 merged would make my life easier 👍

[ZedThree]

Yes, that's sensible. We know that there's an issue with either LaplaceXZ or the forward operator in next, and the 3D metrics don't really change that.

Closing this in favour of #2359