Index: /issm/trunk-jpl/test/NightlyRun/test248.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test248.m (revision 26438) +++ /issm/trunk-jpl/test/NightlyRun/test248.m (revision 26438) @@ -0,0 +1,24 @@ +%Test Name: SquareShelfStressMLHO2d +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md=setflowequation(md,'MLHO','all'); +md.cluster=generic('name',oshostname(),'np',3); +md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'}; +md=solve(md,'Stressbalance'); + +%Fields and tolerances to track changes +field_names ={'Vx','Vy','Vel','Pressure','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.StressbalanceSolution.Vx),... + (md.results.StressbalanceSolution.Vy),... + (md.results.StressbalanceSolution.Vel),... + (md.results.StressbalanceSolution.Pressure),... + (md.results.StressbalanceSolution.VxSurface),... + (md.results.StressbalanceSolution.VySurface),... + (md.results.StressbalanceSolution.VxShear),... + (md.results.StressbalanceSolution.VyShear),... + (md.results.StressbalanceSolution.VxBase),... + (md.results.StressbalanceSolution.VyBase),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test248.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test248.py (revision 26438) +++ /issm/trunk-jpl/test/NightlyRun/test248.py (revision 26438) @@ -0,0 +1,32 @@ +#Test Name: SquareShelfStressMLHO2d +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +from ContourToMesh import * + +md = triangle(model(), '../Exp/Square.exp', 150000) +md = setmask(md, 'all', '') +md = parameterize(md, '../Par/SquareShelf.py') +md = setflowequation(md, 'MLHO', 'all') +md.cluster = generic('name', gethostname(), 'np', 3) +md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase'] +md = solve(md, 'Stressbalance') + +#Fields and tolerances to track changes +field_names = ['Vx', 'Vy', 'Vel', 'Pressure', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase'] +field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13] +field_values = [md.results.StressbalanceSolution.Vx, + md.results.StressbalanceSolution.Vy, + md.results.StressbalanceSolution.Vel, + md.results.StressbalanceSolution.Pressure, + md.results.StressbalanceSolution.VxSurface, + md.results.StressbalanceSolution.VySurface, + md.results.StressbalanceSolution.VxShear, + md.results.StressbalanceSolution.VyShear, + md.results.StressbalanceSolution.VxBase, + md.results.StressbalanceSolution.VyBase] Index: /issm/trunk-jpl/test/NightlyRun/test249.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test249.m (revision 26438) +++ /issm/trunk-jpl/test/NightlyRun/test249.m (revision 26438) @@ -0,0 +1,46 @@ +%Test Name: SquareShelfTranMLHO2d +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md=setflowequation(md,'MLHO','all'); +md.cluster=generic('name',oshostname(),'np',3); +md.transient.requested_outputs={'default','FloatingArea','GroundedArea','TotalGroundedBmb','TotalFloatingBmb'}; +md.basalforcings.floatingice_melting_rate(:)=1; +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','TotalGroundedBmb1','TotalFloatingBmb1',... + 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','TotalGroundedBmb2','TotalFloatingBmb2',... + 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','TotalGroundedBmb3','TotalFloatingBmb3'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).TotalGroundedBmb),... + (md.results.TransientSolution(1).TotalFloatingBmb),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).TotalGroundedBmb),... + (md.results.TransientSolution(2).TotalFloatingBmb),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).TotalGroundedBmb),... + (md.results.TransientSolution(3).TotalFloatingBmb),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test249.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test249.py (revision 26438) +++ /issm/trunk-jpl/test/NightlyRun/test249.py (revision 26438) @@ -0,0 +1,62 @@ +#Test Name: SquareShelfTranMLHO2d + +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 150000) +md = setmask(md, 'all', '') +md = parameterize(md, '../Par/SquareShelf.py') +md.basalforcings.floatingice_melting_rate[:] = 1. +md = setflowequation(md, 'MLHO', 'all') +md.cluster = generic('name', gethostname(), 'np', 3) +md.transient.requested_outputs = ['default', 'FloatingArea', 'GroundedArea', 'TotalFloatingBmb', 'TotalGroundedBmb'] +md = solve(md, 'Transient') + + +# Fields and tolerances to track changes + +field_names = ['Vx1', 'Vy1', 'Vel1', 'Pressure1', + 'Bed1', 'Surface1', 'Thickness1', 'TotalGroundedBmb1', 'TotalFloatingBmb1', + 'Vx2', 'Vy2', 'Vel2', 'Pressure2', + 'Bed2', 'Surface2', 'Thickness2', 'TotalGroundedBmb2', 'TotalFloatingBmb2', + 'Vx3', 'Vy3', 'Vel3', 'Pressure3', + 'Bed3', 'Surface3', 'Thickness3', 'TotalGroundedBmb3', 'TotalFloatingBmb3'] +field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13] +field_values = [md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].TotalGroundedBmb, + md.results.TransientSolution[0].TotalFloatingBmb, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].TotalGroundedBmb, + md.results.TransientSolution[1].TotalFloatingBmb, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].TotalGroundedBmb, + md.results.TransientSolution[2].TotalFloatingBmb] Index: /issm/trunk-jpl/test/NightlyRun/test254.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test254.m (revision 26438) +++ /issm/trunk-jpl/test/NightlyRun/test254.m (revision 26438) @@ -0,0 +1,72 @@ +%Test Name: SquareShelfConstrainedMLHO +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md=setflowequation(md,'MLHO','all'); +md.cluster=generic('name',oshostname(),'np',3); + +%redo the parameter file for this special shelf. +%constant thickness, constrained (vy=0) flow into an icefront, +%from 0 m/yr at the grounding line. + +%tighten +md.stressbalance.restol=10^-4; + +%needed later +ymin=min(md.mesh.y); +ymax=max(md.mesh.y); +xmin=min(md.mesh.x); +xmax=max(md.mesh.x); + +di=md.materials.rho_ice/md.materials.rho_water; + +h=1000; +md.geometry.thickness=h*ones(md.mesh.numberofvertices,1); +md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; +md.geometry.surface=md.geometry.base+md.geometry.thickness; + +%Initial velocity and pressure +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.vz=zeros(md.mesh.numberofvertices,1); +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); + +%Materials +md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1); +md.materials.rheology_B=paterson(md.initialization.temperature); +md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + +%Boundary conditions: +md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1); +md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1); +md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1); + +%constrain flanks to 0 normal velocity +pos=find(md.mesh.x==xmin | md.mesh.x==xmax); +md.stressbalance.spcvx(pos)=0; +md.stressbalance.spcvz(pos)=NaN; + +%constrain grounding line to 0 velocity +pos=find(md.mesh.y==ymin); +md.stressbalance.spcvx(pos)=0; +md.stressbalance.spcvy(pos)=0; + +%icefront +nodeonicefront=zeros(md.mesh.numberofvertices,1); +pos=find(md.mesh.y==ymax); nodeonicefront(pos)=1; +md.mask.ice_levelset=-1+nodeonicefront; + +md.stressbalance.requested_outputs={'default','VySurface','VyShear','VyBase'}; +md=solve(md,'Stressbalance'); + +%create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292. +%ey_c=(md.materials.rho_ice*md.constants.g*(1-di)*md.geometry.thickness./(4*md.materials.rheology_B)).^3; +%vy_c=ey_c.*md.mesh.y*md.constants.yts; + +%Fields and tolerances to track changes +field_names ={'Vy','VySurface','VyShear','VyBase'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13}; +field_values={(md.results.StressbalanceSolution.Vy),... + (md.results.StressbalanceSolution.VySurface),... + (md.results.StressbalanceSolution.VyShear),... + (md.results.StressbalanceSolution.VyBase)}; Index: /issm/trunk-jpl/test/NightlyRun/test254.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test254.py (revision 26438) +++ /issm/trunk-jpl/test/NightlyRun/test254.py (revision 26438) @@ -0,0 +1,81 @@ +#Test Name: SquareShelfConstrainedMLHO +from model import * +from socket import gethostname +import numpy as np +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from paterson import * +from solve import * +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 150000) +md = setmask(md, 'all', '') +md = parameterize(md, '../Par/SquareShelf.py') +md = setflowequation(md, 'MLHO', 'all') +md.cluster = generic('name', gethostname(), 'np', 3) + +# redo the parameter file for this special shelf. +# constant thickness, constrained (vy = 0) flow into an icefront, +# from 0 m / yr at the grounding line. + +# tighten +md.stressbalance.restol = 1e-4 + +# needed later +ymin = min(md.mesh.y) +ymax = max(md.mesh.y) +xmin = min(md.mesh.x) +xmax = max(md.mesh.x) + +di = md.materials.rho_ice / md.materials.rho_water + +h = 1000. +md.geometry.thickness = h * np.ones((md.mesh.numberofvertices)) +md.geometry.base = -md.materials.rho_ice / md.materials.rho_water * md.geometry.thickness +md.geometry.surface = md.geometry.base + md.geometry.thickness + +# Initial velocity and pressure +md.initialization.vx = np.zeros((md.mesh.numberofvertices)) +md.initialization.vy = np.zeros((md.mesh.numberofvertices)) +md.initialization.vz = np.zeros((md.mesh.numberofvertices)) +md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) + +# Materials +md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices)) +md.materials.rheology_B = paterson(md.initialization.temperature) +md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements)) + +# Boundary conditions: +md.stressbalance.spcvx = np.nan * np.ones((md.mesh.numberofvertices)) +md.stressbalance.spcvy = np.nan * np.ones((md.mesh.numberofvertices)) +md.stressbalance.spcvz = np.nan * np.ones((md.mesh.numberofvertices)) + +# constrain flanks to 0 normal velocity +pos = np.where(np.logical_or(md.mesh.x == xmin, md.mesh.x == xmax)) +md.stressbalance.spcvx[pos] = 0 +md.stressbalance.spcvz[pos] = np.nan + +# constrain grounding line to 0 velocity +pos = np.where(md.mesh.y == ymin) +md.stressbalance.spcvx[pos] = 0 +md.stressbalance.spcvy[pos] = 0 + +# icefront +nodeonicefront = np.zeros(md.mesh.numberofvertices) +pos = np.where(md.mesh.y == ymax) +nodeonicefront[pos] = 1 +md.mask.ice_levelset = -1 + nodeonicefront + +md.stressbalance.requested_outputs = ['default', 'VySurface', 'VyShear', 'VyBase'] +md = solve(md, 'Stressbalance') + +# create analytical solution: strain rate is constant = ((rho_ice * g * h) / 4B)^3 (Paterson, 4th Edition, page 292. +# ey_c=(md.materials.rho_ice * md.constants.g * (1 - di) * md.geometry.thickness. / (4 * md.materials.rheology_B)).^3 +# vy_c = ey_c. * md.mesh.y * md.constants.yts + +# Fields and tolerances to track changes +field_names = ['Vy', 'VySurface', 'VyShear', 'VyBase'] +field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13] +field_values = [md.results.StressbalanceSolution.Vy, md.results.StressbalanceSolution.VySurface, md.results.StressbalanceSolution.VyShear, md.results.StressbalanceSolution.VyBase]