Index: /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt =================================================================== --- /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt (revision 13808) +++ /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt (revision 13809) @@ -11,5 +11,4 @@ test418 needs Dakota test420 needs Dakota -test423 needs roundmesh (and not straightforward) test511 needs extract test613 needs extract Index: /issm/trunk-jpl/test/NightlyRun/test423.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test423.m (revision 13808) +++ /issm/trunk-jpl/test/NightlyRun/test423.m (revision 13809) @@ -1,4 +1,4 @@ -radius=1e6; -shelfextent=2e5; +radius=1.e6; +shelfextent=2.e5; md=roundmesh(model(),radius,50000.); @@ -6,7 +6,7 @@ rad=sqrt(md.mesh.x.^2+md.mesh.y.^2); pos=find(rad==min(rad)); -md.mesh.x(pos)=0; md.mesh.y(pos)=0; %the closest node to the center is changed to be exactly at the center -xelem=md.mesh.x(md.mesh.elements)*[1;1;1]/3; -yelem=md.mesh.y(md.mesh.elements)*[1;1;1]/3; +md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center +xelem=md.mesh.x(md.mesh.elements)*[1;1;1]/3.; +yelem=md.mesh.y(md.mesh.elements)*[1;1;1]/3.; rad=sqrt(xelem.^2+yelem.^2); flags=zeros(md.mesh.numberofelements,1); Index: /issm/trunk-jpl/test/NightlyRun/test423.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test423.py (revision 13809) +++ /issm/trunk-jpl/test/NightlyRun/test423.py (revision 13809) @@ -0,0 +1,48 @@ +import numpy +from model import * +from roundmesh import * +from setmask import * +from parameterize import * +from setflowequation import * +from EnumDefinitions import * +from solve import * +from MatlabFuncs import * + +radius=1.e6 +shelfextent=2.e5 + +md=roundmesh(model(),radius,50000.) +#fix center node to 0,0 +rad=numpy.sqrt(md.mesh.x**2+md.mesh.y**2) +pos=numpy.argmin(rad) +md.mesh.x[pos]=0. +md.mesh.y[pos]=0. #the closest node to the center is changed to be exactly at the center +xelem=numpy.mean(md.mesh.x[md.mesh.elements.astype(int)-1],axis=1) +yelem=numpy.mean(md.mesh.y[md.mesh.elements.astype(int)-1],axis=1) +rad=numpy.sqrt(xelem**2+yelem**2) +flags=numpy.zeros(md.mesh.numberofelements) +pos=numpy.nonzero(rad>=(radius-shelfextent)) +flags[pos]=1 +md=setmask(md,flags,'') +md=parameterize(md,'../Par/RoundSheetShelf.py') +md=setflowequation(md,'macayeal','all') +md.cluster=generic('name',oshostname(),'np',3) + +md.transient.isthermal=0 +md.transient.isprognostic=0 +md.transient.isdiagnostic=0 +md.transient.isgroundingline=1 + +#test different grounding line dynamics. +md.groundingline.migration='AgressiveMigration' +md=solve(md,TransientSolutionEnum()) +element_on_iceshelf_agressive=md.results['TransientSolution'][1]['MaskElementonfloatingice'] + +md.groundingline.migration='SoftMigration' +md=solve(md,TransientSolutionEnum()) +element_on_iceshelf_soft=md.results['TransientSolution'][1]['MaskElementonfloatingice'] + +#Fields and tolerances to track changes +field_names =['ElementOnIceShelfAgressive','ElementOnIceShelfSoft'] +field_tolerances=[1e-13,1e-13] +field_values=[element_on_iceshelf_agressive,element_on_iceshelf_soft] Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 13808) +++ /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 13809) @@ -2,9 +2,9 @@ di=md.materials.rho_ice/md.materials.rho_water; -rad=1e6; -shelfextent=2e5; +rad=1.e6; +shelfextent=2.e5; %Geometry -hmin=300; -hmax=1000; +hmin=300.; +hmax=1000.; radius=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); ymin=min(radius); @@ -14,21 +14,20 @@ pos=find(md.mask.vertexongroundedice); -md.geometry.bed(pos)=md.geometry.bed(pos)-300*(radius(pos)-(rad-shelfextent))/(rad-shelfextent); +md.geometry.bed(pos)=md.geometry.bed(pos)-300.*(radius(pos)-(rad-shelfextent))/(rad-shelfextent); md.geometry.surface=md.geometry.bed+md.geometry.thickness; pos=find(radius<200000); -md.geometry.thickness(pos)=100; -md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; +md.geometry.thickness(pos)=100.; +md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20.; md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); - -pos=find(md.mesh.x<.2*1e6 & md.mesh.x>-.2*1e6 & md.mesh.y>0); -md.geometry.thickness(pos)=100; -md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; +pos=find(md.mesh.x<0.2*1.e6 & md.mesh.x>-0.2*1.e6 & md.mesh.y>0.); +md.geometry.thickness(pos)=100.; +md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20.; md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); -pos=find(md.mesh.x<.1*1e6 & md.mesh.x>-.1*1e6 & md.mesh.y<-.5*1e6 & md.mesh.y>-.6*1e6); -md.geometry.thickness(pos)=100; -md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; +pos=find(md.mesh.x<0.1*1.e6 & md.mesh.x>-0.1*1.e6 & md.mesh.y<-0.5*1.e6 & md.mesh.y>-0.6*1.e6); +md.geometry.thickness(pos)=100.; +md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20.; md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); @@ -36,18 +35,19 @@ di=md.materials.rho_ice/md.materials.rho_water; rad=sqrt(md.mesh.x.^2+md.mesh.y.^2); -pos=find(rad<200000); -md.geometry.thickness(pos)=100; -md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; +pos=find(rad<200000.); +md.geometry.thickness(pos)=100.; +md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20.; md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); -pos=find(md.mesh.x<.2*1e6 & md.mesh.x>-.2*1e6 & md.mesh.y>0); -md.geometry.thickness(pos)=100; -md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; +pos=find(md.mesh.x<0.2*1.e6 & md.mesh.x>-0.2*1.e6 & md.mesh.y>0.); +md.geometry.thickness(pos)=100.; +md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20.; md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); -pos=find(md.mesh.x<.1*1e6 & md.mesh.x>-.1*1e6 & md.mesh.y<-.5*1e6 & md.mesh.y>-.6*1e6); -md.geometry.thickness(pos)=100; -md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; +pos=find(md.mesh.x<0.1*1.e6 & md.mesh.x>-0.1*1.e6 & md.mesh.y<-0.5*1.e6 & md.mesh.y>-0.6*1.e6); +md.geometry.thickness(pos)=100.; +md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20.; md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); + %Initial velocity md.initialization.vx=zeros(md.mesh.numberofvertices,1); @@ -57,17 +57,17 @@ %Materials -md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1); +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); +md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); %Surface mass balance and basal melting -md.surfaceforcings.mass_balance=-10*ones(md.mesh.numberofvertices,1); +md.surfaceforcings.mass_balance=-10.*ones(md.mesh.numberofvertices,1); md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1); -pos=find(md.mask.vertexonfloatingice);md.basalforcings.melting_rate(pos)=10; +pos=find(md.mask.vertexonfloatingice);md.basalforcings.melting_rate(pos)=10.; md.basalforcings.geothermalflux=ones(md.mesh.numberofvertices,1); %Friction pos=find(md.mask.elementonfloatingice); -md.friction.coefficient=20*ones(md.mesh.numberofvertices,1); +md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); md.friction.coefficient(md.mesh.elements(pos,:))=0; md.friction.p=ones(md.mesh.numberofelements,1); @@ -84,11 +84,11 @@ md.steadystate.reltol=0.05; md.diagnostic.abstol=NaN; -md.timestepping.time_step=5; -md.timestepping.final_time=5; +md.timestepping.time_step=5.; +md.timestepping.final_time=5.; %bathymetry and grounding line migration: md.groundingline.migration='AgressiveMigration'; md.geometry.bathymetry=md.geometry.bed; -pos=find(md.mask.vertexonfloatingice); md.geometry.bathymetry(pos)=md.geometry.bed(pos)-900; +pos=find(md.mask.vertexonfloatingice); md.geometry.bathymetry(pos)=md.geometry.bed(pos)-900.; %Deal with boundary conditions: @@ -107,5 +107,5 @@ md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1); md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6); -md.thermal.spctemperature=737*ones(md.mesh.numberofvertices,1); +md.thermal.spctemperature=737.*ones(md.mesh.numberofvertices,1); %Change name so that no test have the same name Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetShelf.py (revision 13809) +++ /issm/trunk-jpl/test/Par/RoundSheetShelf.py (revision 13809) @@ -0,0 +1,122 @@ +import os.path +import numpy +import copy +import inspect +from paterson import * +from verbose import * + +#Start defining model parameters here + +di=md.materials.rho_ice/md.materials.rho_water +rad=1.e6 +shelfextent=2.e5 +#Geometry +hmin=300. +hmax=1000. +radius=numpy.sqrt(md.mesh.x*md.mesh.x+md.mesh.y*md.mesh.y) +ymin=numpy.min(radius) +ymax=numpy.max(radius) +md.geometry.thickness=hmax+(hmin-hmax)*(radius-ymin)/(ymax-ymin) +md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness + +pos=numpy.nonzero(md.mask.vertexongroundedice) +md.geometry.bed[pos]=md.geometry.bed[pos]-300.*(radius[pos]-(rad-shelfextent))/(rad-shelfextent) +md.geometry.surface=md.geometry.bed+md.geometry.thickness + +pos=numpy.nonzero(radius<200000.) +md.geometry.thickness[pos]=100. +md.geometry.bed[pos]=-di*md.geometry.thickness[pos]-20. +md.geometry.surface[pos]=md.geometry.bed[pos]+md.geometry.thickness[pos] + +pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.2*1.e6,md.mesh.x>-0.2*1.e6),md.mesh.y>0.)) +md.geometry.thickness[pos]=100. +md.geometry.bed[pos]=-di*md.geometry.thickness[pos]-20. +md.geometry.surface[pos]=md.geometry.bed[pos]+md.geometry.thickness[pos] + +pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.1*1.e6,md.mesh.x>-0.1*1.e6),numpy.logical_and(md.mesh.y<-0.5*1.e6,md.mesh.y>-0.6*1.e6))) +md.geometry.thickness[pos]=100. +md.geometry.bed[pos]=-di*md.geometry.thickness[pos]-20. +md.geometry.surface[pos]=md.geometry.bed[pos]+md.geometry.thickness[pos] + +#plug holes into the ice sheet, to test for grounding line migration. +di=md.materials.rho_ice/md.materials.rho_water +rad=numpy.sqrt(md.mesh.x**2+md.mesh.y**2) +pos=numpy.nonzero(rad<200000.) +md.geometry.thickness[pos]=100. +md.geometry.bed[pos]=-di*md.geometry.thickness[pos]-20. +md.geometry.surface[pos]=md.geometry.bed[pos]+md.geometry.thickness[pos] + +pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.2*1.e6,md.mesh.x>-0.2*1.e6),md.mesh.y>0.)) +md.geometry.thickness[pos]=100. +md.geometry.bed[pos]=-di*md.geometry.thickness[pos]-20. +md.geometry.surface[pos]=md.geometry.bed[pos]+md.geometry.thickness[pos] + +pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x<0.1*1.e6,md.mesh.x>-0.1*1.e6),numpy.logical_and(md.mesh.y<-0.5*1.e6,md.mesh.y>-0.6*1.e6))) +md.geometry.thickness[pos]=100. +md.geometry.bed[pos]=-di*md.geometry.thickness[pos]-20. +md.geometry.surface[pos]=md.geometry.bed[pos]+md.geometry.thickness[pos] + +#Initial velocity +md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1)) +md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1)) +md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1)) +md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1)) + +#Materials +md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1)) +md.materials.rheology_B=paterson(md.initialization.temperature) +md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) + +#Surface mass balance and basal melting +md.surfaceforcings.mass_balance=-10.*numpy.ones((md.mesh.numberofvertices,1)) +md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1)) +pos=numpy.nonzero(md.mask.vertexonfloatingice) +md.basalforcings.melting_rate[pos]=10. +md.basalforcings.geothermalflux=numpy.ones((md.mesh.numberofvertices,1)) + +#Friction +pos=numpy.nonzero(md.mask.elementonfloatingice) +md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1)) +md.friction.coefficient[md.mesh.elements[pos,:].astype(int)-1]=0 +md.friction.p=numpy.ones((md.mesh.numberofelements,1)) +md.friction.q=numpy.ones((md.mesh.numberofelements,1)) + +#Numerical parameters +md.diagnostic.viscosity_overshoot=0.0 +md.prognostic.stabilization=1 +md.thermal.stabilization=1 +md.verbose=verbose(0) +md.settings.waitonlock=30 +md.diagnostic.restol=0.05 +md.diagnostic.reltol=0.05 +md.steadystate.reltol=0.05 +md.diagnostic.abstol=float('nan') +md.timestepping.time_step=5. +md.timestepping.final_time=5. + +#bathymetry and grounding line migration: +md.groundingline.migration='AgressiveMigration' +md.geometry.bathymetry=copy.deepcopy(md.geometry.bed) +pos=numpy.nonzero(md.mask.vertexonfloatingice) +md.geometry.bathymetry[pos]=md.geometry.bed[pos]-900. + +#Deal with boundary conditions: +md.diagnostic.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) +md.diagnostic.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) +md.diagnostic.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + +pos=numpy.nonzero(numpy.logical_and(md.mesh.x==0,md.mesh.y==0)) +md.diagnostic.spcvx[pos]=0 +md.diagnostic.spcvy[pos]=0 + +pressureload=copy.deepcopy(md.mesh.segments) +pressureload=numpy.hstack((pressureload,(1*md.mask.elementonfloatingice[pressureload[:,-1].astype(int)-1] + 0*md.mask.elementongroundedice[pressureload[:,-1].astype(int)-1]).reshape(-1,1))) +md.diagnostic.icefront=pressureload +md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) +md.prognostic.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) +md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6)) +md.thermal.spctemperature=737.*numpy.ones((md.mesh.numberofvertices,1)) + +#Change name so that no test have the same name +if len(inspect.stack()) > 2: + md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]