Changeset 24214 for issm/trunk-jpl/test/NightlyRun/test1108.py
- Timestamp:
- 10/11/19 00:27:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/test/NightlyRun/test1108.py
r23793 r24214 13 13 14 14 """ 15 This test is a test from the ISMP -HOM Intercomparison project.15 This test is a test from the ISMP - HOM Intercomparison project. 16 16 Pattyn and Payne 2006 17 17 """ … … 22 22 23 23 for L in L_list: 24 nx = 30 24 nx = 30 #numberof nodes in x direction 25 25 ny = 30 26 26 md = model() 27 27 md = squaremesh(md, L, L, nx, ny) 28 md = setmask(md, '', '') 28 md = setmask(md, '', '') #ice sheet test 29 29 md = parameterize(md, '../Par/ISMIPD.py') 30 30 md.extrude(10, 1.) … … 32 32 md = setflowequation(md, 'HO', 'all') 33 33 34 34 #We need one grd on dirichlet: the 4 corners are set to zero 35 35 md.stressbalance.spcvx = np.nan * np.ones((md.mesh.numberofvertices)) 36 36 md.stressbalance.spcvy = np.nan * np.ones((md.mesh.numberofvertices)) … … 48 48 md.stressbalance.spcvz[pos] = 0. 49 49 50 50 #Create MPCs to have periodic boundary conditions 51 51 posx = np.nonzero(md.mesh.x == 0.)[0] 52 52 posx2 = np.nonzero(md.mesh.x == np.max(md.mesh.x))[0] 53 53 54 posy = np.intersect1d(np.intersect1d(np.where(md.mesh.y == 0.), np.where(md.mesh.x != 0.)), np.where(md.mesh.x != np.max(md.mesh.x)))[0] 54 posy = np.intersect1d(np.intersect1d(np.where(md.mesh.y == 0.), np.where(md.mesh.x != 0.)), np.where(md.mesh.x != np.max(md.mesh.x)))[0] #Don't take the same nodes two times 55 55 posy2 = np.intersect1d(np.intersect1d(np.where(md.mesh.y == np.max(md.mesh.y)), np.where(md.mesh.x != 0.)), np.where(md.mesh.x != np.max(md.mesh.x)))[0] 56 56 57 md.stressbalance.vertex_pairing = np.vstack((np.hstack((posx.reshape(- 1, 1) + 1, posx2.reshape(-1, 1) + 1)), np.hstack((posy.reshape(-1, 1) + 1, posy2.reshape(-1, 1) + 1))))57 md.stressbalance.vertex_pairing = np.vstack((np.hstack((posx.reshape(- 1, 1) + 1, posx2.reshape(- 1, 1) + 1)), np.hstack((posy.reshape(- 1, 1) + 1, posy2.reshape(- 1, 1) + 1)))) 58 58 59 59 #Compute the stressbalance 60 60 md.cluster = generic('name', gethostname(), 'np', 8) 61 61 md.verbose = verbose('convergence', True) … … 64 64 md.stressbalance.abstol = np.nan 65 65 md.stressbalance.vertex_pairing = np.empty((0, 2)) 66 66 #We need one grid on dirichlet: the 4 corners are set to zero 67 67 md.stressbalance.spcvx = np.nan * np.ones((md.mesh.numberofvertices)) 68 68 md.stressbalance.spcvy = np.nan * np.ones((md.mesh.numberofvertices)) 69 69 md.stressbalance.spcvz = np.nan * np.ones((md.mesh.numberofvertices)) 70 pos = np.nonzero( logical_or.reduce_n(md.mesh.y == 0., md.mesh.x == 0., md.mesh.x == np.max(md.mesh.x), md.mesh.y == np.max(md.mesh.y)))#Don't take the same nodes two times71 md.stressbalance.spcvx[pos] = md.results.StressbalanceSolution.Vx[pos]72 md.stressbalance.spcvy[pos] = md.results.StressbalanceSolution.Vy[pos]70 pos = np.nonzero(np.logical_or.reduce((md.mesh.y == 0., md.mesh.x == 0., md.mesh.x == np.max(md.mesh.x), md.mesh.y == np.max(md.mesh.y)))) #Don't take the same nodes two times 71 md.stressbalance.spcvx[pos] = np.squeeze(md.results.StressbalanceSolution.Vx[pos]) 72 md.stressbalance.spcvy[pos] = np.squeeze(md.results.StressbalanceSolution.Vy[pos]) 73 73 md = setflowequation(md, 'FS', 'all') 74 74 md = solve(md, 'Stressbalance') 75 75 76 76 #Plot the results and save them 77 77 vx = md.results.StressbalanceSolution.Vx 78 78 vy = md.results.StressbalanceSolution.Vy … … 80 80 results.append(md.results.StressbalanceSolution) 81 81 82 # plotmodel(md, 'data', vx, 'data', vy, 'data', vz, 'layer #all', md.mesh.numberoflayers)82 # plotmodel(md, 'data', vx, 'data', vy, 'data', vz, 'layer #all', md.mesh.numberoflayers) 83 83 84 84 #Fields and tolerances to track changes
Note:
See TracChangeset
for help on using the changeset viewer.