clawpack / geoclaw

Version of Clawpack for geophysical waves and flows
http://www.clawpack.org/geoclaw
BSD 3-Clause "New" or "Revised" License
76 stars 87 forks source link

Conservation in GeoClaw #415

Open donnaaboise opened 4 years ago

donnaaboise commented 4 years ago

@mjberger has observed a relative conservation error of about 1e-5 in GeoClaw. The test case she used was a modified the Chile example, initialized it with a single hump of water in the center of a single grid and solid wall boundary conditions.

I took this example, and modified it slightly to get the following results, which also show an apparent loss of conservation. I compared the conservation on a flat domain with results on the latlong grid. I adjust the radius of the earth (r = 63) so that the total mass in both the flat and latlong case were the same, and use a dam-break type initial condition (disk of height 0.1 over bathymetry=-1)

I've pushed my example here. (Branch : donnaaboise/conservation_test), under examples/tsunami/conservation_test.

Conservation error on the flat domain

   original total mass ...
time t =  0.00000E+00,  total mass =  0.363135709869270E+04  diff =  0.0000E+00
time t =  0.20000E-01,  total mass =  0.363135709869271E+04  diff =  0.1364E-11
time t =  0.24831E+00,  total mass =  0.363135709869271E+04  diff =  0.6366E-11
time t =  0.47645E+00,  total mass =  0.363135709869272E+04  diff =  0.1728E-10
time t =  0.70119E+00,  total mass =  0.363135709869272E+04  diff =  0.1819E-10
time t =  0.92531E+00,  total mass =  0.363135709869273E+04  diff =  0.2137E-10
time t =  0.11492E+01,  total mass =  0.363135709869272E+04  diff =  0.1819E-10
time t =  0.13730E+01,  total mass =  0.363135709869273E+04  diff =  0.2274E-10
time t =  0.15968E+01,  total mass =  0.363135709869273E+04  diff =  0.2137E-10
time t =  0.18206E+01,  total mass =  0.363135709869273E+04  diff =  0.2183E-10
time t =  0.20446E+01,  total mass =  0.363135709869273E+04  diff =  0.2910E-10

The relative conservation error is approximately 1e-15 at the final time.

Conservation on the lat/long domain

   original total mass ...
time t =  0.00000E+00,  total mass =  0.363271169391532E+04  diff =  0.0000E+00
time t =  0.20000E-01,  total mass =  0.363271169391532E+04  diff = -0.1819E-11
time t =  0.15364E+00,  total mass =  0.363271146718917E+04  diff = -0.2267E-03
time t =  0.28728E+00,  total mass =  0.363270972617553E+04  diff = -0.1968E-02
time t =  0.42092E+00,  total mass =  0.363270647552242E+04  diff = -0.5218E-02
time t =  0.55456E+00,  total mass =  0.363270171647609E+04  diff = -0.9977E-02
time t =  0.68820E+00,  total mass =  0.363269545018405E+04  diff = -0.1624E-01
time t =  0.82184E+00,  total mass =  0.363268767774769E+04  diff = -0.2402E-01
time t =  0.95549E+00,  total mass =  0.363267840025930E+04  diff = -0.3329E-01
time t =  0.10891E+01,  total mass =  0.363266761881659E+04  diff = -0.4408E-01
time t =  0.12228E+01,  total mass =  0.363265533452964E+04  diff = -0.5636E-01

The relative conservation error is approximately 1e-5 at the final time.

Here are the two views of the solution.

flat latlong

With no initial perturbation, there is no conservation error, so boundary conditions are working correctly. I've also checked to see that the solution does not reach the boundaries within the 10 time steps shown.

Has anyone else seen this behavior? Or have a possible explanation?

mandli commented 4 years ago

I assume this is also demonstrable in amrclaw?

donnaaboise commented 4 years ago

Are you referring to a simple SWE example in AMRClaw? I haven't tried that. It would be interesting to see if using a `rpn2_shallow_bathymetry_fwave.f90' has the same issue.

mandli commented 4 years ago

I guess I was thinking that this might be replicated in general for anything with a non-flat mapping and qad?

donnaaboise commented 4 years ago

The problem shouldn't occur in AMRClaw on a single grid, for say, a simple SWE solver. But with AMR, on a non-flat grid, there could be a loss of conservation related to how qad is set up.

mandli commented 4 years ago

So it's completely based on a qad problem?

donnaaboise commented 4 years ago

On a curved mapped grid with AMR, the way that QAD currently calls the Riemann solver to get flux differences can lead to a conservation error. But the apparent loss of conservation in GeoClaw is on a single grid, so QAD isn't involved at all.

mandli commented 4 years ago

Ah ok. I think I am getting confused by multiple conversations.