| 1 | #Test Name: ISMIPAHO
|
|---|
| 2 | import numpy as np
|
|---|
| 3 | from model import *
|
|---|
| 4 | from socket import gethostname
|
|---|
| 5 | from squaremesh import *
|
|---|
| 6 | from setmask import *
|
|---|
| 7 | from parameterize import *
|
|---|
| 8 | from setflowequation import *
|
|---|
| 9 | from solve import *
|
|---|
| 10 |
|
|---|
| 11 | """
|
|---|
| 12 | This test is a test from the ISMP - HOM Intercomparison project.
|
|---|
| 13 | Pattyn and Payne 2006
|
|---|
| 14 | """
|
|---|
| 15 |
|
|---|
| 16 | printingflag = False
|
|---|
| 17 |
|
|---|
| 18 | L_list = [80000.]
|
|---|
| 19 | results = []
|
|---|
| 20 | minvx = []
|
|---|
| 21 | maxvx = []
|
|---|
| 22 |
|
|---|
| 23 | for L in L_list:
|
|---|
| 24 | nx = 20 #numberof nodes in x direction
|
|---|
| 25 | ny = 20
|
|---|
| 26 | md = model()
|
|---|
| 27 | md = squaremesh(md, L, L, nx, ny)
|
|---|
| 28 | md = setmask(md, '', '') #ice sheet test
|
|---|
| 29 | md = parameterize(md, '../Par/ISMIPA.py')
|
|---|
| 30 | md.extrude(9, 1.)
|
|---|
| 31 |
|
|---|
| 32 | md = setflowequation(md, 'HO', 'all')
|
|---|
| 33 |
|
|---|
| 34 | #Create dirichlet on the bed only
|
|---|
| 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))
|
|---|
| 38 |
|
|---|
| 39 | pos = np.where(md.mesh.vertexonbase)
|
|---|
| 40 | md.stressbalance.spcvx[pos] = 0.
|
|---|
| 41 | md.stressbalance.spcvy[pos] = 0.
|
|---|
| 42 |
|
|---|
| 43 | #Create MPCs to have periodic boundary conditions
|
|---|
| 44 | posx = np.where(md.mesh.x == 0.)[0]
|
|---|
| 45 | posx2 = np.where(md.mesh.x == np.max(md.mesh.x))[0]
|
|---|
| 46 |
|
|---|
| 47 | 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
|
|---|
| 48 | 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]
|
|---|
| 49 |
|
|---|
| 50 | md.stressbalance.vertex_pairing = np.vstack((np.vstack((posx + 1, posx2 + 1)).T, np.vstack((posy + 1, posy2 + 1)).T))
|
|---|
| 51 |
|
|---|
| 52 | #Compute the stressbalance
|
|---|
| 53 | md.cluster = generic('name', gethostname(), 'np', 8)
|
|---|
| 54 | md = solve(md, 'Stressbalance')
|
|---|
| 55 |
|
|---|
| 56 | #Plot the results and save them
|
|---|
| 57 | vx = md.results.StressbalanceSolution.Vx
|
|---|
| 58 | vy = md.results.StressbalanceSolution.Vy
|
|---|
| 59 | vz = md.results.StressbalanceSolution.Vz
|
|---|
| 60 | results.append(md.results.StressbalanceSolution)
|
|---|
| 61 | minvx.append(np.min(vx[-md.mesh.numberofvertices2d:]))
|
|---|
| 62 | maxvx.append(np.max(vx[-md.mesh.numberofvertices2d:]))
|
|---|
| 63 |
|
|---|
| 64 | #Now plot vx, vy, vz and vx on a cross section
|
|---|
| 65 | # plotmodel(md, 'data', vx, 'layer #all', md.mesh.numberoflayers, 'xlim', [0 L / 1.0e3], 'ylim', [0 L / 1.0e3], 'unit', 'km')
|
|---|
| 66 | if printingflag:
|
|---|
| 67 | pass
|
|---|
| 68 | # set(gcf, 'Color', 'w')
|
|---|
| 69 | # printmodel(['ismipaHOvx' num2str(L)], 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off')
|
|---|
| 70 | # shutil.move("ismipaHOvx%d.png" % L, ISSM_DIR + '/website/doc_pdf/validation/Images/ISMIP/TestA')
|
|---|
| 71 | # plotmodel(md, 'data', vy, 'layer #all', md.mesh.numberoflayers, 'xlim', [0 L / 1.0e3], 'ylim', [0 L / 1.0e3], 'unit', 'km')
|
|---|
| 72 | if printingflag:
|
|---|
| 73 | pass
|
|---|
| 74 | # set(gcf, 'Color', 'w')
|
|---|
| 75 | # printmodel(['ismipaHOvy' num2str(L)], 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off')
|
|---|
| 76 | # shutil.move("ismipaHOvy%d.png" % L, ISSM_DIR + '/website/doc_pdf/validation/Images/ISMIP/TestA')
|
|---|
| 77 | # plotmodel(md, 'data', vz, 'layer #all', md.mesh.numberoflayers, 'xlim', [0 L / 1.0e3], 'ylim', [0 L / 1.0e3], 'unit', 'km')
|
|---|
| 78 | if printingflag:
|
|---|
| 79 | pass
|
|---|
| 80 | # set(gcf, 'Color', 'w')
|
|---|
| 81 | # printmodel(['ismipaHOvz' num2str(L)], 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off')
|
|---|
| 82 | # shutil.move("ismipaHOvz%d.png" % L, ISSM_DIR + '/website/doc_pdf/validation/Images/ISMIP/TestA')
|
|---|
| 83 |
|
|---|
| 84 | if (L == 5000.):
|
|---|
| 85 | pass
|
|---|
| 86 | # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP5000.exp', 'layer', md.mesh.numberoflayers, ...
|
|---|
| 87 | # 'resolution', [10 10], 'ylim', [10 18], 'xlim', [0 5000], 'title', '', 'xlabel', '')
|
|---|
| 88 | elif (L == 10000.):
|
|---|
| 89 | pass
|
|---|
| 90 | # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP10000.exp', 'layer', md.mesh.numberoflayers, ...
|
|---|
| 91 | # 'resolution', [10 10], 'ylim', [10 30], 'xlim', [0 10000], 'title', '', 'xlabel', '')
|
|---|
| 92 | elif (L == 20000.):
|
|---|
| 93 | pass
|
|---|
| 94 | # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP20000.exp', 'layer', md.mesh.numberoflayers, ...
|
|---|
| 95 | # 'resolution', [10 10], 'ylim', [0 50], 'xlim', [0 20000], 'title', '', 'xlabel', '')
|
|---|
| 96 | elif (L == 40000.):
|
|---|
| 97 | pass
|
|---|
| 98 | # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP40000.exp', 'layer', md.mesh.numberoflayers, ...
|
|---|
| 99 | # 'resolution', [10 10], 'ylim', [0 80], 'xlim', [0 40000], 'title', '', 'xlabel', '')
|
|---|
| 100 | elif (L == 80000.):
|
|---|
| 101 | pass
|
|---|
| 102 | # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP80000.exp', 'layer', md.mesh.numberoflayers, ...
|
|---|
| 103 | # 'resolution', [10 10], 'ylim', [0 100], 'xlim', [0 80000], 'title', '', 'xlabel', '')
|
|---|
| 104 | elif (L == 160000.):
|
|---|
| 105 | pass
|
|---|
| 106 | # plotmodel(md, 'data', vx, 'sectionvalue', '../Exp/ISMIP160000.exp', 'layer', md.mesh.numberoflayers, ...
|
|---|
| 107 | # 'resolution', [10 10], 'ylim', [0 120], 'xlim', [0 160000], 'title', '', 'xlabel', '')
|
|---|
| 108 | if printingflag:
|
|---|
| 109 | pass
|
|---|
| 110 | # set(gcf, 'Color', 'w')
|
|---|
| 111 | # printmodel(['ismipaHOvxsec' num2str(L)], 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off')
|
|---|
| 112 | # shutil.move("ismipaHOvxsec%d.png" % L, ISSM_DIR + '/website/doc_pdf/validation/Images/ISMIP/TestA')
|
|---|
| 113 |
|
|---|
| 114 | #Now plot the min and max values of vx for each size of the square
|
|---|
| 115 | #plot([5 10 20 40 80 160], minvx)ylim([0 18])xlim([0 160])
|
|---|
| 116 | if printingflag:
|
|---|
| 117 | pass
|
|---|
| 118 | # set(gcf, 'Color', 'w')
|
|---|
| 119 | # printmodel('ismipaHOminvx', 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off')
|
|---|
| 120 | # shutil.move('ismipaHOminvx.png', ISSM_DIR + '/website/doc_pdf/validation/Images/ISMIP/TestA')
|
|---|
| 121 | #plot([5 10 20 40 80 160], maxvx)ylim([0 120])xlim([0 160])
|
|---|
| 122 | if printingflag:
|
|---|
| 123 | pass
|
|---|
| 124 | # set(gcf, 'Color', 'w')
|
|---|
| 125 | # printmodel('ismipaHOmaxvx', 'png', 'margin', 'on', 'marginsize', 25, 'frame', 'off', 'resolution', 1.5, 'hardcopy', 'off')
|
|---|
| 126 | # shutil.move('ismipaHOmaxvx.png', ISSM_DIR + '/website/doc_pdf/validation/Images/ISMIP/TestA')
|
|---|
| 127 |
|
|---|
| 128 | #Fields and tolerances to track changes
|
|---|
| 129 | field_names = ['Vx80km', 'Vy80km', 'Vz80km']
|
|---|
| 130 | field_tolerances = [1e-08, 1e-08, 1e-08]
|
|---|
| 131 | field_values = []
|
|---|
| 132 | for result in results:
|
|---|
| 133 | field_values = field_values + [result.Vx, result.Vy, result.Vz]
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
