Index: /issm/trunk-jpl/src/m/materials/nye.py =================================================================== --- /issm/trunk-jpl/src/m/materials/nye.py (revision 24095) +++ /issm/trunk-jpl/src/m/materials/nye.py (revision 24096) @@ -29,11 +29,11 @@ if ice_type == 1: # CO2 ice A_const = 10**(10.8) # s^-1 MPa - Q = 63000 # J mol^-1 - n = 7 # Glen's exponent + Q = 63000. # J mol^-1 + n = 7. # Glen's exponent if ice_type == 2: # H2O ice A_const = 9*10**4 # s^-1 MPa - Q = 60000 # J mol^-1 - n = 3 # Glen's exponent + Q = 60000. # J mol^-1 + n = 3. # Glen's exponent A = A_const*np.exp(-1*Q/(T*Rg)) # s^-1 MPa Index: /issm/trunk-jpl/test/NightlyRun/test346.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test346.m (revision 24095) +++ /issm/trunk-jpl/test/NightlyRun/test346.m (revision 24096) @@ -1,3 +1,3 @@ -%Test Name: SquareSheetConstrainedTherTranNyeCO2 +%Test Name: SquareSheetConstrainedTherTranNyeH2O md=triangle(model(),'../Exp/Square.exp',180000.); md=setmask(md,'',''); @@ -16,8 +16,11 @@ %Fields and tolerances to track changes -field_names ={'Temperature','BasalforcingsGroundediceMeltingRate'}; -field_tolerances={1e-13,1e-13}; +field_names ={'Temperature1','BasalforcingsGroundediceMeltingRate1',... + 'Temperature3','BasalforcingsGroundediceMeltingRate3'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13}; field_values={... (md.results.TransientSolution(1).Temperature),... (md.results.TransientSolution(1).BasalforcingsGroundediceMeltingRate),... + (md.results.TransientSolution(3).Temperature),... + (md.results.TransientSolution(3).BasalforcingsGroundediceMeltingRate),... }; Index: /issm/trunk-jpl/test/NightlyRun/test346.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test346.py (revision 24096) +++ /issm/trunk-jpl/test/NightlyRun/test346.py (revision 24096) @@ -0,0 +1,35 @@ +#Test Name: SquareSheetConstrainedTherTranNyeCO2 +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(), '../Exp/Square.exp', 180000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareSheetConstrainedCO2.py') + +md = md.extrude(3, 1.) +md = setflowequation(md, 'SSA', 'all') + +#Transient options +md.cluster = generic('name', gethostname(), 'np', 3) +md.materials.rheology_law='NyeCO2' +md.transient.isstressbalance=0 +md.transient.ismasstransport=0 +md.transient.issmb=1 +md.transient.isthermal=1 +md.transient.isgroundingline=0 +md = solve(md, 'Transient') + +#Fields and tolerances to track changes +field_names = ['Temperature1','BasalforcingsGroundediceMeltingRate1', + 'Temperature3','BasalforcingsGroundediceMeltingRate3'] +field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13] +field_values = [md.results.TransientSolution[0].Temperature, + md.results.TransientSolution[0].BasalforcingsGroundediceMeltingRate, + md.results.TransientSolution[2].Temperature, + md.results.TransientSolution[2].BasalforcingsGroundediceMeltingRate] Index: /issm/trunk-jpl/test/NightlyRun/test347.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test347.m (revision 24095) +++ /issm/trunk-jpl/test/NightlyRun/test347.m (revision 24096) @@ -3,5 +3,4 @@ md=setmask(md,'',''); md=parameterize(md,'../Par/SquareSheetConstrained.par'); -md.materials.rheology_B=nye(md.initialization.temperature,2); md=extrude(md,3,1.); md=setflowequation(md,'SSA','all'); @@ -9,4 +8,6 @@ md.verbose=verbose('convergence',true,'solution',true); md.materials.rheology_law = 'NyeH2O'; +md.materials.rheology_B=nye(md.initialization.temperature,2) + md.transient.isstressbalance=0; md.transient.ismasstransport=0; @@ -17,8 +18,11 @@ %Fields and tolerances to track changes -field_names ={'Temperature','BasalforcingsGroundediceMeltingRate'}; -field_tolerances={1e-13,1e-13}; +field_names ={'Temperature1','BasalforcingsGroundediceMeltingRate1',... + 'Temperature3','BasalforcingsGroundediceMeltingRate3'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13}; field_values={... (md.results.TransientSolution(1).Temperature),... (md.results.TransientSolution(1).BasalforcingsGroundediceMeltingRate),... + (md.results.TransientSolution(3).Temperature),... + (md.results.TransientSolution(3).BasalforcingsGroundediceMeltingRate),... }; Index: /issm/trunk-jpl/test/NightlyRun/test347.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test347.py (revision 24096) +++ /issm/trunk-jpl/test/NightlyRun/test347.py (revision 24096) @@ -0,0 +1,38 @@ +#Test Name: SquareSheetConstrainedTherTranNyeH2O +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from nye import * + +md = triangle(model(), '../Exp/Square.exp', 180000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareSheetConstrained.py') + +md = md.extrude(3, 1.) +md = setflowequation(md, 'SSA', 'all') + +#Transient options +md.cluster = generic('name', gethostname(), 'np', 3) +md.materials.rheology_law='NyeH2O' +md.materials.rheology_B=nye(md.initialization.temperature,2) + +md.transient.isstressbalance=0 +md.transient.ismasstransport=0 +md.transient.issmb=1 +md.transient.isthermal=1 +md.transient.isgroundingline=0 +md = solve(md, 'Transient') + +#Fields and tolerances to track changes +field_names = ['Temperature1','BasalforcingsGroundediceMeltingRate1', + 'Temperature3','BasalforcingsGroundediceMeltingRate3'] +field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13] +field_values = [md.results.TransientSolution[0].Temperature, + md.results.TransientSolution[0].BasalforcingsGroundediceMeltingRate, + md.results.TransientSolution[2].Temperature, + md.results.TransientSolution[2].BasalforcingsGroundediceMeltingRate] Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 24095) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 24096) @@ -24,7 +24,7 @@ x = np.array(archread('../Data/SquareSheetConstrained.arch','x')) y = np.array(archread('../Data/SquareSheetConstrained.arch','y')) -vx = np.array(archread('../Data/SquareSheetConstrained.arch','vx')); -vy = np.array(archread('../Data/SquareSheetConstrained.arch','vy')); -index = archread('../Data/SquareSheetConstrained.arch','index').astype(int); +vx = np.array(archread('../Data/SquareSheetConstrained.arch','vx')) +vy = np.array(archread('../Data/SquareSheetConstrained.arch','vy')) +index = archread('../Data/SquareSheetConstrained.arch','index').astype(int) md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)[0] Index: /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.par (revision 24095) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.par (revision 24096) @@ -20,7 +20,6 @@ xmax=max(md.mesh.x); md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin); -md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; +md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20.; md.geometry.surface=md.geometry.base+md.geometry.thickness; -md.geometry.bed=md.geometry.base-10; %Initial velocity @@ -69,4 +68,5 @@ md.timestepping.time_step=1; md.timestepping.final_time=3; +md.groundingline.migration='None'; %Deal with boundary conditions: Index: /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.py (revision 24096) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrainedCO2.py (revision 24096) @@ -0,0 +1,85 @@ +import os.path +from arch import archread +import numpy as np +import inspect +from verbose import verbose +from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d +from nye import nye +from SetIceSheetBC import SetIceSheetBC + +#Start defining model parameters here +CO2_temp = 150. +CO2_n = 7. +CO2_meltingPoint = 195. +CO2_latentHeat = 199000. +CO2_rhoIce = 1562. +CO2_heatCapacity = 700. +CO2_thermalCond = 0.5 +CO2_dynViscosity = 13.72*10**-6 +CO2_rhoLiquidZeroDeg = 929. +md.materials.rho_ice = CO2_rhoIce +md.materials.rho_freshwater = CO2_rhoLiquidZeroDeg +md.materials.thermalconductivity = CO2_thermalCond +md.materials.heatcapacity = CO2_heatCapacity +md.materials.meltingpoint = CO2_meltingPoint +md.materials.latentheat = CO2_latentHeat +md.materials.mu_water = CO2_dynViscosity + +#Geometry +hmin=300. +hmax=1000. +ymin=np.min(md.mesh.y) +ymax=np.max(md.mesh.y) +xmin=min(md.mesh.x) +xmax=max(md.mesh.x) +md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin) +md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20. +md.geometry.surface=md.geometry.base+md.geometry.thickness + +#Initial velocity +x = np.array(archread('../Data/SquareSheetConstrained.arch','x')) +y = np.array(archread('../Data/SquareSheetConstrained.arch','y')) +vx = np.array(archread('../Data/SquareSheetConstrained.arch','vx')) +vy = np.array(archread('../Data/SquareSheetConstrained.arch','vy')) +index = archread('../Data/SquareSheetConstrained.arch','index').astype(int) + +[md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y) +[md.initialization.vy]=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y) +md.initialization.vz=np.zeros((md.mesh.numberofvertices)) +md.initialization.pressure=np.zeros((md.mesh.numberofvertices)) + +#Materials +md.initialization.temperature=CO2_temp*np.ones((md.mesh.numberofvertices)) +md.materials.rheology_B=nye(md.initialization.temperature,1) +md.materials.rheology_n=CO2_n*np.ones((md.mesh.numberofelements)) + +#Surface mass balance and basal melting +md.smb.mass_balance=10.*np.ones((md.mesh.numberofvertices)) +md.basalforcings.groundedice_melting_rate=5.*np.ones((md.mesh.numberofvertices)) +md.basalforcings.floatingice_melting_rate=5.*np.ones((md.mesh.numberofvertices)) + +#Friction +md.friction.coefficient=20.*np.ones((md.mesh.numberofvertices)) +md.friction.coefficient[np.nonzero(md.mask.groundedice_levelset<0.)[0]]=0. +md.friction.p=np.ones((md.mesh.numberofelements)) +md.friction.q=np.ones((md.mesh.numberofelements)) + +#Numerical parameters +md.masstransport.stabilization=1 +md.thermal.stabilization=1 +md.verbose = verbose(0) +md.settings.waitonlock=30 +md.stressbalance.restol=0.05 +md.stressbalance.reltol=0.05 +md.steadystate.reltol=0.05 +md.stressbalance.abstol=float('nan') +md.timestepping.time_step=1. +md.timestepping.final_time=3. +md.groundingline.migration='None' + +#Deal with boundary conditions: +md = SetIceSheetBC(md) + +#Change name so that no tests have the same name +if len(inspect.stack()) > 2: + md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]