Index: /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par (revision 14105) +++ /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par (revision 14106) @@ -1,10 +1,10 @@ %Ok, start defining model parameters here disp(' creating thickness'); -md.geometry.thickness=10*ones(md.mesh.numberofvertices,1); +md.geometry.thickness=10.*ones(md.mesh.numberofvertices,1); md.geometry.bed=zeros(md.mesh.numberofvertices,1); md.geometry.surface=md.geometry.bed+md.geometry.thickness; disp(' creating drag'); -md.friction.coefficient=20*ones(md.mesh.numberofvertices,1); %q=1. %no drag is specified in the analytical solution +md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); %q=1. no drag is specified in the analytical solution md.friction.p=ones(md.mesh.numberofelements,1); md.friction.q=ones(md.mesh.numberofelements,1); @@ -12,23 +12,23 @@ disp(' creating temperatures'); tmin=238.15; %K -st=1.67*10^-2/1000; %k/m; +st=1.67*10^-2/1000.; %k/m radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2); md.initialization.temperature=(tmin+st*radius); md.basalforcings.geothermalflux=4.2*10^-2*ones(md.mesh.numberofvertices,1); -disp(' creating flow law paramter'); -md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +disp(' creating flow law parameter'); +md.materials.rheology_B=6.81*10^7*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution +md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); disp(' creating surface mass balance'); smb_max=0.5; %m/yr -sb=10^-2/1000; %m/yr/m -rel=450*1000; %m +sb=10^-2/1000.; %m/yr/m +rel=450.*1000.; %m md.surfaceforcings.mass_balance=min(smb_max,sb*(rel-radius)); disp(' creating velocities'); constant=0.3; -md.inversion.vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; -md.inversion.vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; +md.inversion.vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1; +md.inversion.vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1; md.inversion.vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2)); md.initialization.vx=zeros(md.mesh.numberofvertices,1); @@ -38,22 +38,22 @@ %Deal with boundary conditions: -disp(' boundary conditions for diagnostic model: '); +disp(' boundary conditions for diagnostic model:'); md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp'); -radius=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); +radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2); pos=find(radius==min(radius)); -md.mesh.x(pos)=0; md.mesh.y(pos)=0; %the closest node to the center is changed to be exactly at the center +md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center -md.diagnostic.spcvx(pos)=0; -md.diagnostic.spcvy(pos)=0; -md.diagnostic.spcvz(pos)=0; +md.diagnostic.spcvx(pos)=0.; +md.diagnostic.spcvy(pos)=0.; +md.diagnostic.spcvz(pos)=0.; %parallel options -md.timestepping.final_time=50000; +md.timestepping.final_time=50000.; %Constants -md.materials.rho_ice=910; +md.materials.rho_ice=910.; md.materials.thermalconductivity=2.1; md.materials.latentheat=3.35*10^5; md.materials.beta=8.66*10^-4/(md.materials.rho_ice*md.constants.g); %conversion from K/m to K/Pa -md.constants.yts=31556926; +md.constants.yts=31556926.; Index: /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py (revision 14106) +++ /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py (revision 14106) @@ -0,0 +1,63 @@ +import numpy +from SetMarineIceSheetBC import * + +#Ok, start defining model parameters here +print " creating thickness" +md.geometry.thickness=10.*numpy.ones((md.mesh.numberofvertices,1)) +md.geometry.bed=numpy.zeros((md.mesh.numberofvertices,1)) +md.geometry.surface=md.geometry.bed+md.geometry.thickness + +print " creating drag" +md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1)) #q=1. no drag is specified in the analytical solution +md.friction.p=numpy.ones((md.mesh.numberofelements,1)) +md.friction.q=numpy.ones((md.mesh.numberofelements,1)) + +print " creating temperatures" +tmin=238.15 #K +st=1.67*10**-2/1000. #k/m +radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2).reshape(-1,1) +md.initialization.temperature=tmin+st*radius +md.basalforcings.geothermalflux=4.2*10**-2*numpy.ones((md.mesh.numberofvertices,1)) + +print " creating flow law parameter" +md.materials.rheology_B=6.81*10**7*numpy.ones((md.mesh.numberofvertices,1)) #to have the same B as the analytical solution +md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) + +print " creating surface mass balance" +smb_max=0.5 #m/yr +sb=10**-2/1000. #m/yr/m +rel=450.*1000. #m +md.surfaceforcings.mass_balance=numpy.minimum(smb_max*numpy.ones_like(radius),sb*(rel-radius)) + +print " creating velocities" +constant=0.3 +md.inversion.vx_obs=constant/2.*md.mesh.x.reshape(-1,1)*(md.geometry.thickness)**-1 +md.inversion.vy_obs=constant/2.*md.mesh.y.reshape(-1,1)*(md.geometry.thickness)**-1 +md.inversion.vel_obs=numpy.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2) +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)) + +#Deal with boundary conditions: +print " boundary conditions for diagnostic model:" +md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp') + +radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2) +pos=numpy.nonzero(radius==numpy.min(radius))[0] +md.mesh.x[pos]=0. +md.mesh.y[pos]=0. #the closest node to the center is changed to be exactly at the center + +md.diagnostic.spcvx[pos]=0. +md.diagnostic.spcvy[pos]=0. +md.diagnostic.spcvz[pos]=0. + +#parallel options +md.timestepping.final_time=50000. + +#Constants +md.materials.rho_ice=910. +md.materials.thermalconductivity=2.1 +md.materials.latentheat=3.35*10**5 +md.materials.beta=8.66*10**-4/(md.materials.rho_ice*md.constants.g) #conversion from K/m to K/Pa +md.constants.yts=31556926. Index: /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par (revision 14105) +++ /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par (revision 14106) @@ -2,15 +2,15 @@ hmin=0.01; hmax=2756.7; -radius=(sqrt((md.mesh.x).^2+(md.mesh.y).^2)); +radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2); radiusmax=max(radius); -md.geometry.thickness=hmin*ones(size(md.mesh.x,1),1)+hmax*(4*((1/2)^(4/3)*ones(size(md.mesh.x,1),1)-((radius)./(2*radiusmax)).^(4/3))).^(3/8); -md.geometry.bed=0*md.geometry.thickness; +md.geometry.thickness=hmin*ones(size(md.mesh.x,1),1)+hmax*(4.*((1./2.)^(4./3.)*ones(size(md.mesh.x,1),1)-((radius)./(2.*radiusmax)).^(4./3.))).^(3./8.); +md.geometry.bed=0.*md.geometry.thickness; md.geometry.surface=md.geometry.bed+md.geometry.thickness; disp(' creating drag'); -md.friction.coefficient=20*ones(md.mesh.numberofvertices,1); %q=1. %no drag is specified in the analytical solution +md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); %q=1. no drag is specified in the analytical solution %Take care of iceshelves: no basal drag pos=find(md.mask.elementonfloatingice); -md.friction.coefficient(md.mesh.elements(pos,:))=0; +md.friction.coefficient(md.mesh.elements(pos,:))=0.; md.friction.p=ones(md.mesh.numberofelements,1); md.friction.q=ones(md.mesh.numberofelements,1); @@ -18,23 +18,23 @@ disp(' creating temperatures'); tmin=238.15; %K -st=1.67*10^-2/1000; %k/m; -md.initialization.temperature=(tmin+st*radius); +st=1.67*10^-2/1000.; %k/m +md.initialization.temperature=tmin+st*radius; md.basalforcings.geothermalflux=4.2*10^-2*ones(md.mesh.numberofvertices,1); -disp(' creating flow law paramter'); -md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +disp(' creating flow law parameter'); +md.materials.rheology_B=6.81*10^7*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution +md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); disp(' creating surface mass balance'); smb_max=0.5; %m/yr -sb=10^-2/1000; %m/yr/m -rel=450*1000; %m +sb=10^-2/1000.; %m/yr/m +rel=450.*1000.; %m md.surfaceforcings.mass_balance=min(smb_max,sb*(rel-radius)); disp(' creating velocities'); constant=0.3; -md.inversion.vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; -md.inversion.vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; -md.inversion.vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2)); +md.inversion.vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1; +md.inversion.vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1; +md.inversion.vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2); md.initialization.vx=zeros(md.mesh.numberofvertices,1); md.initialization.vy=zeros(md.mesh.numberofvertices,1); @@ -43,12 +43,12 @@ %Deal with boundary conditions: -disp(' boundary conditions for diagnostic model: '); +disp(' boundary conditions for diagnostic model:'); md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp'); -radius=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); +radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2); pos=find(radius==min(radius)); -md.mesh.x(pos)=0; md.mesh.y(pos)=0; %the closest node to the center is changed to be exactly at the center +md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center -md.diagnostic.spcvx(pos)=0; -md.diagnostic.spcvy(pos)=0; -md.diagnostic.spcvz(pos)=0; +md.diagnostic.spcvx(pos)=0.; +md.diagnostic.spcvy(pos)=0.; +md.diagnostic.spcvz(pos)=0.; Index: /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py (revision 14106) +++ /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py (revision 14106) @@ -0,0 +1,58 @@ +import numpy +from SetMarineIceSheetBC import * + +print " creating thickness" +hmin=0.01 +hmax=2756.7 +radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2).reshape(-1,1) +radiusmax=numpy.max(radius) +md.geometry.thickness=hmin*numpy.ones((numpy.size(md.mesh.x),1))+hmax*(4.*((1./2.)**(4./3.)*numpy.ones((numpy.size(md.mesh.x),1))-((radius)/(2.*radiusmax))**(4./3.)))**(3./8.) +md.geometry.bed=0.*md.geometry.thickness +md.geometry.surface=md.geometry.bed+md.geometry.thickness + +print " creating drag" +md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1)) #q=1. no drag is specified in the analytical solution +#Take care of iceshelves: no basal drag +pos=numpy.nonzero(md.mask.elementonfloatingice)[0] +md.friction.coefficient[md.mesh.elements[pos,:]-1]=0. +md.friction.p=numpy.ones((md.mesh.numberofelements,1)) +md.friction.q=numpy.ones((md.mesh.numberofelements,1)) + +print " creating temperatures" +tmin=238.15 #K +st=1.67*10**-2/1000. #k/m +md.initialization.temperature=tmin+st*radius +md.basalforcings.geothermalflux=4.2*10**-2*numpy.ones((md.mesh.numberofvertices,1)) + +print " creating flow law parameter" +md.materials.rheology_B=6.81*10**7*numpy.ones((md.mesh.numberofvertices,1)) #to have the same B as the analytical solution +md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) + +print " creating surface mass balance" +smb_max=0.5 #m/yr +sb=10**-2/1000. #m/yr/m +rel=450.*1000. #m +md.surfaceforcings.mass_balance=numpy.minimum(smb_max*numpy.ones_like(radius),sb*(rel-radius)) + +print " creating velocities" +constant=0.3 +md.inversion.vx_obs=constant/2.*md.mesh.x.reshape(-1,1)*(md.geometry.thickness)**-1 +md.inversion.vy_obs=constant/2.*md.mesh.y.reshape(-1,1)*(md.geometry.thickness)**-1 +md.inversion.vel_obs=numpy.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2) +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)) + +#Deal with boundary conditions: +print " boundary conditions for diagnostic model:" +md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp') + +radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2) +pos=numpy.nonzero(radius==numpy.min(radius))[0] +md.mesh.x[pos]=0. +md.mesh.y[pos]=0. #the closest node to the center is changed to be exactly at the center + +md.diagnostic.spcvx[pos]=0. +md.diagnostic.spcvy[pos]=0. +md.diagnostic.spcvz[pos]=0. Index: /issm/trunk-jpl/test/Par/SquareEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareEISMINT.par (revision 14105) +++ /issm/trunk-jpl/test/Par/SquareEISMINT.par (revision 14106) @@ -4,18 +4,18 @@ ymin=min(md.mesh.y); ymax=max(md.mesh.y); -md.geometry.thickness=500*ones(md.mesh.numberofvertices,1); +md.geometry.thickness=500.*ones(md.mesh.numberofvertices,1); md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; md.geometry.surface=md.geometry.bed+md.geometry.thickness; disp(' creating drag'); -md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1. +md.friction.coefficient=200.*ones(md.mesh.numberofvertices,1); %q=1. %Take care of iceshelves: no basal drag pos=find(md.mask.elementonfloatingice); -md.friction.coefficient(md.mesh.elements(pos,:))=0; +md.friction.coefficient(md.mesh.elements(pos,:))=0.; md.friction.p=ones(md.mesh.numberofelements,1); md.friction.q=ones(md.mesh.numberofelements,1); disp(' creating initial values'); -md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1); +md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1); md.initialization.vx=zeros(md.mesh.numberofvertices,1); md.initialization.vy=zeros(md.mesh.numberofvertices,1); @@ -24,23 +24,23 @@ md.initialization.pressure=zeros(md.mesh.numberofvertices,1); -disp(' creating flow law paramter'); +disp(' creating flow law parameter'); md.materials.rheology_B=1.7687*10^8*ones(md.mesh.numberofvertices,1); -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); disp(' creating surface mass balance'); md.surfaceforcings.mass_balance=0.2*ones(md.mesh.numberofvertices,1); %0m/a -md.basalforcings.melting_rate=0*ones(md.mesh.numberofvertices,1); %0m/a +md.basalforcings.melting_rate=0.*ones(md.mesh.numberofvertices,1); %0m/a -disp(' boundary conditions '); +disp(' boundary conditions'); md=SetMarineIceSheetBC(md,'../Exp/SquareFrontEISMINT.exp'); %Evolution of the ice shelf -pos=find(md.mesh.y==200000); %nodes on the upper boundary condition +pos=find(md.mesh.y==200000.); %nodes on the upper boundary condition md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1); -md.balancethickness.spcthickness(pos)=500; +md.balancethickness.spcthickness(pos)=500.; md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1); -md.prognostic.spcthickness(pos)=500; +md.prognostic.spcthickness(pos)=500.; md.prognostic.stabilization=0; %Better result with no artificial diffusivity -md.thermal.stabilization=0; -md.timestepping.final_time=500; +md.thermal.stabilization=0; +md.timestepping.final_time=500.; md.timestepping.time_step=1; Index: /issm/trunk-jpl/test/Par/SquareEISMINT.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareEISMINT.py (revision 14106) +++ /issm/trunk-jpl/test/Par/SquareEISMINT.py (revision 14106) @@ -0,0 +1,49 @@ +import numpy +from SetMarineIceSheetBC import * + +#Ok, start defining model parameters here + +print " creating thickness" +ymin=numpy.min(md.mesh.y) +ymax=numpy.max(md.mesh.y) +md.geometry.thickness=500.*numpy.ones((md.mesh.numberofvertices,1)) +md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness +md.geometry.surface=md.geometry.bed+md.geometry.thickness + +print " creating drag" +md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1)) #q=1. +#Take care of iceshelves: no basal drag +pos=numpy.nonzero(md.mask.elementonfloatingice)[0] +md.friction.coefficient[md.mesh.elements[pos,:]-1]=0. +md.friction.p=numpy.ones((md.mesh.numberofelements,1)) +md.friction.q=numpy.ones((md.mesh.numberofelements,1)) + +print " creating initial values" +md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1)) +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.vel=numpy.zeros((md.mesh.numberofvertices,1)) +md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1)) + +print " creating flow law parameter" +md.materials.rheology_B=1.7687*10**8*numpy.ones((md.mesh.numberofvertices,1)) +md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) + +print " creating surface mass balance" +md.surfaceforcings.mass_balance=0.2*numpy.ones((md.mesh.numberofvertices,1)) #0m/a +md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices,1)) #0m/a + +print " boundary conditions" +md=SetMarineIceSheetBC(md,'../Exp/SquareFrontEISMINT.exp') + +#Evolution of the ice shelf +pos=numpy.nonzero(md.mesh.y==200000.) #nodes on the upper boundary condition +md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) +md.balancethickness.spcthickness[pos]=500. +md.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) +md.prognostic.spcthickness[pos]=500. +md.prognostic.stabilization=0 #Better result with no artificial diffusivity +md.thermal.stabilization=0 +md.timestepping.final_time=500. +md.timestepping.time_step=1