Ignore:
Timestamp:
05/19/17 14:48:02 (8 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 21727

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk

  • issm/trunk/test

  • issm/trunk/test/NightlyRun/test1207.py

    r21341 r21729  
    11#Test Name: EISMINTRoundIceSheetStaticFS
    2 import numpy
     2import numpy as np
    33from model import *
     4from socket import gethostname
    45from roundmesh import *
    56from setmask import *
     
    78from setflowequation import *
    89from solve import *
    9 from MatlabFuncs import *
     10
    1011
    1112"""
     
    2627#Calculation of the analytical 2d velocity field
    2728constant=0.3
    28 vx_obs=constant/2.*md.mesh.x.reshape(-1,1)*(md.geometry.thickness)**-1
    29 vy_obs=constant/2.*md.mesh.y.reshape(-1,1)*(md.geometry.thickness)**-1
    30 vel_obs=numpy.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2)
     29vx_obs=constant/2.*md.mesh.x*(md.geometry.thickness)**-1
     30vy_obs=constant/2.*md.mesh.y*(md.geometry.thickness)**-1
     31vel_obs=np.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2)
    3132
    3233#We extrude the model to have a 3d model
     
    3536
    3637#Spc the nodes on the bed
    37 pos=numpy.nonzero(md.mesh.vertexonbase)
     38pos=np.where(md.mesh.vertexonbase)
    3839md.stressbalance.spcvx[pos]=0.
    3940md.stressbalance.spcvy[pos]=0.
     
    4142
    4243#Now we can solve the problem
    43 md.cluster=generic('name',oshostname(),'np',8)
     44md.cluster=generic('name',gethostname(),'np',8)
    4445md=solve(md,'Stressbalance')
    4546
     
    4748vx=md.results.StressbalanceSolution.Vx
    4849vy=md.results.StressbalanceSolution.Vy
    49 vel=numpy.zeros((md.mesh.numberofvertices2d,1))
     50vel=np.zeros((md.mesh.numberofvertices2d))
    5051
    5152for i in xrange(0,md.mesh.numberofvertices2d):
     
    5354        for j in xrange(0,md.mesh.numberoflayers-1):
    5455                node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*\
    55                         (numpy.sqrt(vx[i+(j+1)*md.mesh.numberofvertices2d,0]**2+vy[i+(j+1)*md.mesh.numberofvertices2d,0]**2)+\
    56                         numpy.sqrt(vx[i+j*md.mesh.numberofvertices2d,0]**2+vy[i+j*md.mesh.numberofvertices2d,0]**2))
    57         vel[i,0]=node_vel
     56                        (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
    5859
    5960#Plot of the velocity from the exact and calculated solutions
     
    6162#subplot(2,2,1)
    6263#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')
    6465#title('Modelled velocity','FontSize',14,'FontWeight','bold')
    65 #colorbar;
    66 #caxis([0 200]);
     66#colorbar
     67#caxis([0 200])
    6768   
    6869#subplot(2,2,2)
    6970#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')
    7172#title('Analytical velocity','FontSize',14,'FontWeight','bold')
    72 #colorbar;
    73 #caxis([0 200]);
     73#colorbar
     74#caxis([0 200])
    7475
    7576#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
    8586
    8687#subplot(2,2,4)
    8788#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')
    8990#title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
    90 #colorbar;
    91 #caxis([0 100]);
     91#colorbar
     92#caxis([0 100])
    9293
    9394if printingflag:
    9495        pass
    9596#       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'])
    9899
    99100#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 ]
     101field_names     =['Vx','Vy','Vel']
     102field_tolerances=[1e-12,1e-12,1e-12]
     103field_values=[vx,vy,vel]
Note: See TracChangeset for help on using the changeset viewer.