ttest_ssafd.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
       ---
       ttest_ssafd.py (2178B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 from PISMNC import PISMDataset as NC
            4 from optparse import OptionParser
            5 import numpy as np
            6 import subprocess
            7 import shlex
            8 import sys
            9 
           10 parser = OptionParser()
           11 
           12 parser.usage = "%prog [options]"
           13 parser.description = "Test the SSAFD solver using various geometric configurations."
           14 parser.add_option("-o", dest="output_file_prefix", default="ssafd_test",
           15                   help="output file prefix")
           16 parser.add_option("-L", dest="L", type=float, default=10.0,
           17                   help="horizontal domain dimensions, km")
           18 parser.add_option("-H", dest="H", type=float, default=500.0,
           19                   help="ice thickness in icy areas")
           20 
           21 (options, args) = parser.parse_args()
           22 
           23 M = 3
           24 
           25 
           26 def generate_input(N):
           27     """Generate a PISM bootstrapping file for a neighborhood number N."""
           28     output_filename = options.output_file_prefix + "_%08d.nc" % N
           29     pism_output_filename = options.output_file_prefix + "_o_%08d.nc" % N
           30     nc = NC(output_filename, 'w')
           31 
           32     format_string = "{0:0%db}" % (M * M)
           33     string = format_string.format(N)
           34 
           35     x = np.linspace(-options.L * 1000, options.L * 1000, M)
           36 
           37     zeros = np.zeros((M, M))
           38     thk = np.zeros(M * M)
           39 
           40     topg = zeros.copy() - 0.1 * options.H
           41     topg[1, 1] = -options.H
           42 
           43     for j in range(M * M):
           44         if string[j] == '1':
           45             thk[j] = options.H
           46 
           47     nc.create_dimensions(x, x)
           48     nc.write("topg", topg, attrs={"units": "m", "long_name": "bed_topography"})
           49     nc.write("climatic_mass_balance", zeros, attrs={"units": "kg m-2 year-1"})
           50     nc.write("ice_surface_temp", zeros, attrs={"units": "Celsius"})
           51     nc.write("thk", thk.reshape(M, M),
           52              attrs={"units": "m", "long_name": "land_ice_thickness"})
           53 
           54     nc.close()
           55 
           56     return output_filename, pism_output_filename
           57 
           58 
           59 def run_pismr(input_filename, output_filename):
           60     command = "pismr -i %s -bootstrap -o %s -Mx %d -My %d -Lz 1000 -Mz 5 -stress_balance ssa+sia -cfbc -y 0.001 -verbose 1" % (
           61         input_filename, output_filename, M, M)
           62     print("Running %s" % command)
           63     subprocess.call(shlex.split(command))
           64 
           65 
           66 for k in range(2 ** (M * M)):
           67     input_file, output_file = generate_input(k)
           68     run_pismr(input_file, output_file)