Changeset 24214 for issm/trunk-jpl/test/Par/RoundSheetShelf.py
- Timestamp:
- 10/11/19 00:27:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/test/Par/RoundSheetShelf.py
r22575 r24214 1 1 import os.path 2 import numpy 2 import numpy as np 3 3 import copy 4 4 import inspect … … 8 8 #Start defining model parameters here 9 9 10 di =md.materials.rho_ice/md.materials.rho_water11 rad =1.e612 shelfextent =2.e510 di = md.materials.rho_ice / md.materials.rho_water 11 rad = 1.e6 12 shelfextent = 2.e5 13 13 #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.thickness14 hmin = 300. 15 hmax = 1000. 16 radius = np.sqrt(md.mesh.x * md.mesh.x + md.mesh.y * md.mesh.y.reshape(- 1)) 17 ymin = np.min(radius) 18 ymax = np.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 21 21 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.thickness22 pos = np.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 25 25 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]26 pos = np.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] 30 30 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]31 pos = 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.)) 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] 35 35 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]36 pos = 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))) 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] 40 40 41 #plug holes into the ice sheet, to test for grounding line migration. 42 di =md.materials.rho_ice/md.materials.rho_water43 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. 42 di = md.materials.rho_ice / md.materials.rho_water 43 rad = np.sqrt(md.mesh.x**2 + md.mesh.y**2) 44 pos = np.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] 48 48 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]49 pos = 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.)) 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] 53 53 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]54 pos = 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))) 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] 58 58 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 60 md.initialization.vx = np.zeros((md.mesh.numberofvertices)) 61 md.initialization.vy = np.zeros((md.mesh.numberofvertices)) 62 md.initialization.vz = np.zeros((md.mesh.numberofvertices)) 63 md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) 64 64 65 65 #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))66 md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices)) 67 md.materials.rheology_B = paterson(md.initialization.temperature) 68 md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements)) 69 69 70 70 #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))71 md.smb.mass_balance = -10. * np.ones((md.mesh.numberofvertices)) 72 md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices)) 73 pos = np.nonzero(md.mask.groundedice_levelset > 0.)[0] 74 md.basalforcings.groundedice_melting_rate[pos] = 10. 75 md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) 76 md.basalforcings.geothermalflux = np.ones((md.mesh.numberofvertices)) 77 77 78 78 #Friction 79 radius =1.e680 shelfextent =2.e581 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))79 radius = 1.e6 80 shelfextent = 2.e5 81 md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices)) 82 xelem = np.mean(md.mesh.x[md.mesh.elements.astype(int) - 1], axis=1) 83 yelem = np.mean(md.mesh.y[md.mesh.elements.astype(int) - 1], axis=1) 84 rad = np.sqrt(xelem**2 + yelem**2) 85 flags = np.zeros(md.mesh.numberofelements) 86 pos = np.nonzero(rad >= (radius - shelfextent)) 87 md.friction.coefficient[md.mesh.elements[pos, :] - 1] = 0. 88 md.friction.p = np.ones((md.mesh.numberofelements)) 89 md.friction.q = np.ones((md.mesh.numberofelements)) 90 90 91 91 #Numerical parameters 92 md.masstransport.stabilization =193 md.thermal.stabilization =194 md.verbose =verbose(0)95 md.settings.waitonlock =3096 md.stressbalance.restol =0.0597 md.stressbalance.reltol =0.0598 md.steadystate.reltol =0.0599 md.stressbalance.abstol =float('nan')100 md.timestepping.time_step =5.101 md.timestepping.final_time =5.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. 102 102 103 103 #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.104 md.groundingline.migration = 'AggressiveMigration' 105 md.geometry.bed = copy.deepcopy(md.geometry.base) 106 pos = np.nonzero(md.mask.groundedice_levelset < 0.)[0] 107 md.geometry.bed[pos] = md.geometry.base[pos] - 900. 108 108 109 109 #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))110 md.stressbalance.spcvx = float('nan') * np.ones((md.mesh.numberofvertices)) 111 md.stressbalance.spcvy = float('nan') * np.ones((md.mesh.numberofvertices)) 112 md.stressbalance.spcvz = float('nan') * np.ones((md.mesh.numberofvertices)) 113 113 114 pos =numpy.nonzero(numpy.logical_and(md.mesh.x==0,md.mesh.y==0))115 md.stressbalance.spcvx[pos] =0116 md.stressbalance.spcvy[pos] =0114 pos = np.nonzero(np.logical_and(md.mesh.x == 0, md.mesh.y == 0)) 115 md.stressbalance.spcvx[pos] = 0 116 md.stressbalance.spcvy[pos] = 0 117 117 118 pos =numpy.nonzero(md.mesh.vertexonboundary)119 md.mask.ice_levelset[pos] =0120 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))118 pos = np.nonzero(md.mesh.vertexonboundary) 119 md.mask.ice_levelset[pos] = 0 120 md.balancethickness.spcthickness = float('nan') * np.ones((md.mesh.numberofvertices)) 121 md.masstransport.spcthickness = float('nan') * np.ones((md.mesh.numberofvertices)) 122 md.stressbalance.referential = float('nan') * np.ones((md.mesh.numberofvertices, 6)) 123 md.stressbalance.loadingforce = 0 * np.ones((md.mesh.numberofvertices, 3)) 124 md.thermal.spctemperature = 737. * np.ones((md.mesh.numberofvertices)) 125 125 126 126 #Change name so that no test have the same name 127 127 if len(inspect.stack()) > 2: 128 128 md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]
Note:
See TracChangeset
for help on using the changeset viewer.