Bug in pyclaw/src/pyclaw/classic/flux3.f90 with capacity function scaling when transverse_waves==21 or 22

[weslowrie]

At line 418 of flux3.f90 the gadd() variable is multiplied by dtdx1d and dtdz. dtdxtd has been divided by the capacity function in step3.f90, while dtdz has not. I believe this is an error.

gadd(m,2,0,i) = gadd(m,2,0,i) + (1.d0/6.d0)*dtdx1d(i)*dtdz * (bpcpapdq(m,i) - bpcmapdq(m,i))

I noticed that when the capacity function is large, the solver fails, but when the computational to physical mapping is identity, it works fine. The failure only occurs when transverse_waves==21 or 22 (m4 > 0 in flux3.f90).

If I modify the scaling such that it does not have a large capacity function, but is non-unity, it does not fail, but the answer is incorrect. After modifying flux3.f90 by dividing dtdz by the capacity function, it does not fail. This is also true for the H fluxes and dtdy needs to be modified.

I could be missing something here?

I suspect this is also an issue in amrclaw?

[mandli]

@gradylemoine I don't suppose you could look at this and confirm? You have more recent experience with this code and could more quickly confirm the bug and bug fix.

[gradylemoine]

I'm a little busy lately with my job and various other distractions. I might be able to handle this more efficiently (though it's not obvious I could), but by the time you wait for me to get to it, somebody who's never seen the code before could probably run it down.

On Sat, Aug 23, 2014 at 10:04 AM, Kyle Mandli [email protected] wrote:

@gradylemoine https://github.com/gradylemoine I don't suppose you could look at this and confirm? You have more recent experience with this code and could more quickly confirm the bug and bug fix.

— Reply to this email directly or view it on GitHub https://github.com/clawpack/pyclaw/issues/457#issuecomment-53158836.

[mandli]

Thanks Grady for getting back so quickly.

After looking at this myself, I think that the current code is correct. Looking at the update formulas I think that the capacity of each cell should only be incorporated once into each update flux which is done through dtdx1d. This conclusion is based on a rereading of @rjleveque's red book, chapter 23 where the capacities of the cells are there only included once which makes some sense although I may be misinterpreting what you think might be the problem.

[weslowrie]

@mandli I agree that if you look at equation 23.6 in @rjleveque red book that for each term the capacity function is included once. But it is unclear to me how the high-resolution correction terms work in the case with capacity functions. Is there a reference where the equation: G += 1/6*dt/dx/kappa * dt/dz * ((B^+C^+A^+\Delta Q) - (B^+C^-A^+\Delta Q)) OR line 418 of flux3.f90: gadd(m,2,0,i) = gadd(m,2,0,i) + (1.d0/6.d0)*dtdx1d(i)*dtdz * (bpcpapdq(m,i) - bpcmapdq(m,i)) is derived?

[mandli]

That's a great question. I took a look through [1] and [2] but could not find a definitive answer.

1. LeVeque, R. J. Wave Propagation Algorithms for Multidimensional Hyperbolic Systems. Journal of Computational Physics 131, 327–353 (1997).
2. Langseth, J. O. & LeVeque, R. J. A Wave Propagation Method for Three-Dimensional Hyperbolic Conservation Laws. Journal of Computational Physics 165, 126–166 (2000).

[gradylemoine]

If it's too complicated to easily figure out, why not take a test case with a known, smooth, analytic solution, and do a convergence study each way on a smooth mapped grid? One way should ruin your order of accuracy. It's not at all an elegant way to decide how to proceed, but it's a way to cut the knot if the problem is too confusing.

--Grady

On Mon, Aug 25, 2014 at 9:42 AM, Kyle Mandli [email protected] wrote:

That's a great question. I took a look through [1] and [2] but could not find a definitive answer.

1. LeVeque, R. J. Wave Propagation Algorithms for Multidimensional Hyperbolic Systems. Journal of Computational Physics 131, 327–353 (1997).
2. Langseth, J. O. & LeVeque, R. J. A Wave Propagation Method for Three-Dimensional Hyperbolic Conservation Laws. Journal of Computational Physics 165, 126–166 (2000).

— Reply to this email directly or view it on GitHub https://github.com/clawpack/pyclaw/issues/457#issuecomment-53289541.

[weslowrie]

I'm not sure that would expose the potential problem. I have done something similar...

I computed a solution to a problem without the capacity-form differencing and did this for several solver settings (transverse_waves = 20, 21, and 22). Then I compute the same problem with capacity-form differencing, but with an identity mapping to make sure I get the same exact solution (to machine precision) for the different solver settings. Then finally, I compute the solution with capacity-form differencing with a non-identity mapping, but yielding the same physical grid as the original problem. This is what exposed the error. I see the same exact solution for transverse_waves = 20, and a different solution with 21, and 22. If the mapping is aggressive enough, the solution fails on the first time-step. The aggressive mapping starts with a computational grid of 0.0 to 1.0, and has a physical grid from 0.0 to 2.0e8.

The fact that the solution is exactly the same for transverse_waves=20, but not for 21 and 22 I thought was pretty definitive that something was up with the algorithm. These things can be tricky though, so there still could be another explanation?