![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 23793 for issm/trunk-jpl/test/NightlyRun/test1107.py
- Timestamp:
- 03/13/19 03:17:46 (6 years ago)
- File:
-
- 1 edited
-
issm/trunk-jpl/test/NightlyRun/test1107.py (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/test/NightlyRun/test1107.py
r21411 r23793 1 1 #Test Name: ISMIPDHO 2 2 import numpy as np 3 import shutil4 3 from model import * 5 4 from socket import gethostname … … 15 14 """ 16 15 17 printingflag =False16 printingflag = False 18 17 19 L_list =[80000.]20 results =[]21 minvx =[]22 maxvx =[]18 L_list = [80000.] 19 results = [] 20 minvx = [] 21 maxvx = [] 23 22 24 23 for L in L_list: 25 nx=30 #numberof nodes in x direction26 ny=3027 md=model()28 md=squaremesh(md,L,L,nx,ny)29 md=setmask(md,'','') #ice sheet test30 md=parameterize(md,'../Par/ISMIPD.py')31 md.extrude(10,1.)24 nx = 30 #numberof nodes in x direction 25 ny = 30 26 md = model() 27 md = squaremesh(md, L, L, nx, ny) 28 md = setmask(md, '', '') #ice sheet test 29 md = parameterize(md, '../Par/ISMIPD.py') 30 md.extrude(10, 1.) 32 31 33 md=setflowequation(md,'HO','all')32 md = setflowequation(md, 'HO', 'all') 34 33 35 36 md.stressbalance.spcvx=np.nan*np.ones((md.mesh.numberofvertices))37 md.stressbalance.spcvy=np.nan*np.ones((md.mesh.numberofvertices))38 md.stressbalance.spcvz=np.nan*np.ones((md.mesh.numberofvertices))34 #We need one grd on dirichlet: the 4 corners are set to zero 35 md.stressbalance.spcvx = np.nan * np.ones((md.mesh.numberofvertices)) 36 md.stressbalance.spcvy = np.nan * np.ones((md.mesh.numberofvertices)) 37 md.stressbalance.spcvz = np.nan * np.ones((md.mesh.numberofvertices)) 39 38 40 41 # posx=find(md.mesh.x==0. & ~(md.mesh.y==0. & md.mesh.vertexonbase) & ~(md.mesh.y==L & md.mesh.vertexonbase))42 posx=np.where(np.logical_and.reduce((md.mesh.x==0.,np.logical_not(np.logical_and(md.mesh.y==0.,md.mesh.vertexonbase)),np.logical_not(np.logical_and(md.mesh.y==L,md.mesh.vertexonbase)))))[0]43 # posx2=find(md.mesh.x==max(md.mesh.x) & ~(md.mesh.y==0. & md.mesh.vertexonbase) & ~(md.mesh.y==L & md.mesh.vertexonbase))44 posx2=np.where(np.logical_and.reduce((md.mesh.x==np.max(md.mesh.x),np.logical_not(np.logical_and(md.mesh.y==0.,md.mesh.vertexonbase)),np.logical_not(np.logical_and(md.mesh.y==L,md.mesh.vertexonbase)))))[0]39 #Create MPCs to have periodic boundary conditions 40 # posx = find(md.mesh.x = 0. & ~(md.mesh.y = 0. & md.mesh.vertexonbase) & ~(md.mesh.y = L & md.mesh.vertexonbase)) 41 posx = np.where(np.logical_and.reduce((md.mesh.x == 0., np.logical_not(np.logical_and(md.mesh.y == 0., md.mesh.vertexonbase)), np.logical_not(np.logical_and(md.mesh.y == L, md.mesh.vertexonbase)))))[0] 42 # posx2 = find(md.mesh.x = max(md.mesh.x) & ~(md.mesh.y = 0. & md.mesh.vertexonbase) & ~(md.mesh.y = L & md.mesh.vertexonbase)) 43 posx2 = np.where(np.logical_and.reduce((md.mesh.x == np.max(md.mesh.x), np.logical_not(np.logical_and(md.mesh.y == 0., md.mesh.vertexonbase)), np.logical_not(np.logical_and(md.mesh.y == L, md.mesh.vertexonbase)))))[0] 45 44 46 posy=np.where(np.logical_and.reduce((md.mesh.y==0.,md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x))))[0] #Don't take the same nodes two times47 posy2=np.where(np.logical_and.reduce((md.mesh.y==np.max(md.mesh.y),md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x))))[0]45 posy = np.where(np.logical_and.reduce((md.mesh.y == 0., md.mesh.x != 0., md.mesh.x != np.max(md.mesh.x))))[0] #Don't take the same nodes two times 46 posy2 = np.where(np.logical_and.reduce((md.mesh.y == np.max(md.mesh.y), md.mesh.x != 0., md.mesh.x != np.max(md.mesh.x))))[0] 48 47 49 md.stressbalance.vertex_pairing=np.vstack((np.vstack((posx+1,posx2+1)).T,np.vstack((posy+1,posy2+1)).T))48 md.stressbalance.vertex_pairing = np.vstack((np.vstack((posx + 1, posx2 + 1)).T, np.vstack((posy + 1, posy2 + 1)).T)) 50 49 51 #Add spc on the corners 52 pos=np.where(np.logical_and.reduce((np.logical_or(md.mesh.x==0.,md.mesh.x==L),np.logical_or(md.mesh.y==0.,md.mesh.y==L),md.mesh.vertexonbase))) 53 md.stressbalance.spcvy[:]=0. 54 md.stressbalance.spcvx[pos]=0. 55 if (L==5000.): 56 md.stressbalance.spcvx[pos]=16.0912 57 elif (L==10000.): 58 md.stressbalance.spcvx[pos]=16.52 59 elif (L==20000.): 60 md.stressbalance.spcvx[pos]=17.77 61 elif (L==40000.): 62 md.stressbalance.spcvx[pos]=19.88 63 elif (L==80000.): 64 md.stressbalance.spcvx[pos]=18.65 65 elif (L==160000.): 66 md.stressbalance.spcvx[pos]=16.91 67 68 #Spc the bed at zero for vz 69 pos=np.where(md.mesh.vertexonbase) 70 md.stressbalance.spcvz[pos]=0. 50 #Add spc on the corners 51 pos = np.where(np.logical_and.reduce((np.logical_or(md.mesh.x == 0., md.mesh.x == L), np.logical_or(md.mesh.y == 0., md.mesh.y == L), md.mesh.vertexonbase))) 52 md.stressbalance.spcvy[:] = 0. 53 md.stressbalance.spcvx[pos] = 0. 54 if (L == 5000.): 55 md.stressbalance.spcvx[pos] = 16.0912 56 elif (L == 10000.): 57 md.stressbalance.spcvx[pos] = 16.52 58 elif (L == 20000.): 59 md.stressbalance.spcvx[pos] = 17.77 60 elif (L == 40000.): 61 md.stressbalance.spcvx[pos] = 19.88 62 elif (L == 80000.): 63 md.stressbalance.spcvx[pos] = 18.65 64 elif (L == 160000.): 65 md.stressbalance.spcvx[pos] = 16.91 71 66 72 #Compute the stressbalance 73 md.cluster=generic('name',gethostname(),'np',8)74 md=solve(md,'Stressbalance') 67 #Spc the bed at zero for vz 68 pos = np.where(md.mesh.vertexonbase) 69 md.stressbalance.spcvz[pos] = 0. 75 70 76 #Plot the results and save them 77 vx=md.results.StressbalanceSolution.Vx 78 vy=md.results.StressbalanceSolution.Vy 79 vz=md.results.StressbalanceSolution.Vz 80 results.append(md.results.StressbalanceSolution) 81 minvx.append(np.min(vx[-md.mesh.numberofvertices2d:])) 82 maxvx.append(np.max(vx[-md.mesh.numberofvertices2d:])) 71 #Compute the stressbalance 72 md.cluster = generic('name', gethostname(), 'np', 8) 73 md = solve(md, 'Stressbalance') 83 74 84 #Now plot vx, vy, vz and vx on a cross section 85 # plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2) 86 if printingflag: 87 pass 88 # set(gcf,'Color','w') 89 # printmodel(['ismipdHOvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off') 90 # shutil.move("ismipdHOvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 91 # plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3) 92 if printingflag: 93 pass 94 # set(gcf,'Color','w') 95 # printmodel(['ismipdHOvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off') 96 # shutil.move("ismipdHOvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 75 #Plot the results and save them 76 vx = md.results.StressbalanceSolution.Vx 77 vy = md.results.StressbalanceSolution.Vy 78 vz = md.results.StressbalanceSolution.Vz 79 results.append(md.results.StressbalanceSolution) 80 minvx.append(np.min(vx[-md.mesh.numberofvertices2d:])) 81 maxvx.append(np.max(vx[-md.mesh.numberofvertices2d:])) 97 82 98 if (L==5000.): 99 pass 100 # plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,... 101 # 'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',4) 102 elif (L==10000.): 103 pass 104 # plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,... 105 # 'resolution',[10 10],'ylim',[0 20],'xlim',[0 10000],'title','','xlabel','','figure',4) 106 elif (L==20000.): 107 pass 108 # plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,... 109 # 'resolution',[10 10],'ylim',[0 30],'xlim',[0 20000],'title','','xlabel','','figure',4) 110 elif (L==40000.): 111 pass 112 # plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,... 113 # 'resolution',[10 10],'ylim',[10 60],'xlim',[0 40000],'title','','xlabel','','figure',4) 114 elif (L==80000.): 115 pass 116 # plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,... 117 # 'resolution',[10 10],'ylim',[0 200],'xlim',[0 80000],'title','','xlabel','','figure',4) 118 elif (L==160000.): 119 pass 120 # plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,... 121 # 'resolution',[10 10],'ylim',[0 400],'xlim',[0 160000],'title','','xlabel','','figure',4) 122 if printingflag: 123 pass 124 # set(gcf,'Color','w') 125 # printmodel(['ismipdHOvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off') 126 # shutil.move("ismipdHOvxsec%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 83 #Now plot vx, vy, vz and vx on a cross section 84 # plotmodel(md, 'data', vx, 'layer#all', md.mesh.numberoflayers, 'xlim',[0 L/10^3], 'ylim',[0 L/10^3], 'unit', 'km', 'figure', 2) 85 if printingflag: 86 pass 87 # set(gcf, 'Color', 'w') 88 # printmodel(['ismipdHOvx' num2str(L)], 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off') 89 # shutil.move("ismipdHOvx%d.png" % L, ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 90 # plotmodel(md, 'data', vz, 'layer#all', md.mesh.numberoflayers, 'xlim',[0 L/10^3], 'ylim',[0 L/10^3], 'unit', 'km', 'figure', 3) 91 if printingflag: 92 pass 93 # set(gcf, 'Color', 'w') 94 # printmodel(['ismipdHOvz' num2str(L)], 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off') 95 # shutil.move("ismipdHOvz%d.png" % L, ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 96 97 if (L == 5000.): 98 pass 99 # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP5000.exp', 'layer', md.mesh.numberoflayers,... 100 # 'resolution',[10 10], 'ylim',[0 20], 'xlim',[0 5000], 'title', '', 'xlabel', '', 'figure', 4) 101 elif (L == 10000.): 102 pass 103 # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP10000.exp', 'layer', md.mesh.numberoflayers,... 104 # 'resolution',[10 10], 'ylim',[0 20], 'xlim',[0 10000], 'title', '', 'xlabel', '', 'figure', 4) 105 elif (L == 20000.): 106 pass 107 # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP20000.exp', 'layer', md.mesh.numberoflayers,... 108 # 'resolution',[10 10], 'ylim',[0 30], 'xlim',[0 20000], 'title', '', 'xlabel', '', 'figure', 4) 109 elif (L == 40000.): 110 pass 111 # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP40000.exp', 'layer', md.mesh.numberoflayers,... 112 # 'resolution',[10 10], 'ylim',[10 60], 'xlim',[0 40000], 'title', '', 'xlabel', '', 'figure', 4) 113 elif (L == 80000.): 114 pass 115 # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP80000.exp', 'layer', md.mesh.numberoflayers,... 116 # 'resolution',[10 10], 'ylim',[0 200], 'xlim',[0 80000], 'title', '', 'xlabel', '', 'figure', 4) 117 elif (L == 160000.): 118 pass 119 # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP160000.exp', 'layer', md.mesh.numberoflayers,... 120 # 'resolution',[10 10], 'ylim',[0 400], 'xlim',[0 160000], 'title', '', 'xlabel', '', 'figure', 4) 121 if printingflag: 122 pass 123 # set(gcf, 'Color', 'w') 124 # printmodel(['ismipdHOvxsec' num2str(L)], 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off') 125 # shutil.move("ismipdHOvxsec%d.png" % L, ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 127 126 128 127 #Now plot the min and max values of vx for each size of the square 129 #plot([5 10 20 40 80 160], minvx)ylim([2 18])xlim([0 160])128 #plot([5 10 20 40 80 160], minvx)ylim([2 18])xlim([0 160]) 130 129 if printingflag: 131 132 # set(gcf,'Color','w')133 # printmodel('ismipdHOminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')134 # shutil.move('ismipdHOminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD')135 #plot([5 10 20 40 80 160], maxvx)ylim([0 300])xlim([0 160])130 pass 131 # set(gcf, 'Color', 'w') 132 # printmodel('ismipdHOminvx', 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off') 133 # shutil.move('ismipdHOminvx.png', ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 134 #plot([5 10 20 40 80 160], maxvx)ylim([0 300])xlim([0 160]) 136 135 if printingflag: 137 138 # set(gcf,'Color','w')139 # printmodel('ismipdHOmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')140 # shutil.move('ismipdHOmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD')136 pass 137 # set(gcf, 'Color', 'w') 138 # printmodel('ismipdHOmaxvx', 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off') 139 # shutil.move('ismipdHOmaxvx.png', ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD') 141 140 142 141 #Fields and tolerances to track changes 143 field_names =['Vx80km','Vy80km','Vz80km']144 field_tolerances =[1e-08,1e-08,1e-07]145 field_values =[]142 field_names = ['Vx80km', 'Vy80km', 'Vz80km'] 143 field_tolerances = [1e-08, 1e-08, 1e-07] 144 field_values = [] 146 145 for result in results: 147 field_values=field_values+[result.Vx,result.Vy,result.Vz]146 field_values = field_values + [result.Vx, result.Vy, result.Vz]
Note:
See TracChangeset
for help on using the changeset viewer.
