tfracture_model.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
       ---
       tfracture_model.py (2409B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 """This script shows how to use the fracture density model in isolation, i.e. separately
            4 from the rest of PISM.
            5 
            6 Given an input file containing ice thickness, bed topography, x and y components of the
            7 SSA velocity (u_ssa_bc, v_ssa_bc), and the mask bc_mask (ones where fracture density is
            8 set, zero elsewhere) this code will run the fracture density model for a few time steps
            9 and save its outputs to a file.
           10 """
           11 
           12 import PISM
           13 ctx = PISM.Context()
           14 
           15 filename = ctx.config.get_string("input.file")
           16 
           17 # create the grid using the input file
           18 grid = PISM.IceGrid.FromFile(ctx.ctx, filename, ["thk"], PISM.CELL_CENTER)
           19 
           20 # initialize geometric data
           21 geometry = PISM.Geometry(grid)
           22 geometry.ice_thickness.regrid(filename)
           23 geometry.bed_elevation.regrid(filename)
           24 geometry.sea_level_elevation.set(0.0)
           25 geometry.ensure_consistency(0)
           26 
           27 # allocate the fracture density model
           28 flow_law = PISM.FlowLawFactory("stress_balance.ssa.", ctx.config, ctx.enthalpy_converter).create()
           29 fracture = PISM.FractureDensity(grid, flow_law)
           30 
           31 # initialize it using zero fracture age and density
           32 fracture.initialize()
           33 
           34 # read in the velocity field
           35 velocity = PISM.IceModelVec2V(grid, "_ssa_bc", PISM.WITHOUT_GHOSTS)
           36 velocity.set_attrs("", "x-component of the ice velocity", "m s-1", "m s-1", "", 0)
           37 velocity.set_attrs("", "y-component of the ice velocity", "m s-1", "m s-1", "", 1)
           38 velocity.regrid(filename)
           39 
           40 # find the longest time step we can take with this velocity field
           41 data = PISM.max_timestep_cfl_2d(geometry.ice_thickness, geometry.cell_type, velocity)
           42 dt = data.dt_max.value()
           43 
           44 # read in the BC mask
           45 bc_mask = PISM.IceModelVec2Int(grid, "bc_mask", PISM.WITHOUT_GHOSTS)
           46 bc_mask.regrid(filename)
           47 
           48 # Set hardness to a constant corresponding to -30C at 0 Pa
           49 hardness = PISM.IceModelVec2S(grid, "hardness", PISM.WITHOUT_GHOSTS)
           50 T = -30
           51 P = 0.0
           52 E = ctx.enthalpy_converter.enthalpy(T + 273.15, 0.0, P)
           53 hardness.set(flow_law.hardness(E, P))
           54 
           55 # take a few steps
           56 N = 100                         # arbitrary
           57 for k in range(N):
           58     fracture.update(dt, geometry, velocity, hardness, bc_mask)
           59 
           60 # save results
           61 output_filename = ctx.config.get_string("output.file_name")
           62 
           63 f = PISM.util.prepare_output(output_filename, append_time=True)
           64 fracture.density().write(f)
           65 fracture.age().write(f)
           66 fracture.growth_rate().write(f)
           67 fracture.flow_enhancement().write(f)
           68 fracture.healing_rate().write(f)
           69 fracture.toughness().write(f)
           70 f.close()