PETSC Laplacian Boundary condition

[totork]

Hi all!

I am currently fighting with some simulations that include laplace inversions with PETSC, specifically with frequent inversion fails when I want to run with the inner flag 2 (no local gradient) and the outer flag 16 (setting from value). I tried to look a little bit through the code and found the following code

BOUT-dev/src/invert/laplace/impls/petsc/petsc_laplace.cxx

Lines 331 to 341 in 9eb94c7

	// Set boundary conditions
	if (!isInnerBoundaryFlagSet(INVERT_RHS)) {
	BOUT_FOR_SERIAL(index, indexer->getRegionInnerX()) {
	rhs(index) = isInnerBoundaryFlagSet(INVERT_SET) ? x0[index] : 0.0;
	}
	}
	if (!isOuterBoundaryFlagSet(INVERT_RHS)) {
	BOUT_FOR_SERIAL(index, indexer->getRegionOuterX()) {
	rhs(index) = isInnerBoundaryFlagSet(INVERT_SET) ? x0[index] : 0.0;
	}
	}

Does this work as intended? I have to look a little bit more through the code to understand if there is a reason for this, but I found it counter-intuitive that i check the inner boundary when setting the outer boundary

[ZedThree]

Oh no, that looks like my fault, sorry! Yes, it should be isOuterBoundaryFlagSet -- clearly our tests are not good enough on this

[totork]

On another note, besides this bug I still observe that the solver (pretty much always) breaks down running inner neumann boundary conditions. I tried to compare our "old" bout++ implementation that we currently use for the FCI branch of Hermes-3. There the inner boundary is implemented via

BOUT-dev/src/invert/laplace/impls/petsc/petsc_laplace.cxx

Lines 411 to 415 in 6a9b24c

	Element(i, x, z, 0, 0,
	-1.0 / coords->dx(x, y, z) / sqrt(coords->g_11(x, y, z)), MatA);
	Element(i, x, z, 1, 0,
	1.0 / coords->dx(x, y, z) / sqrt(coords->g_11(x, y, z)), MatA);
	Element(i, x, z, 2, 0, 0.0, MatA);

While the current next branch of BOUT++ looks like

BOUT-dev/src/invert/laplace/impls/petsc/petsc_laplace.cxx

Lines 489 to 503 in 673665d

	if (inner_X_neumann) {
	const auto dx = sliceXZ(coords->dx, y);
	const auto g11 = sliceXZ(coords->g11, y);
	BOUT_FOR_SERIAL(i, indexer->getRegionInnerX()) {
	const auto factor = 1. / dx[i] / std::sqrt(g11[i]);
	operator2D(i, i) = -factor;
	operator2D(i, i.xp()) = factor;
	}
	} else {
	BOUT_FOR_SERIAL(i, indexer->getRegionInnerX()) {
	operator2D(i, i) = 0.5;
	operator2D(i, i.xp()) = 0.5;
	}
	}

They look similar at first glance, but the major difference i see is that next uses g11, while the old implementation uses g_11. Maybe @dschwoerer knows more details about this.