Changeset 21729 for issm/trunk/test/NightlyRun/test1207.py
- Timestamp:
- 05/19/17 14:48:02 (8 years ago)
- Location:
- issm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk
-
issm/trunk/test
- Property svn:mergeinfo changed
-
issm/trunk/test/NightlyRun/test1207.py
r21341 r21729 1 1 #Test Name: EISMINTRoundIceSheetStaticFS 2 import numpy 2 import numpy as np 3 3 from model import * 4 from socket import gethostname 4 5 from roundmesh import * 5 6 from setmask import * … … 7 8 from setflowequation import * 8 9 from solve import * 9 from MatlabFuncs import * 10 10 11 11 12 """ … … 26 27 #Calculation of the analytical 2d velocity field 27 28 constant=0.3 28 vx_obs=constant/2.*md.mesh.x .reshape(-1,1)*(md.geometry.thickness)**-129 vy_obs=constant/2.*md.mesh.y .reshape(-1,1)*(md.geometry.thickness)**-130 vel_obs=n umpy.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2)29 vx_obs=constant/2.*md.mesh.x*(md.geometry.thickness)**-1 30 vy_obs=constant/2.*md.mesh.y*(md.geometry.thickness)**-1 31 vel_obs=np.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2) 31 32 32 33 #We extrude the model to have a 3d model … … 35 36 36 37 #Spc the nodes on the bed 37 pos=n umpy.nonzero(md.mesh.vertexonbase)38 pos=np.where(md.mesh.vertexonbase) 38 39 md.stressbalance.spcvx[pos]=0. 39 40 md.stressbalance.spcvy[pos]=0. … … 41 42 42 43 #Now we can solve the problem 43 md.cluster=generic('name', oshostname(),'np',8)44 md.cluster=generic('name',gethostname(),'np',8) 44 45 md=solve(md,'Stressbalance') 45 46 … … 47 48 vx=md.results.StressbalanceSolution.Vx 48 49 vy=md.results.StressbalanceSolution.Vy 49 vel=n umpy.zeros((md.mesh.numberofvertices2d,1))50 vel=np.zeros((md.mesh.numberofvertices2d)) 50 51 51 52 for i in xrange(0,md.mesh.numberofvertices2d): … … 53 54 for j in xrange(0,md.mesh.numberoflayers-1): 54 55 node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*\ 55 (n umpy.sqrt(vx[i+(j+1)*md.mesh.numberofvertices2d,0]**2+vy[i+(j+1)*md.mesh.numberofvertices2d,0]**2)+\56 n umpy.sqrt(vx[i+j*md.mesh.numberofvertices2d,0]**2+vy[i+j*md.mesh.numberofvertices2d,0]**2))57 vel[i ,0]=node_vel56 (np.sqrt(vx[i+(j+1)*md.mesh.numberofvertices2d,0]**2+vy[i+(j+1)*md.mesh.numberofvertices2d,0]**2)+\ 57 np.sqrt(vx[i+j*md.mesh.numberofvertices2d,0]**2+vy[i+j*md.mesh.numberofvertices2d,0]**2)) 58 vel[i]=node_vel 58 59 59 60 #Plot of the velocity from the exact and calculated solutions … … 61 62 #subplot(2,2,1) 62 63 #p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 63 #vel,'FaceColor','interp','EdgeColor','none') ;64 #vel,'FaceColor','interp','EdgeColor','none') 64 65 #title('Modelled velocity','FontSize',14,'FontWeight','bold') 65 #colorbar ;66 #caxis([0 200]) ;66 #colorbar 67 #caxis([0 200]) 67 68 68 69 #subplot(2,2,2) 69 70 #p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 70 #vel_obs,'FaceColor','interp','EdgeColor','none') ;71 #vel_obs,'FaceColor','interp','EdgeColor','none') 71 72 #title('Analytical velocity','FontSize',14,'FontWeight','bold') 72 #colorbar ;73 #caxis([0 200]) ;73 #colorbar 74 #caxis([0 200]) 74 75 75 76 #subplot(2,2,3) 76 #hold on ;77 #plot(sqrt((md.mesh.x(1:md.mesh.numberofvertices2d)).^2+(md.mesh.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.') ;78 #plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.') ;79 #title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold') ;80 #xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold') ;81 #ylabel('velocity [m/yr]','FontSize',14,'FontWeight','bold') ;82 #legend('calculated velocity','exact velocity') ;83 #axis([0 750000 0 200]) ;84 #hold off ;77 #hold on 78 #plot(sqrt((md.mesh.x(1:md.mesh.numberofvertices2d)).^2+(md.mesh.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.') 79 #plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.') 80 #title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold') 81 #xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold') 82 #ylabel('velocity [m/yr]','FontSize',14,'FontWeight','bold') 83 #legend('calculated velocity','exact velocity') 84 #axis([0 750000 0 200]) 85 #hold off 85 86 86 87 #subplot(2,2,4) 87 88 #p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 88 #abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none') ;89 #abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none') 89 90 #title('Relative misfit [%]','FontSize',14,'FontWeight','bold') 90 #colorbar ;91 #caxis([0 100]) ;91 #colorbar 92 #caxis([0 100]) 92 93 93 94 if printingflag: 94 95 pass 95 96 # set(gcf,'Color','w') 96 # printmodel('FSstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off') ;97 # system(['mv FSstatic.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT/IceSheet']) ;97 # printmodel('FSstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off') 98 # system(['mv FSstatic.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT/IceSheet']) 98 99 99 100 #Fields and tolerances to track changes 100 field_names =[ \ 101 'Vx','Vy','Vel', \ 102 ] 103 field_tolerances=[ \ 104 1e-12,1e-12,1e-12, \ 105 ] 106 field_values=[ \ 107 vx,vy,vel, \ 108 ] 101 field_names =['Vx','Vy','Vel'] 102 field_tolerances=[1e-12,1e-12,1e-12] 103 field_values=[vx,vy,vel]
Note:
See TracChangeset
for help on using the changeset viewer.