ttemp_continuity.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
 (HTM) git clone git://src.adamsgaard.dk/pism
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       ttemp_continuity.py (1928B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 from sys import exit, argv, stderr
            4 from os import system
            5 from numpy import squeeze, abs, diff
            6 
            7 try:
            8     from netCDF4 import Dataset as NC
            9 except:
           10     print("netCDF4 is not installed!")
           11     sys.exit(1)
           12 
           13 pism_path = argv[1]
           14 mpiexec = argv[2]
           15 
           16 stderr.write("Testing: temperature continuity at ice-bed interface (polythermal case).\n")
           17 
           18 cmd = "%s %s/pismv -test F -y 10 -verbose 1 -o bar-temp-continuity.nc" % (mpiexec, pism_path)
           19 stderr.write(cmd + '\n')
           20 
           21 e = system(cmd)
           22 if e != 0:
           23     exit(1)
           24 
           25 deltas = []
           26 dts = [200, 100]                # FIXME: this is fragile and the test fails if I add smaller dt like 50 here
           27 for dt in dts:
           28     cmd = "%s %s/pisms -eisII B -y 2400 -Mx 16 -My 16 -Mz 21 -Lbz 1000 -Mbz 11 -energy enthalpy -regrid_file bar-temp-continuity.nc -regrid_vars thk -verbose 1 -max_dt %f -o foo-temp-continuity.nc -o_size big" % (
           29         mpiexec, pism_path, dt)
           30     stderr.write(cmd + '\n')
           31 
           32     e = system(cmd)
           33     if e != 0:
           34         exit(1)
           35 
           36     e = system("ncks -O -v temp -d z,0 foo-temp-continuity.nc temp-temp-continuity.nc")
           37     if e != 0:
           38         exit(1)
           39 
           40     e = system("ncks -O -v litho_temp -d zb,10 foo-temp-continuity.nc litho_temp-temp-continuity.nc")
           41     if e != 0:
           42         exit(1)
           43 
           44     nc1 = NC("temp-temp-continuity.nc")
           45     nc2 = NC("litho_temp-temp-continuity.nc")
           46 
           47     temp = squeeze(nc1.variables['temp'][:])
           48     litho_temp = squeeze(nc2.variables['litho_temp'][:])
           49 
           50     deltas.append(abs(temp - litho_temp).max())
           51 
           52 # these deltas are observed to decrease O(dt^1) approximately, which is expected from theory
           53 for (dt, delta) in zip(dts, deltas):
           54     stderr.write("dt = %f, delta = %f\n" % (dt, delta))
           55 
           56 # the only test is whether they decrease; no rate measured
           57 if any(diff(deltas) > 0):
           58     print(diff(deltas))
           59     exit(1)
           60 
           61 system("rm foo-temp-continuity.nc foo-temp-continuity.nc~ bar-temp-continuity.nc temp-temp-continuity.nc litho_temp-temp-continuity.nc")
           62 exit(0)