| 1 | #Test Name: MISMIP3DHO
|
|---|
| 2 | import numpy as np
|
|---|
| 3 | from model import *
|
|---|
| 4 | from socket import gethostname
|
|---|
| 5 | from triangle import *
|
|---|
| 6 | from setmask import *
|
|---|
| 7 | from parameterize import *
|
|---|
| 8 | from setflowequation import *
|
|---|
| 9 | from solve import *
|
|---|
| 10 |
|
|---|
| 11 | md = triangle(model(), '../Exp/Square.exp', 100000.)
|
|---|
| 12 | md = setmask(md, '../Exp/SquareShelf.exp', '')
|
|---|
| 13 | md = parameterize(md, '../Par/SquareSheetShelf.py')
|
|---|
| 14 | md.initialization.vx[:] = 1.
|
|---|
| 15 | md.initialization.vy[:] = 1.
|
|---|
| 16 | md.geometry.thickness[:] = 500. - md.mesh.x / 10000.
|
|---|
| 17 | md.geometry.bed = -100. - md.mesh.x / 1000.
|
|---|
| 18 | md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water
|
|---|
| 19 | md.mask.ocean_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed
|
|---|
| 20 | pos = np.where(md.mask.ocean_levelset >= 0)
|
|---|
| 21 | md.geometry.base[pos] = md.geometry.bed[pos]
|
|---|
| 22 | md.geometry.surface = md.geometry.base + md.geometry.thickness
|
|---|
| 23 | md = md.extrude(4, 1.)
|
|---|
| 24 | md = setflowequation(md, 'HO', 'all')
|
|---|
| 25 |
|
|---|
| 26 | #Boundary conditions:
|
|---|
| 27 | md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices, ))
|
|---|
| 28 | md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0.
|
|---|
| 29 | md.stressbalance.spcvx[:] = float('Nan')
|
|---|
| 30 | md.stressbalance.spcvy[:] = float('Nan')
|
|---|
| 31 | md.stressbalance.spcvz[:] = float('Nan')
|
|---|
| 32 | posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)), np.array(np.where(md.mesh.y > 999999.9)))
|
|---|
| 33 | posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)), np.array(np.where(md.mesh.y > -0.1)))
|
|---|
| 34 | pos = np.unique(np.concatenate((posA, posB)))
|
|---|
| 35 | md.stressbalance.spcvy[pos] = 0.
|
|---|
| 36 | pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1)))
|
|---|
| 37 | md.stressbalance.spcvx[pos2] = 0.
|
|---|
| 38 | md.stressbalance.spcvy[pos2] = 0.
|
|---|
| 39 |
|
|---|
| 40 | md.materials.rheology_B = 1. / ((10**-25)**(1. / 3.)) * np.ones((md.mesh.numberofvertices, ))
|
|---|
| 41 | md.materials.rheology_law = 'None'
|
|---|
| 42 | md.friction.coefficient[:] = np.sqrt(1e7) * np.ones((md.mesh.numberofvertices, ))
|
|---|
| 43 | md.friction.p = 3. * np.ones((md.mesh.numberofelements, ))
|
|---|
| 44 | md.smb.mass_balance[:] = 1.
|
|---|
| 45 | md.basalforcings.groundedice_melting_rate[:] = 0.
|
|---|
| 46 | md.basalforcings.floatingice_melting_rate[:] = 30.
|
|---|
| 47 | md.transient.isthermal = 0
|
|---|
| 48 | md.transient.isstressbalance = 1
|
|---|
| 49 | md.transient.isgroundingline = 1
|
|---|
| 50 | md.transient.ismasstransport = 1
|
|---|
| 51 | md.transient.issmb = 1
|
|---|
| 52 | md.transient.requested_outputs = ['default', 'BasalforcingsFloatingiceMeltingRate']
|
|---|
| 53 | md.groundingline.migration = 'SubelementMigration'
|
|---|
| 54 | md.groundingline.melt_interpolation = 'SubelementMelt1'
|
|---|
| 55 | md.timestepping.final_time = 30
|
|---|
| 56 | md.timestepping.time_step = 10
|
|---|
| 57 |
|
|---|
| 58 | md.cluster = generic('name', gethostname(), 'np', 3)
|
|---|
| 59 | md = solve(md, 'Transient')
|
|---|
| 60 |
|
|---|
| 61 | #Fields and tolerances to track changes
|
|---|
| 62 | field_names = ['Bed1', 'Surface1', 'Thickness1', 'Floatingice1', 'Vx1', 'Vy1', 'Vz1', 'Pressure1', 'FloatingiceMeltingrate1',
|
|---|
| 63 | 'Bed2', 'Surface2', 'Thickness2', 'Floatingice2', 'Vx2', 'Vy2', 'Vz2', 'Pressure2', 'FloatingiceMeltingrate2',
|
|---|
| 64 | 'Bed3', 'Surface3', 'Thickness3', 'Floatingice3', 'Vx3', 'Vy3', 'Vz3', 'Pressure3', 'FloatingiceMeltingrate3']
|
|---|
| 65 | field_tolerances = [2e-11, 5e-12, 2e-11, 1e-11, 7e-10, 3e-08, 6e-10, 1e-13, 1e-13,
|
|---|
| 66 | 3e-11, 3e-11, 9e-10, 7e-11, 9e-09, 2e-07, 1e-09, 1e-10, 1e-13,
|
|---|
| 67 | 1e-9, 2e-08, 7e-09, 2e-7, 1e-03, 8e-04, 2e-09, 1e-10, 1e-13]
|
|---|
| 68 | field_values = [md.results.TransientSolution[0].Base,
|
|---|
| 69 | md.results.TransientSolution[0].Surface,
|
|---|
| 70 | md.results.TransientSolution[0].Thickness,
|
|---|
| 71 | md.results.TransientSolution[0].MaskOceanLevelset,
|
|---|
| 72 | md.results.TransientSolution[0].Vx,
|
|---|
| 73 | md.results.TransientSolution[0].Vy,
|
|---|
| 74 | md.results.TransientSolution[0].Vz,
|
|---|
| 75 | md.results.TransientSolution[0].Pressure,
|
|---|
| 76 | md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate,
|
|---|
| 77 | md.results.TransientSolution[1].Base,
|
|---|
| 78 | md.results.TransientSolution[1].Surface,
|
|---|
| 79 | md.results.TransientSolution[1].Thickness,
|
|---|
| 80 | md.results.TransientSolution[1].MaskOceanLevelset,
|
|---|
| 81 | md.results.TransientSolution[1].Vx,
|
|---|
| 82 | md.results.TransientSolution[1].Vy,
|
|---|
| 83 | md.results.TransientSolution[1].Vz,
|
|---|
| 84 | md.results.TransientSolution[1].Pressure,
|
|---|
| 85 | md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate,
|
|---|
| 86 | md.results.TransientSolution[2].Base,
|
|---|
| 87 | md.results.TransientSolution[2].Surface,
|
|---|
| 88 | md.results.TransientSolution[2].Thickness,
|
|---|
| 89 | md.results.TransientSolution[2].MaskOceanLevelset,
|
|---|
| 90 | md.results.TransientSolution[2].Vx,
|
|---|
| 91 | md.results.TransientSolution[2].Vy,
|
|---|
| 92 | md.results.TransientSolution[2].Vz,
|
|---|
| 93 | md.results.TransientSolution[2].Pressure,
|
|---|
| 94 | md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate]
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
