Index: /issm/trunk-jpl/test/NightlyRun/test2001.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 21533) @@ -1,3 +1,3 @@ -%Test Name: SquareSheetConstrainedGia2d +%Test Name: SquareSheetConstrainedGiaIvins2d %GIA test, inspired on test101 md=triangle(model(),'../Exp/Square.exp',100000.); @@ -6,5 +6,5 @@ %% indicate what you want to compute -md.gia.cross_section_shape=1; % for square-edged x-section +md.giaivins.cross_section_shape=1; % for square-edged x-section %% define loading history @@ -20,11 +20,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2001.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.py (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2001.py (revision 21533) @@ -1,3 +1,3 @@ -#Test Name: SquareSheetConstrainedGia2d +#Test Name: SquareSheetConstrainedGiaIvins2d import numpy as np import copy @@ -18,5 +18,5 @@ #Indicate what you want to compute -md.gia.cross_section_shape=1 # for square-edged x-section +md.giaivins.cross_section_shape=1 # for square-edged x-section #Define loading history (see test2001.m for the description) @@ -31,10 +31,10 @@ #Solve for GIA deflection md.cluster=generic('name',gethostname(),'np',3) -md=solve(md,'Gia') +md=solve(md,'GiaIvins') #Fields and tolerances to track changes -field_names =['GiaW','GiadWdt'] +field_names =['GiaIvinsW','GiaIvinsdWdt'] field_tolerances=[1e-13,1e-13] -field_values =[md.results.GiaSolution.GiaW, - md.results.GiaSolution.GiadWdt] +field_values =[md.results.GiaIvinsSolution.GiaIvinsW, + md.results.GiaIvinsSolution.GiaIvinsdWdt] Index: /issm/trunk-jpl/test/NightlyRun/test2051.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2051.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2051.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dA1 +%Test Name: GiaIvinsBenchmarksAB2dA1 % Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksAB.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksAB.par'); %% indicate what you want to compute -md.gia.cross_section_shape=1; % for square-edged x-section +md.giaivins.cross_section_shape=1; % for square-edged x-section %% define loading history @@ -22,11 +22,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2052.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2052.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2052.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dA2 +%Test Name: GiaIvinsBenchmarksAB2dA2 % Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksAB.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksAB.par'); %% indicate what you want to compute -md.gia.cross_section_shape=1; % for square-edged x-section +md.giaivins.cross_section_shape=1; % for square-edged x-section %% define loading history @@ -22,11 +22,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2053.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2053.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2053.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dA3 +%Test Name: GiaIvinsBenchmarksAB2dA3 % Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksAB.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksAB.par'); %% indicate what you want to compute -md.gia.cross_section_shape=1; % for square-edged x-section +md.giaivins.cross_section_shape=1; % for square-edged x-section %% define loading history @@ -22,11 +22,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2071.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2071.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2071.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dC1 +%Test Name: GiaIvinsBenchmarksAB2dC1 % Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksCD.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par'); %% indicate what you want to compute -md.gia.cross_section_shape=1; % for square-edged x-section +md.giaivins.cross_section_shape=1; % for square-edged x-section %% define loading history @@ -21,11 +21,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2072.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2072.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2072.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dC2 +%Test Name: GiaIvinsBenchmarksAB2dC2 % Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksCD.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par'); %% indicate what you want to compute -md.gia.cross_section_shape=1; % for square-edged x-section +md.giaivins.cross_section_shape=1; % for square-edged x-section %% define loading history @@ -21,11 +21,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2073.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2073.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2073.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dC3 +%Test Name: GiaIvinsBenchmarksAB2dC3 % Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksCD.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par'); %% indicate what you want to compute -md.gia.cross_section_shape=1; % for square-edged x-section +md.giaivins.cross_section_shape=1; % for square-edged x-section %% define loading history @@ -21,11 +21,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2081.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2081.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2081.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dD1 +%Test Name: GiaIvinsBenchmarksAB2dD1 % Benchmark experiments (Figure A2d Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksCD.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par'); %% indicate what you want to compute -md.gia.cross_section_shape=2; % for elliptical edge +md.giaivins.cross_section_shape=2; % for elliptical edge %% define loading history @@ -21,11 +21,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2082.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2082.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2082.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dD2 +%Test Name: GiaIvinsBenchmarksAB2dD2 % Benchmark experiments (Figure A2d Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksCD.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par'); %% indicate what you want to compute -md.gia.cross_section_shape=2; % for elliptical edge +md.giaivins.cross_section_shape=2; % for elliptical edge %% define loading history @@ -21,11 +21,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: /issm/trunk-jpl/test/NightlyRun/test2083.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2083.m (revision 21532) +++ /issm/trunk-jpl/test/NightlyRun/test2083.m (revision 21533) @@ -1,10 +1,10 @@ -%Test Name: GiaBenchmarksAB2dD3 +%Test Name: GiaIvinsBenchmarksAB2dD3 % Benchmark experiments (Figure A2d Ivins and James, 1999, Geophys. J. Int.) md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.); md=setmask(md,'',''); -md=parameterize(md,'../Par/GiaBenchmarksCD.par'); +md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par'); %% indicate what you want to compute -md.gia.cross_section_shape=2; % for elliptical edge +md.giaivins.cross_section_shape=2; % for elliptical edge %% define loading history @@ -21,11 +21,11 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('1111111'); -md=solve(md,'Gia'); +md=solve(md,'GiaIvins'); %Fields and tolerances to track changes -field_names ={'GiaW','GiadWdt'}; +field_names ={'GiaIvinsW','GiaIvinsdWdt'}; field_tolerances={1e-13,1e-13}; field_values={... - (md.results.GiaSolution.GiaW),... - (md.results.GiaSolution.GiadWdt),... + (md.results.GiaIvinsSolution.GiaIvinsW),... + (md.results.GiaIvinsSolution.GiaIvinsdWdt),... }; Index: sm/trunk-jpl/test/Par/GiaBenchmarksAB.par =================================================================== --- /issm/trunk-jpl/test/Par/GiaBenchmarksAB.par (revision 21532) +++ (revision ) @@ -1,77 +1,0 @@ -%Geometry specific to Experiments A and B -rad=600000; -nv=md.mesh.numberofvertices; -for i=1:nv, - dist=sqrt(md.mesh.x(i).^2+md.mesh.y(i)^2); - if (dist <= rad) - md.geometry.thickness(i)=2000.0; - else - md.geometry.thickness(i)=1.0; % non-zero thickness - end -end -md.geometry.thickness=md.geometry.thickness'; -md.geometry.base=zeros(md.mesh.numberofvertices,1); -md.geometry.surface=md.geometry.thickness+md.geometry.base; - -%Ice density used for benchmarking, not 917 kg/m^3 -md.materials.rho_ice=1000; %kg m^3 - -%GIA parameters specific to Experiments A and B -md.gia.mantle_viscosity=10^21*ones(md.mesh.numberofvertices,1); %in Pa.s -md.gia.lithosphere_thickness=100*ones(md.mesh.numberofvertices,1); %in km -md.materials.lithosphere_shear_modulus=6.7*10^10; %in Pa -md.materials.lithosphere_density=3.36; %in g/cm^3 -md.materials.mantle_shear_modulus=1.45*10^11; %in Pa -md.materials.mantle_density=3.38; %in g/cm^3 - - -%Initial velocity -x = archread('../Data/SquareSheetConstrained.arch','x'); -y = archread('../Data/SquareSheetConstrained.arch','y'); -vx = archread('../Data/SquareSheetConstrained.arch','vx'); -vy = archread('../Data/SquareSheetConstrained.arch','vy'); -index = archread('../Data/SquareSheetConstrained.arch','index'); - -x = x{1}; -y = y{1}; -vx = vx{1}; -vy = vy{1}; -index = index{1}; - -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); -clear vx vy x y index; -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); - -%Friction -md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); -md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.; -md.friction.p=ones(md.mesh.numberofelements,1); -md.friction.q=ones(md.mesh.numberofelements,1); - -%Numerical parameters -md.stressbalance.viscosity_overshoot=0.0; -md.masstransport.stabilization=1.; -md.thermal.stabilization=1.; -md.verbose=verbose(0); -md.settings.waitonlock=30; -md.stressbalance.restol=0.05; -md.steadystate.reltol=0.05; -md.stressbalance.reltol=0.05; -md.stressbalance.abstol=NaN; -md.timestepping.time_step=1.; -md.timestepping.final_time=3.; - -%Boundary conditions: -md=SetIceSheetBC(md); - -%Change name so that no test have the same name -A=dbstack; -if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end - Index: sm/trunk-jpl/test/Par/GiaBenchmarksCD.par =================================================================== --- /issm/trunk-jpl/test/Par/GiaBenchmarksCD.par (revision 21532) +++ (revision ) @@ -1,76 +1,0 @@ -%Geometry specific to Experiments C and D -rad=800000; -nv=md.mesh.numberofvertices; -for i=1:nv, - dist=sqrt(md.mesh.x(i).^2+md.mesh.y(i)^2); - if (dist <= rad) - md.geometry.thickness(i)=3000.0; - else - md.geometry.thickness(i)=1.0; % non-zero thickness - end -end -md.geometry.thickness=md.geometry.thickness'; -md.geometry.base=zeros(md.mesh.numberofvertices,1); -md.geometry.surface=md.geometry.thickness+md.geometry.base; - -%Ice density used for benchmarking, not 917 kg/m^3 -md.materials.rho_ice=1000; %kg m^3 - -%GIA parameters specific to Experiments A and B -md.gia.mantle_viscosity=10^21*ones(md.mesh.numberofvertices,1); %in Pa.s -md.gia.lithosphere_thickness=100*ones(md.mesh.numberofvertices,1); %in km -md.materials.lithosphere_shear_modulus=6.7*10^10; %in Pa -md.materials.lithosphere_density=3.32; %in g/cm^3 -md.materials.mantle_shear_modulus=1.45*10^11; %in Pa -md.materials.mantle_density=3.34; %in g/cm^3 - -%Initial velocity -x = archread('../Data/SquareSheetConstrained.arch','x'); -y = archread('../Data/SquareSheetConstrained.arch','y'); -vx = archread('../Data/SquareSheetConstrained.arch','vx'); -vy = archread('../Data/SquareSheetConstrained.arch','vy'); -index = archread('../Data/SquareSheetConstrained.arch','index'); - -x = x{1}; -y = y{1}; -vx = vx{1}; -vy = vy{1}; -index = index{1}; - -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); -clear vx vy x y index; -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); - -%Friction -md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); -md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.; -md.friction.p=ones(md.mesh.numberofelements,1); -md.friction.q=ones(md.mesh.numberofelements,1); - -%Numerical parameters -md.stressbalance.viscosity_overshoot=0.0; -md.masstransport.stabilization=1.; -md.thermal.stabilization=1.; -md.verbose=verbose(0); -md.settings.waitonlock=30; -md.stressbalance.restol=0.05; -md.steadystate.reltol=0.05; -md.stressbalance.reltol=0.05; -md.stressbalance.abstol=NaN; -md.timestepping.time_step=1.; -md.timestepping.final_time=3.; - -%Boundary conditions: -md=SetIceSheetBC(md); - -%Change name so that no test have the same name -A=dbstack; -if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end - Index: /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.par =================================================================== --- /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.par (revision 21533) +++ /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.par (revision 21533) @@ -0,0 +1,77 @@ +%Geometry specific to Experiments A and B +rad=600000; +nv=md.mesh.numberofvertices; +for i=1:nv, + dist=sqrt(md.mesh.x(i).^2+md.mesh.y(i)^2); + if (dist <= rad) + md.geometry.thickness(i)=2000.0; + else + md.geometry.thickness(i)=1.0; % non-zero thickness + end +end +md.geometry.thickness=md.geometry.thickness'; +md.geometry.base=zeros(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.thickness+md.geometry.base; + +%Ice density used for benchmarking, not 917 kg/m^3 +md.materials.rho_ice=1000; %kg m^3 + +%GIA parameters specific to Experiments A and B +md.giaivins.mantle_viscosity=10^21*ones(md.mesh.numberofvertices,1); %in Pa.s +md.giaivins.lithosphere_thickness=100*ones(md.mesh.numberofvertices,1); %in km +md.materials.lithosphere_shear_modulus=6.7*10^10; %in Pa +md.materials.lithosphere_density=3.36; %in g/cm^3 +md.materials.mantle_shear_modulus=1.45*10^11; %in Pa +md.materials.mantle_density=3.38; %in g/cm^3 + + +%Initial velocity +x = archread('../Data/SquareSheetConstrained.arch','x'); +y = archread('../Data/SquareSheetConstrained.arch','y'); +vx = archread('../Data/SquareSheetConstrained.arch','vx'); +vy = archread('../Data/SquareSheetConstrained.arch','vy'); +index = archread('../Data/SquareSheetConstrained.arch','index'); + +x = x{1}; +y = y{1}; +vx = vx{1}; +vy = vy{1}; +index = index{1}; + +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); +clear vx vy x y index; +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); + +%Friction +md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); +md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.; +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=ones(md.mesh.numberofelements,1); + +%Numerical parameters +md.stressbalance.viscosity_overshoot=0.0; +md.masstransport.stabilization=1.; +md.thermal.stabilization=1.; +md.verbose=verbose(0); +md.settings.waitonlock=30; +md.stressbalance.restol=0.05; +md.steadystate.reltol=0.05; +md.stressbalance.reltol=0.05; +md.stressbalance.abstol=NaN; +md.timestepping.time_step=1.; +md.timestepping.final_time=3.; + +%Boundary conditions: +md=SetIceSheetBC(md); + +%Change name so that no test have the same name +A=dbstack; +if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end + Index: /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.par =================================================================== --- /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.par (revision 21533) +++ /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.par (revision 21533) @@ -0,0 +1,76 @@ +%Geometry specific to Experiments C and D +rad=800000; +nv=md.mesh.numberofvertices; +for i=1:nv, + dist=sqrt(md.mesh.x(i).^2+md.mesh.y(i)^2); + if (dist <= rad) + md.geometry.thickness(i)=3000.0; + else + md.geometry.thickness(i)=1.0; % non-zero thickness + end +end +md.geometry.thickness=md.geometry.thickness'; +md.geometry.base=zeros(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.thickness+md.geometry.base; + +%Ice density used for benchmarking, not 917 kg/m^3 +md.materials.rho_ice=1000; %kg m^3 + +%GIA parameters specific to Experiments A and B +md.giaivins.mantle_viscosity=10^21*ones(md.mesh.numberofvertices,1); %in Pa.s +md.giaivins.lithosphere_thickness=100*ones(md.mesh.numberofvertices,1); %in km +md.materials.lithosphere_shear_modulus=6.7*10^10; %in Pa +md.materials.lithosphere_density=3.32; %in g/cm^3 +md.materials.mantle_shear_modulus=1.45*10^11; %in Pa +md.materials.mantle_density=3.34; %in g/cm^3 + +%Initial velocity +x = archread('../Data/SquareSheetConstrained.arch','x'); +y = archread('../Data/SquareSheetConstrained.arch','y'); +vx = archread('../Data/SquareSheetConstrained.arch','vx'); +vy = archread('../Data/SquareSheetConstrained.arch','vy'); +index = archread('../Data/SquareSheetConstrained.arch','index'); + +x = x{1}; +y = y{1}; +vx = vx{1}; +vy = vy{1}; +index = index{1}; + +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); +clear vx vy x y index; +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); + +%Friction +md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); +md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.; +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=ones(md.mesh.numberofelements,1); + +%Numerical parameters +md.stressbalance.viscosity_overshoot=0.0; +md.masstransport.stabilization=1.; +md.thermal.stabilization=1.; +md.verbose=verbose(0); +md.settings.waitonlock=30; +md.stressbalance.restol=0.05; +md.steadystate.reltol=0.05; +md.stressbalance.reltol=0.05; +md.stressbalance.abstol=NaN; +md.timestepping.time_step=1.; +md.timestepping.final_time=3.; + +%Boundary conditions: +md=SetIceSheetBC(md); + +%Change name so that no test have the same name +A=dbstack; +if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end + Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.par (revision 21532) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.par (revision 21533) @@ -60,6 +60,6 @@ %GIA: -md.gia.lithosphere_thickness=100.*ones(md.mesh.numberofvertices,1); % in km -md.gia.mantle_viscosity=1.0*10^21*ones(md.mesh.numberofvertices,1); % in Pa.s +md.giaivins.lithosphere_thickness=100.*ones(md.mesh.numberofvertices,1); % in km +md.giaivins.mantle_viscosity=1.0*10^21*ones(md.mesh.numberofvertices,1); % in Pa.s md.materials.lithosphere_shear_modulus=6.7*10^10; % in Pa md.materials.lithosphere_density=3.32; % in g/cm^-3 Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 21532) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 21533) @@ -62,6 +62,6 @@ #GIA: -md.gia.lithosphere_thickness=100.*numpy.ones((md.mesh.numberofvertices)); # in km -md.gia.mantle_viscosity=1.*10**21*numpy.ones((md.mesh.numberofvertices)); # in Pa.s +md.giaivins.lithosphere_thickness=100.*numpy.ones((md.mesh.numberofvertices)); # in km +md.giaivins.mantle_viscosity=1.*10**21*numpy.ones((md.mesh.numberofvertices)); # in Pa.s md.materials.lithosphere_shear_modulus=6.7*10**10; # in Pa md.materials.lithosphere_density=3.32; # in g/cm^-3