![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/21724-22754/ISSM-22270-22271.diff@
22755
| Last change on this file since 22755 was 22755, checked in by , 7 years ago | |
|---|---|
| File size: 8.4 KB | |
-
../trunk-jpl/test/NightlyRun/test703.py
1 #Test Name: FlowbandFSsheetshelfTrans2 import numpy as np3 from scipy.interpolate import interp1d4 import copy5 from model import *6 from socket import gethostname7 from triangle import *8 from meshconvert import *9 from setmask import *10 from parameterize import *11 from setflowequation import *12 from solve import *13 from NowickiProfile import *14 from bamgflowband import *15 from paterson import *16 17 #mesh parameters18 x = np.arange(-5,5.5,.5).T19 [b,h,sea] = NowickiProfile(x)20 x = x * 10**321 h = h * 10**322 b = (b-sea) * 10**323 24 #mesh domain25 md = bamgflowband(model(),x,b+h,b,'hmax',150)26 print md.mesh.numberofvertices27 28 #geometry29 md.geometry.surface = interp1d(x,b+h)(md.mesh.x)30 md.geometry.base = interp1d(x,b)(md.mesh.x)31 md.geometry.thickness = md.geometry.surface - md.geometry.base32 33 #mask34 md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))35 md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 036 md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices,)) - 0.537 38 #materials39 md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices,))40 md.materials.rheology_B = paterson(md.initialization.temperature)41 md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))42 43 #friction44 md.friction.coefficient = np.zeros((md.mesh.numberofvertices,))45 md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 2046 md.friction.p = np.ones((md.mesh.numberofelements,))47 md.friction.q = np.ones((md.mesh.numberofelements,))48 49 #boundary conditions50 md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,))51 md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,))52 md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,))53 md.stressbalance.referential = float('NaN') * np.ones((md.mesh.numberofvertices,6))54 md.stressbalance.loadingforce = 0 * np.ones((md.mesh.numberofvertices,3))55 md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 800.56 md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0.57 58 #print md.stressbalance.referential59 #print md.stressbalance.loadingforce60 #print md.stressbalance.spcvx61 #print md.stressbalance.spcvy62 #print md.stressbalance.spcvz63 #print np.shape(md.stressbalance.referential)64 #print np.shape(md.stressbalance.loadingforce)65 #print np.shape(md.stressbalance.spcvx)66 #print np.shape(md.stressbalance.spcvy)67 #print np.shape(md.stressbalance.spcvz)68 69 #Misc70 md = setflowequation(md,'FS','all')71 md.flowequation.fe_FS = 'TaylorHood'72 md.stressbalance.abstol = float('NaN')73 md.miscellaneous.name = 'test703'74 75 #Transient settings76 md.timestepping.time_step = 0.00000177 md.timestepping.final_time = 0.00000578 md.stressbalance.shelf_dampening = 1.79 md.smb.mass_balance = np.zeros((md.mesh.numberofvertices,))80 md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,))81 md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices,))82 md.basalforcings.geothermalflux = np.zeros((md.mesh.numberofvertices,))83 posb = np.intersect1d(np.where(md.mesh.x > 0.), np.where(md.mesh.vertexonbase))84 md.basalforcings.groundedice_melting_rate[posb] = 18.85 md.basalforcings.floatingice_melting_rate[posb] = 18.86 md.initialization.vx = np.zeros((md.mesh.numberofvertices,))87 md.initialization.vy = np.zeros((md.mesh.numberofvertices,))88 md.initialization.pressure = np.zeros((md.mesh.numberofvertices,))89 md.masstransport.spcthickness = float('NaN') * np.ones((md.mesh.numberofvertices,))90 md.thermal.spctemperature = float('NaN') * np.ones((md.mesh.numberofvertices,))91 md.transient.isthermal = 092 md.masstransport.isfreesurface = 193 94 #Go solve95 md.cluster = generic('np',3)96 md.stressbalance.shelf_dampening = 197 md1 = solve(md,'Transient')98 99 md.stressbalance.shelf_dampening = 0100 md = solve(md,'Transient')101 102 #Fields and tolerances to track changes103 field_names = [104 'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',105 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',106 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3',107 'Vx1_damp','Vy1_damp','Vel1_damp','Pressure1_damp','Bed1_damp','Surface1_damp','Thickness1_damp',108 'Vx2_damp','Vy2_damp','Vel2_damp','Pressure2_damp','Bed2_damp','Surface2_damp','Thickness2_damp',109 'Vx3_damp','Vy3_damp','Vel3_damp','Pressure3_damp','Bed3_damp','Surface3_damp','Thickness3_damp']110 field_tolerances = [111 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,112 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,113 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,114 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,115 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,116 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10]117 field_values = [118 md.results.TransientSolution[0].Vx,119 md.results.TransientSolution[0].Vy,120 md.results.TransientSolution[0].Vel,121 md.results.TransientSolution[0].Pressure,122 md.results.TransientSolution[0].Base,123 md.results.TransientSolution[0].Surface,124 md.results.TransientSolution[0].Thickness,125 md.results.TransientSolution[1].Vx,126 md.results.TransientSolution[1].Vy,127 md.results.TransientSolution[1].Vel,128 md.results.TransientSolution[1].Pressure,129 md.results.TransientSolution[1].Base,130 md.results.TransientSolution[1].Surface,131 md.results.TransientSolution[1].Thickness,132 md.results.TransientSolution[2].Vx,133 md.results.TransientSolution[2].Vy,134 md.results.TransientSolution[2].Vel,135 md.results.TransientSolution[2].Pressure,136 md.results.TransientSolution[2].Base,137 md.results.TransientSolution[2].Surface,138 md.results.TransientSolution[2].Thickness,139 md1.results.TransientSolution[0].Vx,140 md1.results.TransientSolution[0].Vy,141 md1.results.TransientSolution[0].Vel,142 md1.results.TransientSolution[0].Pressure,143 md1.results.TransientSolution[0].Base,144 md1.results.TransientSolution[0].Surface,145 md1.results.TransientSolution[0].Thickness,146 md1.results.TransientSolution[1].Vx,147 md1.results.TransientSolution[1].Vy,148 md1.results.TransientSolution[1].Vel,149 md1.results.TransientSolution[1].Pressure,150 md1.results.TransientSolution[1].Base,151 md1.results.TransientSolution[1].Surface,152 md1.results.TransientSolution[1].Thickness,153 md1.results.TransientSolution[2].Vx,154 md1.results.TransientSolution[2].Vy,155 md1.results.TransientSolution[2].Vel,156 md1.results.TransientSolution[2].Pressure,157 md1.results.TransientSolution[2].Base,158 md1.results.TransientSolution[2].Surface,159 md1.results.TransientSolution[2].Thickness,160 ] -
../trunk-jpl/src/m/mesh/bamg.py
1 1 import os.path 2 2 import numpy as np 3 from mesh2d import * 3 4 from mesh2dvertical import * 5 from mesh3dsurface import * 4 6 from collections import OrderedDict 5 7 from pairoptions import pairoptions 6 8 from bamggeom import bamggeom … … 91 93 for i,domaini in enumerate(domain): 92 94 93 95 #Check that the domain is closed 94 if (domaini .x[0]!=domaini.x[-1] or domaini.y[0]!=domaini.y[-1]):96 if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]): 95 97 raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed") 96 98 97 99 #Checks that all holes are INSIDE the principle domain outline 98 100 if i: 99 flags=ContourToNodes(domaini .x,domaini.y,domainfile,0)[0]101 flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0] 100 102 if np.any(np.logical_not(flags)): 101 103 raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") 102 104 103 105 #Add all points to bamg_geometry 104 nods=domaini .nods-1 #the domain are closed 0=end105 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini .x[0:nods],domaini.y[0:nods],np.ones((nods)))).T))106 nods=domaini['nods']-1 #the domain are closed 0=end 107 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T)) 106 108 bamg_geometry.Edges =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T)) 107 109 if i: 108 110 bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,-subdomain_ref]))
Note:
See TracBrowser
for help on using the repository browser.
