Changeset 25307


Ignore:
Timestamp:
07/29/20 00:45:32 (5 years ago)
Author:
jdquinn
Message:

CHG: Fix for WriteData; cleanup

Location:
issm/trunk-jpl
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/SMBgemb.py

    r24793 r25307  
    11import numpy as np
     2
     3from checkfield import checkfield
    24from fielddisplay import fielddisplay
    3 from checkfield import checkfield
     5from project3d import project3d
    46from WriteData import WriteData
    5 from project3d import project3d
    67
    78
    … …  
    1415    """
    1516
    16     def __init__(self):  # {{{
     17    def __init__(self, *args):  # {{{
    1718        #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived
    1819        #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number
    … …  
    3233
    3334        #inputs:
    34         self.Ta = float('NaN')       #2 m air temperature, in Kelvin
    35         self.V = float('NaN')        #wind speed (m/s-1)
    36         self.dswrf = float('NaN')    #downward shortwave radiation flux [W/m^2]
    37         self.dlwrf = float('NaN')    #downward longwave radiation flux [W/m^2]
    38         self.P = float('NaN')        #precipitation [mm w.e. / m^2]
    39         self.eAir = float('NaN')     #screen level vapor pressure [Pa]
    40         self.pAir = float('NaN')     #surface pressure [Pa]
    41         self.Tmean = float('NaN')    #mean annual temperature [K]
    42         self.Vmean = float('NaN')    #mean annual wind velocity [m s-1]
    43         self.C = float('NaN')        #mean annual snow accumulation [kg m-2 yr-1]
    44         self.Tz = float('NaN')       #height above ground at which temperature (T) was sampled [m]
    45         self.Vz = float('NaN')       #height above ground at which wind (V) eas sampled [m]
     35        self.Ta                     = np.nan    # 2 m air temperature, in Kelvin
     36        self.V                      = np.nan    # wind speed (m/s-1)
     37        self.dswrf                  = np.nan    # downward shortwave radiation flux [W/m^2]
     38        self.dlwrf                  = np.nan    # downward longwave radiation flux [W/m^2]
     39        self.P                      = np.nan    # precipitation [mm w.e. / m^2]
     40        self.eAir                   = np.nan    # screen level vapor pressure [Pa]
     41        self.pAir                   = np.nan    # surface pressure [Pa]
     42        self.Tmean                  = np.nan    # mean annual temperature [K]
     43        self.Vmean                  = np.nan    # mean annual wind velocity [m s-1]
     44        self.C                      = np.nan    # mean annual snow accumulation [kg m-2 yr-1]
     45        self.Tz                     = np.nan    # height above ground at which temperature (T) was sampled [m]
     46        self.Vz                     = np.nan    # height above ground at which wind (V) was sampled [m]
    4647
    4748        #optional inputs:
    48         self.aValue = float('NaN')  #Albedo forcing at every element.  Used only if aIdx == 0.
    49         self.teValue = float('NaN')  #Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
     49        self.aValue                 = np.nan    # Albedo forcing at every element.  Used only if aIdx == 0.
     50        self.teValue                = np.nan    # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
    5051
    5152        # Initialization of snow properties
    52         self.Dzini = float('NaN')    #cell depth (m)
    53         self.Dini = float('NaN')     #snow density (kg m-3)
    54         self.Reini = float('NaN')    #effective grain size (mm)
    55         self.Gdnini = float('NaN')   #grain dricity (0-1)
    56         self.Gspini = float('NaN')   #grain sphericity (0-1)
    57         self.ECini = float('NaN')    #evaporation/condensation (kg m-2)
    58         self.Wini = float('NaN')     #Water content (kg m-2)
    59         self.Aini = float('NaN')     #albedo (0-1)
    60         self.Tini = float('NaN')     #snow temperature (K)
    61         self.Sizeini = float('NaN')  #Number of layers
     53        self.Dzini                  = np.nan    # cell depth (m)
     54        self.Dini                   = np.nan    # snow density (kg m-3)
     55        self.Reini                  = np.nan    # effective grain size (mm)
     56        self.Gdnini                 = np.nan    # grain dricity (0-1)
     57        self.Gspini                 = np.nan    # grain sphericity (0-1)
     58        self.ECini                  = np.nan    # evaporation/condensation (kg m-2)
     59        self.Wini                   = np.nan    # Water content (kg m-2)
     60        self.Aini                   = np.nan    # albedo (0-1)
     61        self.Tini                   = np.nan    # snow temperature (K)
     62        self.Sizeini                = np.nan    # Number of layers
    6263
    6364        #settings:
    64         self.aIdx = float('NaN')     #method for calculating albedo and subsurface absorption (default is 1)
    65         self.swIdx = float('NaN')    # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
    66         self.denIdx = float('NaN')   #densification model to use (default is 2):
    67         self.dsnowIdx = float('NaN')  #model for fresh snow accumulation density (default is 1):
    68         self.zTop = float('NaN')     # depth over which grid length is constant at the top of the snopack (default 10) [m]
    69         self.dzTop = float('NaN')    # initial top vertical grid spacing (default .05) [m]
    70         self.dzMin = float('NaN')    # initial min vertical allowable grid spacing (default dzMin/2) [m]
    71         self.zY = float('NaN')       # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
    72         self.zMax = float('NaN')     #initial max model depth (default is min(thickness, 250)) [m]
    73         self.zMin = float('NaN')     #initial min model depth (default is min(thickness, 130)) [m]
    74         self.outputFreq = float('NaN')       #output frequency in days (default is monthly, 30)
     65        self.aIdx                   = np.nan    # method for calculating albedo and subsurface absorption (default is 1)
     66        self.swIdx                  = np.nan    # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
     67        self.denIdx                 = np.nan    # densification model to use (default is 2):
     68        self.dsnowIdx               = np.nan    # model for fresh snow accumulation density (default is 1):
     69        self.zTop                   = np.nan    # depth over which grid length is constant at the top of the snopack (default 10) [m]
     70        self.dzTop                  = np.nan    # initial top vertical grid spacing (default .05) [m]
     71        self.dzMin                  = np.nan    # initial min vertical allowable grid spacing (default dzMin/2) [m]
     72        self.zY                     = np.nan    # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
     73        self.zMax                   = np.nan    # initial max model depth (default is min(thickness, 250)) [m]
     74        self.zMin                   = np.nan    # initial min model depth (default is min(thickness, 130)) [m]
     75        self.outputFreq             = np.nan    # output frequency in days (default is monthly, 30)
    7576
    7677        #specific albedo parameters:
    7778        #Method 1 and 2:
    78         self.aSnow = float('NaN')    # new snow albedo (0.64 - 0.89)
    79         self.aIce = float('NaN')     # range 0.27-0.58 for old snow
     79        self.aSnow                  = np.nan    # new snow albedo (0.64 - 0.89)
     80        self.aIce                   = np.nan    # range 0.27-0.58 for old snow
    8081        #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
    81         self.cldFrac = float('NaN')  # average cloud amount
     82        self.cldFrac                = np.nan    # average cloud amount
    8283        #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
    83         self.t0wet = float('NaN')    # time scale for wet snow (15-21.9)
    84         self.t0dry = float('NaN')    # warm snow timescale (30)
    85         self.K = float('NaN')        # time scale temperature coef. (7)
    86         self.adThresh = float('NaN')  # Apply aIdx method to all areas with densities below this value,
     84        self.t0wet                  = np.nan    # time scale for wet snow (15-21.9)
     85        self.t0dry                  = np.nan    # warm snow timescale (30)
     86        self.K                      = np.nan    # time scale temperature coef. (7)
     87        self.adThresh               = np.nan    # Apply aIdx method to all areas with densities below this value,
    8788        # or else apply direct input value from aValue, allowing albedo to be altered.
    8889        # Default value is rho water (1023 kg m-3).
    8990
    9091        #densities:
    91         self.InitDensityScaling = float('NaN')       #initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
     92        self.InitDensityScaling     = np.nan    # initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
    9293
    9394        #thermo:
    94         self.ThermoDeltaTScaling = float('NaN')  #scaling factor to multiply the thermal diffusion timestep (delta t)
    95 
    96         self.steps_per_step = 1
    97         self.averaging = 0
    98         self.requested_outputs = []
     95        self.ThermoDeltaTScaling    = np.nan    # scaling factor to multiply the thermal diffusion timestep (delta t)
     96
     97        self.steps_per_step         = 1
     98        self.averaging              = 0
     99        self.requested_outputs      = []
    99100
    100101        #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM.
    … …  
    103104        #elev:  this is taken from the ISSM surface itself.
    104105
     106        nargin = len(args)
     107        if nargin:
     108            if nargin == 2:
     109                mesh = args[0]
     110                geometry = args[1]
     111                self.setdefaultparameters(mesh, geometry)
     112            else:
     113                raise Exception('constructor not supported: need mesh and geometry to set defaults')
    105114        #}}}
    106115
    … …  
    281290
    282291        md = checkfield(md, 'fieldname', 'smb.isgraingrowth', 'values', [0, 1])
    283         #md = checkfield(md, 'fieldname', 'smb.issmbgradients', 'values', [0, 1])
     292        md = checkfield(md, 'fieldname', 'smb.issmbgradients', 'values', [0, 1])
    284293        md = checkfield(md, 'fieldname', 'smb.isalbedo', 'values', [0, 1])
    285294        md = checkfield(md, 'fieldname', 'smb.isshortwave', 'values', [0, 1])
    … …  
    352361
    353362    def marshall(self, prefix, md, fid):    # {{{
    354 
    355363        yts = md.constants.yts
    356364

  • issm/trunk-jpl/src/m/solve/WriteData.m

    r24199 r25307  
    3232end
    3333
    34 %Scale data if necesarry
     34%Scale data if necessary
    3535if strcmpi(format,'MatArray')
    3636        for i=1:numel(data)

  • issm/trunk-jpl/src/m/solve/WriteData.py

    r25301 r25307  
    11from copy import deepcopy
    2 from struct import pack, error
     2from struct import error, pack
     3
    34import numpy as np
     5
    46import pairoptions
    57
    … …  
    4749    #       (see https://ross.ics.uci.edu/jenkins/view/All/job/Debian_Linux-Python/1036).
    4850    #
    49     if mattype != 0:
    50         data = deepcopy(data)
    51 
    52     #Scale data if necesarry
    53     if options.exist('scale'):
    54         data = np.array(data)
    55         scale = options.getfieldvalue('scale')
    56         if np.size(data) > 1 and np.ndim(data) > 1 and np.size(data, 0) == timeserieslength:
    57             data[0:-1, :] = scale * data[0:-1, :]
    58         else:
    59             data = scale * data
    60     if np.size(data) > 1 and np.size(data, 0) == timeserieslength:
    61         yts = options.getfieldvalue('yts')
    62         if np.ndim(data) > 1:
     51    # data = deepcopy(data)
     52
     53    #Scale data if necessary
     54    if datatype == 'MatArray':
     55        for i in range(np.size(data)):
     56            if options.exist('scale'):
     57                scale = options.getfieldvalue('scale')
     58                if np.ndim(data[i]) > 1 and data[i].shape[0] == timeserieslength:
     59                    data[i][:-2, :] = scale * data[i][:-2, :]
     60                else:
     61                    data[i] = scale * data[i]
     62            if np.ndim(data) > 1 and data[i].shape[0] == timeserieslength:
     63                yts = options.getfieldvalue('yts')
     64                data[i][-1, :] = yts * data[i][-1, :]
     65    else:
     66        if options.exist('scale'):
     67            scale = options.getfieldvalue('scale')
     68            if np.ndim(data) > 1 and data.shape[0] == timeserieslength:
     69                data[:-2, :] = scale * data[:-2, :]
     70            else:
     71                data = scale * data
     72        if np.ndim(data) > 1 and data.shape[0] == timeserieslength:
     73            yts = options.getfieldvalue('yts')
    6374            data[-1, :] = yts * data[-1, :]
    64         else:
    65             data[-1] = yts * data[-1]
    6675
    6776    #Step 1: write the enum to identify this record uniquely

  • issm/trunk-jpl/test/NightlyRun/test234.py

    r25112 r25307  
    1010from setflowequation import *
    1111from solve import *
    12 from SMBgemb import *
    1312from partitioner import *
    1413from dmeth_params_set import *
    … …  
    3938#partitioning
    4039npart = 20
    41 partition = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on') - 1
     40partition = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on')[0] - 1
    4241
    4342#variables

  • issm/trunk-jpl/test/NightlyRun/test235.py

    r25133 r25307  
    1010from setflowequation import *
    1111from solve import *
    12 from SMBgemb import *
    1312from partitioner import *
    1413from dmeth_params_set import *

  • issm/trunk-jpl/test/NightlyRun/test243.py

    r24795 r25307  
    11#Test Name: SquareShelfSMBGemb
     2from socket import gethostname
     3import sys
     4
    25import numpy as np
    3 import sys
     6
    47from model import *
    5 from socket import gethostname
    6 from triangle import *
    7 from setmask import *
    88from parameterize import *
    99from setflowequation import *
     10from setmask import *
     11from SMBgemb import *
    1012from solve import *
    11 from SMBgemb import *
     13from triangle import *
    1214
    1315md = triangle(model(), '../Exp/Square.exp', 350000.)
    … …  
    1921
    2022#Use of Gemb method for SMB computation
    21 md.smb = SMBgemb()
    22 md.smb.setdefaultparameters(md.mesh, md.geometry)
     23md.smb = SMBgemb(md.mesh, md.geometry)
    2324md.smb.dsnowIdx = 1
    2425
    … …  
    3738md.smb.eAir = np.append(np.tile(np.conjugate(inputs[b'eAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0)
    3839md.smb.pAir = np.append(np.tile(np.conjugate(inputs[b'pAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0)
    39 md.smb.pAir = np.append(np.tile(np.conjugate(inputs[b'pAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0)
    4040md.smb.Vz = np.tile(np.conjugate(inputs[b'LP']['Vz']), (md.mesh.numberofelements, 1)).flatten()
    4141md.smb.Tz = np.tile(np.conjugate(inputs[b'LP']['Tz']), (md.mesh.numberofelements, 1)).flatten()
    … …  
    4444
    4545#smb settings
    46 md.smb.requested_outputs = ['SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbEC', 'SmbA', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF', 'SmbMeanULW', 'SmbNetLW', 'SmbNetSW']
     46md.smb.requested_outputs = [
     47    'SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbEC',
     48    'SmbA', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF',
     49    'SmbMeanULW', 'SmbNetLW', 'SmbNetSW'
     50    ]
    4751
    4852#only run smb core:
    … …  
    5458md.timestepping.start_time = 1965.
    5559md.timestepping.final_time = 1966.
    56 md.timestepping.time_step = 1. / 365.0
     60md.timestepping.time_step = 1.0 / 365
    5761md.timestepping.interp_forcings = 0.
    5862
    … …  
    6872field_tolerances = [1e-12, 2e-12, 1e-12, 1e-11, 1e-11, 2e-11, 1e-11, 1e-12, 1e-11, 1e-12, 1e-12, 1e-12, 1e-11, 2e-11, 2e-11, 1e-11, 9e-10, 2e-11]
    6973#shape is different in python solution (fixed using reshape) which can cause test failure:
    70 field_values = [nlayers,
    71                 md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1),
    72                 md.results.TransientSolution[-1].SmbT[0, 0:nlayers].reshape(1, -1),
    73                 md.results.TransientSolution[-1].SmbD[0, 0:nlayers].reshape(1, -1),
    74                 md.results.TransientSolution[-1].SmbRe[0, 0:nlayers].reshape(1, -1),
    75                 md.results.TransientSolution[-1].SmbGdn[0, 0:nlayers].reshape(1, -1),
    76                 md.results.TransientSolution[-1].SmbGsp[0, 0:nlayers].reshape(1, -1),
    77                 md.results.TransientSolution[-1].SmbA[0, 0:nlayers].reshape(1, -1),
    78                 md.results.TransientSolution[-1].SmbEC[0],
    79                 md.results.TransientSolution[-1].SmbMassBalance[0],
    80                 md.results.TransientSolution[-1].SmbMAdd[0],
    81                 md.results.TransientSolution[-1].SmbDzAdd[0],
    82                 md.results.TransientSolution[-1].SmbFAC[0],
    83                 md.results.TransientSolution[-1].SmbMeanSHF[0],
    84                 md.results.TransientSolution[-1].SmbMeanLHF[0],
    85                 md.results.TransientSolution[-1].SmbMeanULW[0],
    86                 md.results.TransientSolution[-1].SmbNetLW[0],
    87                 md.results.TransientSolution[-1].SmbNetSW[0]]
     74field_values = [
     75    nlayers,
     76    md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1),
     77    md.results.TransientSolution[-1].SmbT[0, 0:nlayers].reshape(1, -1),
     78    md.results.TransientSolution[-1].SmbD[0, 0:nlayers].reshape(1, -1),
     79    md.results.TransientSolution[-1].SmbRe[0, 0:nlayers].reshape(1, -1),
     80    md.results.TransientSolution[-1].SmbGdn[0, 0:nlayers].reshape(1, -1),
     81    md.results.TransientSolution[-1].SmbGsp[0, 0:nlayers].reshape(1, -1),
     82    md.results.TransientSolution[-1].SmbA[0, 0:nlayers].reshape(1, -1),
     83    md.results.TransientSolution[-1].SmbEC[0],
     84    md.results.TransientSolution[-1].SmbMassBalance[0],
     85    md.results.TransientSolution[-1].SmbMAdd[0],
     86    md.results.TransientSolution[-1].SmbDzAdd[0],
     87    md.results.TransientSolution[-1].SmbFAC[0],
     88    md.results.TransientSolution[-1].SmbMeanSHF[0],
     89    md.results.TransientSolution[-1].SmbMeanLHF[0],
     90    md.results.TransientSolution[-1].SmbMeanULW[0],
     91    md.results.TransientSolution[-1].SmbNetLW[0],
     92    md.results.TransientSolution[-1].SmbNetSW[0]
     93    ]

  • issm/trunk-jpl/test/NightlyRun/test244.py

    r25133 r25307  
    2828
    2929# Use of Gemb method for SMB computation
    30 md.smb = SMBgemb()
    31 md.smb.setdefaultparameters(md.mesh, md.geometry)
     30md.smb = SMBgemb(md.mesh, md.geometry)
    3231md.smb.dsnowIdx = 0
    3332

  • issm/trunk-jpl/test/NightlyRun/test250.py

    r25133 r25307  
    99from setflowequation import *
    1010from solve import *
    11 from SMBgemb import *
    1211from partitioner import *
    1312from dmeth_params_set import *

  • issm/trunk-jpl/test/NightlyRun/test251.py

    r25133 r25307  
    99from setflowequation import *
    1010from solve import *
    11 from SMBgemb import *
    1211from partitioner import *
    1312from dmeth_params_set import *

  • issm/trunk-jpl/test/NightlyRun/test252.py

    r24773 r25307  
    11#Test Name: SquareShelfSMBGembClim
     2from socket import gethostname
     3import sys
     4
    25import numpy as np
    3 import sys
     6
    47from model import *
    5 from socket import gethostname
    6 from triangle import *
    7 from setmask import *
    88from parameterize import *
    99from setflowequation import *
     10from setmask import *
     11from SMBgemb import *
    1012from solve import *
    11 from SMBgemb import *
     13from triangle import *
    1214
    1315md = triangle(model(), '../Exp/Square.exp', 350000.)
    … …  
    1921
    2022#Use of Gemb method for SMB computation
    21 md.smb = SMBgemb()
    22 md.smb.setdefaultparameters(md.mesh, md.geometry)
     23md.smb = SMBgemb(md.mesh, md.geometry)
    2324md.smb.dsnowIdx = 4
    2425

  • issm/trunk-jpl/test/NightlyRun/test253.py

    r24760 r25307  
    1919
    2020#Use of Gemb method for SMB computation
    21 md.smb = SMBgemb()
    22 md.smb.setdefaultparameters(md.mesh, md.geometry)
     21md.smb = SMBgemb(md.mesh, md.geometry)
    2322md.smb.dsnowIdx = 3
    2423md.smb.aIdx = 2

  • issm/trunk-jpl/test/NightlyRun/test418.py

    r25133 r25307  
    3737# - Run valgrind and fix the above
    3838#
    39 partition, md = partitioner(md, 'package', 'chaco', 'npart', npart)
    40 partition -= 1
     39partition = partitioner(md, 'package', 'chaco', 'npart', npart)[0] - 1
    4140
    4241vector = np.arange(1, 1 + md.mesh.numberofvertices, 1).reshape(-1, 1)

  • issm/trunk-jpl/test/NightlyRun/test420.py

    r25113 r25307  
    1818#partitioning
    1919npart = 10
    20 partition = partitioner(md, 'package', 'chaco', 'npart', npart) - 1
     20partition = partitioner(md, 'package', 'chaco', 'npart', npart)[0] - 1
    2121md.qmu.isdakota = 1
    2222

  • issm/trunk-jpl/test/NightlyRun/test444.py

    r25176 r25307  
    1313from setflowequation import *
    1414from solve import *
    15 from SMBgemb import *
    1615from partitioner import *
    1716from dmeth_params_set import *

  • issm/trunk-jpl/test/Par/SquareShelf.py

    r24862 r25307  
     1from arch import *
     2import inspect
    13import os.path
    2 import inspect
    3 from arch import *
     4
    45import numpy as np
    5 from verbose import verbose
     6
    67from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
    78from paterson import paterson
    89from SetIceShelfBC import SetIceShelfBC
     10from verbose import verbose
    911
    1012#Start defining model parameters here
    … …  
    3436#dbg - end
    3537
     38# x = x[0][:]
     39# y = y[0][:]
     40# vx = vx[0][:]
     41# vy = vy[0][:]
     42# index = index[0][:]
     43
    3644[md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
    3745[md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
    38 md.initialization.vz = np.zeros((md.mesh.numberofvertices))
    39 md.initialization.pressure = np.zeros((md.mesh.numberofvertices))
     46x = y = vx = vy = index = None
     47md.initialization.vz = np.zeros((md.mesh.numberofvertices, 1))
     48md.initialization.pressure = np.zeros((md.mesh.numberofvertices, 1))
    4049
    4150#dbg - begin
    … …  
    5059#dbg - end
    5160
     61# Materials
     62md.initialization.temperature = (273 - 20) * np.ones((md.mesh.numberofvertices, 1))
     63md.materials.rheology_B = paterson(md.initialization.temperature)
     64md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1))
    5265
    53 #Materials
    54 md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices))
    55 md.materials.rheology_B = paterson(md.initialization.temperature)
    56 md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements))
     66# Friction
     67md.friction.coefficient = 20 * np.ones((md.mesh.numberofvertices, 1))
     68md.friction.coefficient[np.where(md.mask.ocean_levelset < 0)[0]] = 0
     69md.friction.p = np.ones((md.mesh.numberofelements, 1))
     70md.friction.q = np.ones((md.mesh.numberofelements, 1))
    5771
    58 #Friction
    59 md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices))
    60 md.friction.coefficient[np.nonzero(md.mask.ocean_levelset < 0.)[0]] = 0.
    61 md.friction.p = np.ones((md.mesh.numberofelements))
    62 md.friction.q = np.ones((md.mesh.numberofelements))
    63 
    64 #Numerical parameters
    65 md.masstransport.stabilization = 1.
    66 md.thermal.stabilization = 1.
     72# Numerical parameters
     73md.masstransport.stabilization = 1
     74md.thermal.stabilization = 1
    6775md.settings.waitonlock = 30
    6876md.verbose = verbose()
    … …  
    7078md.steadystate.reltol = 0.02
    7179md.stressbalance.reltol = 0.02
    72 md.stressbalance.abstol = float('nan')
    73 md.timestepping.time_step = 1.
    74 md.timestepping.final_time = 3.
     80md.stressbalance.abstol = np.nan
     81md.timestepping.time_step = 1
     82md.timestepping.final_time = 3
    7583md.groundingline.migration = 'None'
    7684
    77 #Boundary conditions:
    78 #  #md = SetIceShelfBC(md)
     85# Boundary conditions:
    7986md = SetIceShelfBC(md, '../Exp/SquareFront.exp')
    8087
Note: See TracChangeset for help on using the changeset viewer.