anuga-community / anuga_core

ANUGA for the simulation of the shallow water equation
https://anuga.anu.edu.au
Other
32 stars 23 forks source link

Errors in relative_time when using a realistic starttime #28

Closed stoiver closed 1 year ago

stoiver commented 1 year ago

There is a problem with catastrophic cancellation when using realistic starttime.

For instance the following code:

import anuga
import numpy as np
import matplotlib.pyplot as plt

from datetime import datetime

%matplotlib inline

print('# Setup domain with realistic starttime')
domain = anuga.rectangular_cross_domain(10,10)
Br = anuga.Reflective_boundary(domain)
domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
domain.set_timezone('Australia/Sydney')
domain.set_starttime(datetime(2022,9,29,12,30))
domain.set_quantity('elevation', lambda x,y : x)
domain.set_quantity('stage', lambda x,y : x+0.5)

print('# Evolve')
for t in domain.evolve(yieldstep=0.1, duration=1.0):
    domain.print_timestepping_statistics()
    print('    Relative time :',domain.relative_time)

produces the following:

# Setup domain with realistic starttime
# Evolve
Time = 1664418600.0000 (sec), steps=0 (0s)
    Relative time : 0.0
Time = 1664418600.1000 (sec), delta t in [0.00465004, 0.00677631] (s), steps=19 (0s)
    Relative time : 0.09999990463256836
Time = 1664418600.2000 (sec), delta t in [0.00362748, 0.00464594] (s), steps=25 (0s)
    Relative time : 0.19999980926513672
Time = 1664418600.3000 (sec), delta t in [0.00351577, 0.00390400] (s), steps=28 (0s)
    Relative time : 0.2999997138977051
Time = 1664418600.4000 (sec), delta t in [0.00390404, 0.00439568] (s), steps=25 (0s)
    Relative time : 0.39999961853027344
Time = 1664418600.5000 (sec), delta t in [0.00433682, 0.00457341] (s), steps=23 (0s)
    Relative time : 0.4999995231628418
Time = 1664418600.6000 (sec), delta t in [0.00439858, 0.00479898] (s), steps=22 (0s)
    Relative time : 0.5999994277954102
Time = 1664418600.7000 (sec), delta t in [0.00452521, 0.00499509] (s), steps=21 (0s)
    Relative time : 0.6999993324279785
Time = 1664418600.8000 (sec), delta t in [0.00499236, 0.00505003] (s), steps=20 (0s)
    Relative time : 0.7999992370605469
Time = 1664418600.9000 (sec), delta t in [0.00505604, 0.00508120] (s), steps=20 (0s)
    Relative time : 0.8999991416931152
Time = 1664418601.0000 (sec), delta t in [0.00504063, 0.00506525] (s), steps=20 (0s)
    Relative time : 0.9999990463256836
Time = 1664418601.0000 (sec), delta t = 0.00503902 (s), steps=1 (0s)
    Relative time : 1.0

Observe that Relative_time has inherited a large numerical error (cause by catastrophic cancellation). This example shows that this can lead to a strange extra step just before the finaltime.

To deal with this we will need to carefully reimplement the evolve procedures to work exclusively with relative_time .

stoiver commented 1 year ago

Resolved with PR #29