Ignore:
Timestamp:
10/11/19 00:27:00 (5 years ago)
Author:
bdef
Message:

CHG: syntax cahnge to meet most of Pep8 requirement

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/test/Par/RoundSheetShelf.py

    r22575 r24214  
    11import os.path
    2 import numpy
     2import numpy as np
    33import copy
    44import inspect
     
    88#Start defining model parameters here
    99
    10 di=md.materials.rho_ice/md.materials.rho_water
    11 rad=1.e6
    12 shelfextent=2.e5
     10di = md.materials.rho_ice / md.materials.rho_water
     11rad = 1.e6
     12shelfextent = 2.e5
    1313#Geometry
    14 hmin=300.
    15 hmax=1000.
    16 radius=numpy.sqrt(md.mesh.x*md.mesh.x+md.mesh.y*md.mesh.y.reshape(-1))
    17 ymin=numpy.min(radius)
    18 ymax=numpy.max(radius)
    19 md.geometry.thickness=hmax+(hmin-hmax)*(radius-ymin)/(ymax-ymin)
    20 md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
     14hmin = 300.
     15hmax = 1000.
     16radius = np.sqrt(md.mesh.x * md.mesh.x + md.mesh.y * md.mesh.y.reshape(- 1))
     17ymin = np.min(radius)
     18ymax = np.max(radius)
     19md.geometry.thickness = hmax + (hmin - hmax) * (radius - ymin) / (ymax - ymin)
     20md.geometry.base = -md.materials.rho_ice / md.materials.rho_water * md.geometry.thickness
    2121
    22 pos=numpy.nonzero(md.mask.groundedice_levelset>0.)[0]
    23 md.geometry.base[pos]=md.geometry.base[pos]-300.*(radius[pos]-(rad-shelfextent))/(rad-shelfextent)
    24 md.geometry.surface=md.geometry.base+md.geometry.thickness
     22pos = np.nonzero(md.mask.groundedice_levelset > 0.)[0]
     23md.geometry.base[pos] = md.geometry.base[pos] - 300. * (radius[pos] - (rad - shelfextent)) / (rad - shelfextent)
     24md.geometry.surface = md.geometry.base + md.geometry.thickness
    2525
    26 pos=numpy.nonzero(radius<200000.)
    27 md.geometry.thickness[pos]=100.
    28 md.geometry.base[pos]=-di*md.geometry.thickness[pos]-20.
    29 md.geometry.surface[pos]=md.geometry.base[pos]+md.geometry.thickness[pos]
     26pos = np.nonzero(radius < 200000.)
     27md.geometry.thickness[pos] = 100.
     28md.geometry.base[pos] = -di * md.geometry.thickness[pos] - 20.
     29md.geometry.surface[pos] = md.geometry.base[pos] + md.geometry.thickness[pos]
    3030
    31 pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.2*1.e6,md.mesh.x>-0.2*1.e6),md.mesh.y>0.))
    32 md.geometry.thickness[pos]=100.
    33 md.geometry.base[pos]=-di*md.geometry.thickness[pos]-20.
    34 md.geometry.surface[pos]=md.geometry.base[pos]+md.geometry.thickness[pos]
     31pos = np.nonzero(np.logical_and(np.logical_and(md.mesh.x < 0.2 * 1.e6, md.mesh.x > - 0.2 * 1.e6), md.mesh.y > 0.))
     32md.geometry.thickness[pos] = 100.
     33md.geometry.base[pos] = -di * md.geometry.thickness[pos] - 20.
     34md.geometry.surface[pos] = md.geometry.base[pos] + md.geometry.thickness[pos]
    3535
    36 pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.1*1.e6,md.mesh.x>-0.1*1.e6),numpy.logical_and(md.mesh.y<-0.5*1.e6,md.mesh.y>-0.6*1.e6)))
    37 md.geometry.thickness[pos]=100.
    38 md.geometry.base[pos]=-di*md.geometry.thickness[pos]-20.
    39 md.geometry.surface[pos]=md.geometry.base[pos]+md.geometry.thickness[pos]
     36pos = np.nonzero(np.logical_and(np.logical_and(md.mesh.x < 0.1 * 1.e6, md.mesh.x > - 0.1 * 1.e6), np.logical_and(md.mesh.y < - 0.5 * 1.e6, md.mesh.y > - 0.6 * 1.e6)))
     37md.geometry.thickness[pos] = 100.
     38md.geometry.base[pos] = -di * md.geometry.thickness[pos] - 20.
     39md.geometry.surface[pos] = md.geometry.base[pos] + md.geometry.thickness[pos]
    4040
    41 #plug holes into the ice sheet, to test for grounding line migration. 
    42 di=md.materials.rho_ice/md.materials.rho_water
    43 rad=numpy.sqrt(md.mesh.x**2+md.mesh.y**2)
    44 pos=numpy.nonzero(rad<200000.)
    45 md.geometry.thickness[pos]=100.
    46 md.geometry.base[pos]=-di*md.geometry.thickness[pos]-20.
    47 md.geometry.surface[pos]=md.geometry.base[pos]+md.geometry.thickness[pos]
     41#plug holes into the ice sheet, to test for grounding line migration.
     42di = md.materials.rho_ice / md.materials.rho_water
     43rad = np.sqrt(md.mesh.x**2 + md.mesh.y**2)
     44pos = np.nonzero(rad < 200000.)
     45md.geometry.thickness[pos] = 100.
     46md.geometry.base[pos] = -di * md.geometry.thickness[pos] - 20.
     47md.geometry.surface[pos] = md.geometry.base[pos] + md.geometry.thickness[pos]
    4848
    49 pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.2*1.e6,md.mesh.x>-0.2*1.e6),md.mesh.y>0.))
    50 md.geometry.thickness[pos]=100.
    51 md.geometry.base[pos]=-di*md.geometry.thickness[pos]-20.
    52 md.geometry.surface[pos]=md.geometry.base[pos]+md.geometry.thickness[pos]
     49pos = np.nonzero(np.logical_and(np.logical_and(md.mesh.x < 0.2 * 1.e6, md.mesh.x > - 0.2 * 1.e6), md.mesh.y > 0.))
     50md.geometry.thickness[pos] = 100.
     51md.geometry.base[pos] = -di * md.geometry.thickness[pos] - 20.
     52md.geometry.surface[pos] = md.geometry.base[pos] + md.geometry.thickness[pos]
    5353
    54 pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.1*1.e6,md.mesh.x>-0.1*1.e6),numpy.logical_and(md.mesh.y<-0.5*1.e6,md.mesh.y>-0.6*1.e6)))
    55 md.geometry.thickness[pos]=100.
    56 md.geometry.base[pos]=-di*md.geometry.thickness[pos]-20.
    57 md.geometry.surface[pos]=md.geometry.base[pos]+md.geometry.thickness[pos]
     54pos = np.nonzero(np.logical_and(np.logical_and(md.mesh.x < 0.1 * 1.e6, md.mesh.x > - 0.1 * 1.e6), np.logical_and(md.mesh.y < - 0.5 * 1.e6, md.mesh.y > - 0.6 * 1.e6)))
     55md.geometry.thickness[pos] = 100.
     56md.geometry.base[pos] = -di * md.geometry.thickness[pos] - 20.
     57md.geometry.surface[pos] = md.geometry.base[pos] + md.geometry.thickness[pos]
    5858
    59 #Initial velocity 
    60 md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
    61 md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
    62 md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
    63 md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
     59#Initial velocity
     60md.initialization.vx = np.zeros((md.mesh.numberofvertices))
     61md.initialization.vy = np.zeros((md.mesh.numberofvertices))
     62md.initialization.vz = np.zeros((md.mesh.numberofvertices))
     63md.initialization.pressure = np.zeros((md.mesh.numberofvertices))
    6464
    6565#Materials
    66 md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
    67 md.materials.rheology_B=paterson(md.initialization.temperature)
    68 md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
     66md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices))
     67md.materials.rheology_B = paterson(md.initialization.temperature)
     68md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements))
    6969
    7070#Surface mass balance and basal melting
    71 md.smb.mass_balance=-10.*numpy.ones((md.mesh.numberofvertices))
    72 md.basalforcings.groundedice_melting_rate=numpy.zeros((md.mesh.numberofvertices))
    73 pos=numpy.nonzero(md.mask.groundedice_levelset>0.)[0]
    74 md.basalforcings.groundedice_melting_rate[pos]=10.
    75 md.basalforcings.floatingice_melting_rate=numpy.zeros((md.mesh.numberofvertices))
    76 md.basalforcings.geothermalflux=numpy.ones((md.mesh.numberofvertices))
     71md.smb.mass_balance = -10. * np.ones((md.mesh.numberofvertices))
     72md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices))
     73pos = np.nonzero(md.mask.groundedice_levelset > 0.)[0]
     74md.basalforcings.groundedice_melting_rate[pos] = 10.
     75md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices))
     76md.basalforcings.geothermalflux = np.ones((md.mesh.numberofvertices))
    7777
    7878#Friction
    79 radius=1.e6
    80 shelfextent=2.e5
    81 md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices))
    82 xelem=numpy.mean(md.mesh.x[md.mesh.elements.astype(int)-1],axis=1)
    83 yelem=numpy.mean(md.mesh.y[md.mesh.elements.astype(int)-1],axis=1)
    84 rad=numpy.sqrt(xelem**2+yelem**2)
    85 flags=numpy.zeros(md.mesh.numberofelements)
    86 pos=numpy.nonzero(rad>=(radius-shelfextent))
    87 md.friction.coefficient[md.mesh.elements[pos,:]-1]=0.
    88 md.friction.p=numpy.ones((md.mesh.numberofelements))
    89 md.friction.q=numpy.ones((md.mesh.numberofelements))
     79radius = 1.e6
     80shelfextent = 2.e5
     81md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices))
     82xelem = np.mean(md.mesh.x[md.mesh.elements.astype(int) - 1], axis=1)
     83yelem = np.mean(md.mesh.y[md.mesh.elements.astype(int) - 1], axis=1)
     84rad = np.sqrt(xelem**2 + yelem**2)
     85flags = np.zeros(md.mesh.numberofelements)
     86pos = np.nonzero(rad >= (radius - shelfextent))
     87md.friction.coefficient[md.mesh.elements[pos, :] - 1] = 0.
     88md.friction.p = np.ones((md.mesh.numberofelements))
     89md.friction.q = np.ones((md.mesh.numberofelements))
    9090
    9191#Numerical parameters
    92 md.masstransport.stabilization=1
    93 md.thermal.stabilization=1
    94 md.verbose=verbose(0)
    95 md.settings.waitonlock=30
    96 md.stressbalance.restol=0.05
    97 md.stressbalance.reltol=0.05
    98 md.steadystate.reltol=0.05
    99 md.stressbalance.abstol=float('nan')
    100 md.timestepping.time_step=5.
    101 md.timestepping.final_time=5.
     92md.masstransport.stabilization = 1
     93md.thermal.stabilization = 1
     94md.verbose = verbose(0)
     95md.settings.waitonlock = 30
     96md.stressbalance.restol = 0.05
     97md.stressbalance.reltol = 0.05
     98md.steadystate.reltol = 0.05
     99md.stressbalance.abstol = float('nan')
     100md.timestepping.time_step = 5.
     101md.timestepping.final_time = 5.
    102102
    103103#bathymetry and grounding line migration:
    104 md.groundingline.migration='AggressiveMigration'
    105 md.geometry.bed=copy.deepcopy(md.geometry.base)
    106 pos=numpy.nonzero(md.mask.groundedice_levelset<0.)[0]
    107 md.geometry.bed[pos]=md.geometry.base[pos]-900.
     104md.groundingline.migration = 'AggressiveMigration'
     105md.geometry.bed = copy.deepcopy(md.geometry.base)
     106pos = np.nonzero(md.mask.groundedice_levelset < 0.)[0]
     107md.geometry.bed[pos] = md.geometry.base[pos] - 900.
    108108
    109109#Deal with boundary conditions:
    110 md.stressbalance.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices))
    111 md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices))
    112 md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices))
     110md.stressbalance.spcvx = float('nan') * np.ones((md.mesh.numberofvertices))
     111md.stressbalance.spcvy = float('nan') * np.ones((md.mesh.numberofvertices))
     112md.stressbalance.spcvz = float('nan') * np.ones((md.mesh.numberofvertices))
    113113
    114 pos=numpy.nonzero(numpy.logical_and(md.mesh.x==0,md.mesh.y==0))
    115 md.stressbalance.spcvx[pos]=0
    116 md.stressbalance.spcvy[pos]=0
     114pos = np.nonzero(np.logical_and(md.mesh.x == 0, md.mesh.y == 0))
     115md.stressbalance.spcvx[pos] = 0
     116md.stressbalance.spcvy[pos] = 0
    117117
    118 pos=numpy.nonzero(md.mesh.vertexonboundary)
    119 md.mask.ice_levelset[pos]=0
    120 md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices))
    121 md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices))
    122 md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
    123 md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
    124 md.thermal.spctemperature=737.*numpy.ones((md.mesh.numberofvertices))
     118pos = np.nonzero(md.mesh.vertexonboundary)
     119md.mask.ice_levelset[pos] = 0
     120md.balancethickness.spcthickness = float('nan') * np.ones((md.mesh.numberofvertices))
     121md.masstransport.spcthickness = float('nan') * np.ones((md.mesh.numberofvertices))
     122md.stressbalance.referential = float('nan') * np.ones((md.mesh.numberofvertices, 6))
     123md.stressbalance.loadingforce = 0 * np.ones((md.mesh.numberofvertices, 3))
     124md.thermal.spctemperature = 737. * np.ones((md.mesh.numberofvertices))
    125125
    126126#Change name so that no test have the same name
    127127if len(inspect.stack()) > 2:
    128         md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]
     128    md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]
Note: See TracChangeset for help on using the changeset viewer.