Index: /issm/trunk-jpl/test/SandBox/Makefile =================================================================== --- /issm/trunk-jpl/test/SandBox/Makefile (revision 25216) +++ /issm/trunk-jpl/test/SandBox/Makefile (revision 25217) @@ -1,2 +1,2 @@ clean: - rm -rf *.tar.gz *.errlog *.outlog *.bin *.qmu.in *.qmu.out *qmu.err *.queue *.toolkits + rm -rf *.tar.gz *.errlog *.outlog *.bin *.qmu.in *.qmu.out *qmu.err *.queue *.toolkits *.out Index: /issm/trunk-jpl/test/SandBox/bp2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/bp2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/bp2002.m (revision 25217) @@ -0,0 +1,74 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %500 km resolution mesh +md.cluster.np=3 +md.verbose=verbose('11111111'); + +%parameterize slr solution: +%slr loading: {{{ +md.slr.deltathickness=zeros(md.mesh.numberofelements,1); +md.slr.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + +%elastic loading from love numbers: +nlov=101; +md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(md.slr.deltathickness); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +pos=find(md.slr.deltathickness); +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.slr.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%New stuff +md.slr.spcthickness = NaN(md.mesh.numberofvertices,1); +md.slr.Ngia = zeros(md.mesh.numberofvertices,1); +md.slr.Ugia = zeros(md.mesh.numberofvertices,1); +md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.slr.maxiter=10; + +%bottom pressures: +pos=find(md.mesh.long<-150 & md.mesh.long > -160 & md.mesh.lat < 30 & md.mesh.lat > 10); +md.dsl.sea_water_pressure_change_at_sea_floor(pos)=1*md.constants.yts; %mm/yr +md.dsl.compute_fingerprints=1; + +%eustatic run: +md.slr.rigid=1; md.slr.elastic=1;md.slr.rotation=1; +md=solve(md,'Sealevelrise'); + Index: /issm/trunk-jpl/test/SandBox/corr2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/corr2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/corr2002.m (revision 25217) @@ -0,0 +1,182 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %500 km resolution mesh +md.cluster.np=4 + +%parameterize slr solution: +%slr loading: {{{ +nt=10; +md.slr.deltathickness=zeros(md.mesh.numberofelements+1,nt); +md.slr.deltathickness(end,:)=2001:1:2001+nt-1; +md.slr.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; +posa=find(late <-80); +md.slr.deltathickness(posa,:)=-100; + +%greenland +posg=find(late > 70 & late < 80 & longe>-60 & longe<-30); +md.slr.deltathickness(posg,:)=-100; + +%elastic loading from love numbers: +nlov=101; +md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(sum(md.slr.deltathickness(1:end-1,:),2)); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.slr.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%New stuff +md.slr.spcthickness = NaN(md.mesh.numberofvertices,1); +md.slr.Ngia = zeros(md.mesh.numberofvertices,1); +md.slr.Ugia = zeros(md.mesh.numberofvertices,1); +md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.slr.maxiter=10; + +%Uncertainty Quantification%{{{ + + %partitioning + npart=2; + partition=-ones(md.mesh.numberofelements,1); + partition(posg)=0; + partition(posa)=1; + + %variables: + md.qmu.variables.deltathickness=normal_uncertain('descriptor','scaled_SealevelriseDeltathickness',... + 'mean',ones(npart,nt),... + 'stddev',.2*ones(npart,nt),... + 'partition',partition,'nsteps',nt); + + md.qmu.correlation_matrix=zeros(npart*nt,npart*nt); + for i=1:npart, + for j=1:nt, + indi=(i-1)*nt+j; + for k=1:npart, + for l=1:nt, + indj=(k-1)*nt+l; + if i~=k, + md.qmu.correlation_matrix(indi,indj)=0; + else + %same partition: + if j==l, + md.qmu.correlation_matrix(indi,indj)=1; + else + if i==1, + md.qmu.correlation_matrix(indi,indj)=.2; + else + md.qmu.correlation_matrix(indi,indj)=.5; + end + end + end + end + end + end + end + + + %responses + md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1'); + md.qmu.responses.sealevel2=response_function('descriptor','Outputdefinition2'); + md.qmu.responses.sealevel3=response_function('descriptor','Outputdefinition3'); + md.qmu.responses.sealevel4=response_function('descriptor','Outputdefinition4'); + md.qmu.responses.sealevel5=response_function('descriptor','Outputdefinition5'); + md.qmu.responses.sealevel6=response_function('descriptor','Outputdefinition6'); + md.qmu.responses.sealevel7=response_function('descriptor','Outputdefinition7'); + md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition8'); + md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition9'); + md.qmu.responses.sealevel10=response_function('descriptor','Outputdefinition10'); + + %output definitions: + for i=1:10, + if i==1, + md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ... + 'model_string','Sealevel','node',i)}; + else + md.outputdefinition.definitions{end+1}=nodalvalue('name','SNode','definitionstring',['Outputdefinition' num2str(i)], ... + 'model_string','Sealevel','node',i); + end + end + %algorithm: + md.qmu.method =dakota_method('nond_samp'); + md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... + 'seed',1234,... + 'samples',40,... + 'sample_type','lhs'); + + %parameters + md.qmu.params.direct=true; + md.qmu.params.interval_type='forward'; + md.qmu.params.analysis_driver='matlab'; + md.qmu.params.evaluation_scheduling='master'; + md.qmu.params.processors_per_evaluation=3; + md.qmu.params.tabular_graphics_data=true; + %}}} +md.qmu.isdakota=1; +md.qmu.output=1; + +%transient: +md.timestepping.start_time=2000; +md.timestepping.time_step=1; +md.timestepping.interp_forcings=0; +md.timestepping.final_time=2011; +md.transient.issmb=0; +md.transient.ismasstransport=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +md.transient.isslr=1; +md.slr.requested_outputs= {'default',... + 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'SealevelNGiaRate', 'SealevelUGiaRate',... + 'SealevelEustaticMask','SealevelEustaticOceanMask'}; + +%hack: +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.base=-ones(md.mesh.numberofvertices,1); +md.geometry.surface=zeros(md.mesh.numberofvertices,1); + +%eustatic + rigid + elastic + rotation run: +md.verbose=verbose('0'); +md.verbose.qmu=1; +md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1; +md=solve(md,'tr'); + Index: /issm/trunk-jpl/test/SandBox/dak2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/dak2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/dak2002.m (revision 25217) @@ -0,0 +1,194 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %500 km resolution mesh +md.cluster.np=2 + +%parameterize slr solution: +%slr loading: {{{ +nt=100; +md.slr.deltathickness=zeros(md.mesh.numberofelements+1,nt); +md.slr.deltathickness(end,:)=2000.5:1:2099.5; +md.slr.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; +posa=find(late <-80); +md.slr.deltathickness(posa,:)=-100; +%greenland +posg=find(late > 70 & late < 80 & longe>-60 & longe<-30); +md.slr.deltathickness(posg,:)=-100; + + +%elastic loading from love numbers: +nlov=101; +md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(sum(md.slr.deltathickness(1:end-1,:),2)); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.slr.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%New stuff +md.slr.spcthickness = NaN(md.mesh.numberofvertices,1); +md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.slr.maxiter=10; + +%initialize GIA: +md.gia=giamme(); +md.gia.Ngia=ones(md.mesh.numberofvertices,10); +md.gia.Ugia=ones(md.mesh.numberofvertices,10); +for i=1:10, + md.gia.Ngia(i)=i; + md.gia.Ugia(i)=i; +end +md.gia.modelid=1; + +%Uncertainty Quantification%{{{ + md.qmu.variables=struct();; + + ns=size(md.gia.Ngia,2); + ids=(1:(ns+1))'; + p=rand(ns,1); + probs=[p(1:ns); 0]; probs=probs/sum(probs); + md.qmu.variables.giamodelid=histogram_bin_uncertain('GiaModelid',ns+1,ids,probs); + + %partitioning + npart=2; + partition=-ones(md.mesh.numberofelements,1); + partition(posg)=0; + partition(posa)=1; + + %variables: + md.qmu.variables.deltathickness=normal_uncertain('descriptor','scaled_SealevelriseDeltathickness',... + 'mean',ones(2,nt),... + 'stddev',[.1*ones(1,nt);.2*ones(1,nt)],... + 'partition',partition,'nsteps',nt); + + md.qmu.correlation_matrix=zeros(npart*nt,npart*nt); + for i=1:npart, + for j=1:nt, + indi=(i-1)*nt+j; + for k=1:npart, + for l=1:nt, + indj=(k-1)*nt+l; + if i~=k, + md.qmu.correlation_matrix(indi,indj)=0; + else + %same partition: + if j==l, + md.qmu.correlation_matrix(indi,indj)=1; + else + md.qmu.correlation_matrix(indi,indj)=.2; + end + end + end + end + end + end + md.qmu.correlation_matrix=[]; + + + + %responses + md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1'); + md.qmu.responses.sealevel2=response_function('descriptor','Outputdefinition2'); + md.qmu.responses.sealevel3=response_function('descriptor','Outputdefinition3'); + md.qmu.responses.sealevel4=response_function('descriptor','Outputdefinition4'); + md.qmu.responses.sealevel5=response_function('descriptor','Outputdefinition5'); + md.qmu.responses.sealevel6=response_function('descriptor','Outputdefinition6'); + md.qmu.responses.sealevel7=response_function('descriptor','Outputdefinition7'); + md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition8'); + md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition9'); + md.qmu.responses.sealevel10=response_function('descriptor','Outputdefinition10'); + + %output definitions: + for i=1:10, + if i==1, + md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ... + 'model_string','Sealevel','node',i)}; + else + md.outputdefinition.definitions{end+1}=nodalvalue('name','SNode','definitionstring',['Outputdefinition' num2str(i)], ... + 'model_string','Sealevel','node',i); + end + end + %algorithm: + md.qmu.method =dakota_method('nond_samp'); + md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... + 'seed',1234,... + 'samples',3,... + 'sample_type','lhs'); + + %parameters + md.qmu.params.direct=true; + md.qmu.params.interval_type='forward'; + md.qmu.params.analysis_driver='matlab'; + md.qmu.params.evaluation_scheduling='master'; + md.qmu.params.processors_per_evaluation=1; + md.qmu.params.tabular_graphics_data=true; + %}}} +md.qmu.output=1; + +%transient: +md.timestepping.start_time=2000; +md.timestepping.interp_forcings=0; +md.timestepping.final_time=2002; +md.transient.issmb=0; +md.transient.ismasstransport=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +md.transient.isslr=1; +md.transient.isgia=1; +md.slr.requested_outputs= {'default',... + 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate',... + 'SealevelEustaticMask','SealevelEustaticOceanMask'}; + +%hack: +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.base=-ones(md.mesh.numberofvertices,1); +md.geometry.surface=zeros(md.mesh.numberofvertices,1); + +%eustatic + rigid + elastic + rotation run: +md.verbose=verbose('11111111111'); +%md.verbose.qmu=1; +md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1; +md.qmu.isdakota=1; +md=solve(md,'slr'); + Index: /issm/trunk-jpl/test/SandBox/geoid2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/geoid2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/geoid2002.m (revision 25217) @@ -0,0 +1,80 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %500 km resolution mesh +md.cluster.np=3 +md.verbose=verbose('11111111'); + +%parameterize slr solution: +%slr loading: {{{ +md.slr.deltathickness=zeros(md.mesh.numberofelements,1); +md.slr.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; +pos=find(late <-80); +md.slr.deltathickness(pos)=-100; +%greenland +pos=find(late > 70 & late < 80 & longe>-60 & longe<-30); +md.slr.deltathickness(pos)=-100; + +%elastic loading from love numbers: +nlov=101; +md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(md.slr.deltathickness); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +pos=find(md.slr.deltathickness); +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.slr.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%New stuff +md.slr.spcthickness = NaN(md.mesh.numberofvertices,1); +md.slr.Ngia = zeros(md.mesh.numberofvertices,1); +md.slr.Ugia = zeros(md.mesh.numberofvertices,1); +md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.slr.maxiter=10; + +md.slr.requested_outputs= {'default',... + 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'SealevelNGiaRate', 'SealevelUGiaRate',... + 'SealevelEustaticMask','SealevelEustaticOceanMask','SealevelUEsa','SealevelNEsa'}; + +md.slr.rigid=1; md.slr.elastic=1;md.slr.rotation=1; +md=solve(md,'Sealevelrise'); Index: /issm/trunk-jpl/test/SandBox/gia2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/gia2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/gia2002.m (revision 25217) @@ -0,0 +1,151 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1500.); %500 km resolution mesh +md.cluster.np=16; + +%parameterize slr solution: +%slr loading: {{{ +md.slr.deltathickness=zeros(md.mesh.numberofelements,1); +md.slr.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; +pos=find(late <-80); +md.slr.deltathickness(pos)=-100; +%greenland +pos=find(late > 70 & late < 80 & longe>-60 & longe<-30); +md.slr.deltathickness(pos)=-100; + +%elastic loading from love numbers: +nlov=101; +md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(md.slr.deltathickness); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +pos=find(md.slr.deltathickness); +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.slr.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%New stuff +md.slr.spcthickness = NaN(md.mesh.numberofvertices,1); +md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); + +%GIA: +md.gia=giamme(); +md.gia.Ngia= ones(md.mesh.numberofvertices,10); +md.gia.Ugia= ones(md.mesh.numberofvertices,10); +load post; +ns=length(p); +for i=1:ns, + md.gia.Ngia(:,i)=i; + md.gia.Ugia(:,i)=i; +end +md.gia.modelid=10; +md.verbose=verbose(0); +md.verbose.qmu=1; + +md.slr.requested_outputs= {'default',... + 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate',... + 'SealevelEustaticMask','SealevelEustaticOceanMask'}; + +ids=(1:(ns+1))'; +probs=[p; 0]; +%probs=[(1:ns)'; 0]; probs=probs/sum(probs); +md.qmu.variables.giamodelid=histogram_bin_uncertain('GiaModelid',ns+1,ids,probs); + +%qmu: +%x=0:.1:100; +%y=lognormal_pdf(x2,0,1); +%md.qmu.variables.surface_mass_balance=histogram_bin_uncertain('scaled_SmbMassBalance',length(x),x,[y 0]); + +%x=1:(ns+1); +%y=[lognormal_pdf(x(1:end-1)+.5,0,1) 0]; y=y/sum(y); +%md.qmu.variables.giamodelid=histogram_bin_uncertain('GiaModelid',length(x),x,y); + +%md.qmu.variables.giamodelid=uniform_uncertain('descriptor','GiaModelid','lower',1,... +% 'upper',10); + +%responses +md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1'); +md.qmu.responses.sealevel2=response_function('descriptor','Outputdefinition2'); +md.qmu.responses.sealevel3=response_function('descriptor','Outputdefinition3'); +md.qmu.responses.sealevel4=response_function('descriptor','Outputdefinition4'); +md.qmu.responses.sealevel5=response_function('descriptor','Outputdefinition5'); +md.qmu.responses.sealevel6=response_function('descriptor','Outputdefinition6'); +md.qmu.responses.sealevel7=response_function('descriptor','Outputdefinition7'); +md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition8'); +md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition9'); +md.qmu.responses.sealevel10=response_function('descriptor','Outputdefinition10'); + +%output definitions: +for i=1:10, + if i==1, + md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ... + 'model_string','Sealevel','node',i)}; + else + md.outputdefinition.definitions{end+1}=nodalvalue('name','SNode','definitionstring',['Outputdefinition' num2str(i)], ... + 'model_string','Sealevel','node',i); + end +end + +%algorithm: +md.qmu.method =dakota_method('nond_samp'); +md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... +'seed',1234,... +'samples',50,... +'sample_type','lhs'); + +%parameters +md.qmu.params.direct=true; +md.qmu.params.interval_type='forward'; +md.qmu.params.analysis_driver='matlab'; +md.qmu.params.evaluation_scheduling='master'; +md.qmu.params.processors_per_evaluation=2; +md.qmu.params.tabular_graphics_data=true; + +md.qmu.isdakota=1; +md.qmu.output=1; + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.slr.maxiter=10; + +%eustatic run: +md.slr.rigid=1; md.slr.elastic=1;md.slr.rotation=1; +md=solve(md,'Sealevelrise'); + Index: /issm/trunk-jpl/test/SandBox/hier2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/hier2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/hier2002.m (revision 25217) @@ -0,0 +1,209 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %500 km resolution mesh +md.cluster.np=16; + +%parameterize slr solution: +%slr loading: {{{ +nt=100; +md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements+1,nt); +md.solidearth.surfaceload.icethicknesschange(end,:)=2000.5:1:2099.5; +md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; + +posant=find(late <-80); +md.solidearth.surfaceload.icethicknesschange(posant,:)=-50; + +%greenland +posgre=find(late > 70 & late < 80 & longe>-60 & longe<-30); +md.solidearth.surfaceload.icethicknesschange(posgre,:)=-100; + +%alaska : +posala=find(late > 62 & late < 68 & longe>-152 & longe<-147); +md.solidearth.surfaceload.icethicknesschange(posala,:)=-150; + +%himalaya : +poshim=find(late > 25 & late < 32 & longe>81 & longe<86); +md.solidearth.surfaceload.icethicknesschange(poshim,:)=-200; + + +%elastic loading from love numbers: +md.solidearth.lovenumbers=lovenumbers('maxdeg',100); + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(sum(md.solidearth.surfaceload.icethicknesschange(1:end-1,:),2)); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.solidearth.settings.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +%Solution parameters +md.solidearth.settings.reltol=NaN; +md.solidearth.settings.abstol=1e-3; +md.solidearth.settings.computesealevelchange=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.solidearth.settings.maxiter=10; + +%initialize GIA: +md.gia=giamme(); +md.gia.Ngia=ones(md.mesh.numberofvertices,1); +md.gia.Ugia=ones(md.mesh.numberofvertices,1); +md.gia.modelid=1; + +%Uncertainty Quantification%{{{ +md.qmu.variables=struct();; + +%partition vector: +npart=4; +partition=-ones(md.mesh.numberofelements,1); +partition(posgre)=0; +partition(posant)=1; +partition(posala)=2; +partition(poshim)=3; + +%prepare arrays: +background=md.solidearth.surfaceload.icethicknesschange; +antarray={}; for i=1:10, antarray{end+1}= (-50+i)*ones(length(posant),nsteps); end +grearray={}; for i=1:20, grearray{end+1}= (-100+i)*ones(length(posgre),nsteps); end +alaarray={}; for i=1:30, alaarray{end+1}= (-150+i)*ones(length(posala),nsteps); end +himarray={-200*ones(length(poshim),nsteps)}; +arrays={antarray,grearray,alaarray,himarray}; + +partarray={antarray,grearray,alaarray,himarray}; +histarray={equiprobable_histogram_uncertain('ant',10), equiprobable_histogram_uncertain('gre',20), equiprobable_histogram_uncertain('ala',30), equiprobable_histogram_uncertain('him',1)}; +pdfant=normal_uncertain('descriptor','pdfant','mean',1*ones(1,nt),'stddev',.1*ones(1,nt)); +pdfgre=uniform_uncertain('descriptor','pdfgre','lower',.8*ones(1,nt),'upper',1.2*ones(1,nt)); +pdfala=dirac_uncertain('descriptor','pdfala',2); +pdfhim=dirac_uncertain('descriptor','pdfhim',1); + +%variables: +md.qmu.variables.deltathickness=normal_uncertain('descriptor','SurfaceloadIceThicknessChange',... + 'partition',partition,'nsteps',nt,... + 'background',md.solidearth.surfaceload.icethicknesschange,... + 'histogram_bin_uncertain_perpartition',histarray,... + 'matarray_perpartition',partarray,... + + ); + +md.qmu.correlation_matrix=zeros(npart*nt,npart*nt); +for i=1:npart, + for j=1:nt, + indi=(i-1)*nt+j; + for k=1:npart, + for l=1:nt, + indj=(k-1)*nt+l; + if i~=k, + md.qmu.correlation_matrix(indi,indj)=0; + else + %same partition: + if j==l, + md.qmu.correlation_matrix(indi,indj)=1; + else + md.qmu.correlation_matrix(indi,indj)=.2; + end + end + end + end + end +end +md.qmu.correlation_matrix=[]; + + + +%responses +md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1'); +md.qmu.responses.sealevel2=response_function('descriptor','Outputdefinition2'); +md.qmu.responses.sealevel3=response_function('descriptor','Outputdefinition3'); +md.qmu.responses.sealevel4=response_function('descriptor','Outputdefinition4'); +md.qmu.responses.sealevel5=response_function('descriptor','Outputdefinition5'); +md.qmu.responses.sealevel6=response_function('descriptor','Outputdefinition6'); +md.qmu.responses.sealevel7=response_function('descriptor','Outputdefinition7'); +md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition8'); +md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition9'); +md.qmu.responses.sealevel10=response_function('descriptor','Outputdefinition10'); + +%output definitions: +for i=1:10, + if i==1, + md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ... + 'model_string','Sealevel','node',i)}; + else + md.outputdefinition.definitions{end+1}=nodalvalue('name','SNode','definitionstring',['Outputdefinition' num2str(i)], ... + 'model_string','Sealevel','node',i); +end +end +%algorithm: +md.qmu.method =dakota_method('nond_samp'); +md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... + 'seed',1234,... + 'samples',3,... + 'sample_type','lhs'); + +%parameters +md.qmu.params.direct=true; +md.qmu.params.interval_type='forward'; +md.qmu.params.analysis_driver='matlab'; +md.qmu.params.evaluation_scheduling='master'; +md.qmu.params.processors_per_evaluation=1; +md.qmu.params.tabular_graphics_data=true; +md.qmu.output=1; +%}}} + +%transient: +md.timestepping.start_time=2000; +md.timestepping.interp_forcings=0; +md.timestepping.final_time=2002; +md.transient.issmb=0; +md.transient.ismasstransport=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +md.transient.isslr=1; +md.transient.isgia=1; +md.slr.requested_outputs= {'default',... + 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate',... + 'SealevelEustaticMask','SealevelEustaticOceanMask'}; + +%hack: +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.base=-ones(md.mesh.numberofvertices,1); +md.geometry.surface=zeros(md.mesh.numberofvertices,1); + +%eustatic + rigid + elastic + rotation run: +md.verbose=verbose('11111111111'); +%md.verbose.qmu=1; +md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1;md.solidearth.settings.rotation=1; +md.qmu.isdakota=1; +md=solve(md,'slr'); + Index: /issm/trunk-jpl/test/SandBox/interpBedmap2.m =================================================================== --- /issm/trunk-jpl/test/SandBox/interpBedmap2.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/interpBedmap2.m (revision 25217) @@ -0,0 +1,47 @@ +function [dataout] = interpBedmap2(X,Y,string), +%INTERPBEDMAP2 - interpolate bedmap2 data +% +% Available data: +% 1. bed is bed height +% 2. surface is surface height +% 3. thickness is ice thickness +% 4. icemask_grounded_and_shelves is a mask file showing the grounding line and the extent of the floating ice shelves +% 5. rockmask is a mask file showing rock outcrops +% 6. lakemask_vostok is a mask file showing the extent of the lake cavity of Lake Vostok +% 7. bed_uncertainty is the bed uncertainty grid shown in figure 12 of the manuscript +% 8. thickness_uncertainty_5km is the thickness uncertainty grid shown in figure 11 of the manuscript +% 9. data_coverage is a binary grid showing the dis tribution of ice thickness data used in the grid of ice thickness +% +% Usage: +% [dataout] = interpBedmap2(X,Y,string) + +%reqad data +%path = '/u/astrid-r1b/ModelData/BedMap2/bedmap2_bin/'; +%path = '~/issm-jpl/proj-group/ModelData/BedMap2/bedmap2_bin/'; +path = '/Users/larour/ModelData/BedMap2/bedmap2_bin/'; +fid=fopen([path '/bedmap2_' string '.flt'],'r','l'); +data=fread(fid,[6667,6667],'float32'); +fclose(fid); + +% define grid +if strcmp(string,'thickness_uncertainty_5km'), + ncols =1361; + nrows =1361; + xll =-3401000; + yll =-3402000; + gridsize =5000; +else + ncols =6667; + nrows =6667; + xll =-3333000; + yll =-3333000; + gridsize =1000; +end +x_m=xll+(0:1:ncols-1)'*gridsize; +y_m=yll+(0:1:nrows-1)'*gridsize; + +%Change default to NaN +data(find(data==-9999))=0; + +%Interpolate +dataout = InterpFromGrid(x_m,y_m,flipud(data'),double(X),double(Y)); Index: /issm/trunk-jpl/test/SandBox/mmedak.m =================================================================== --- /issm/trunk-jpl/test/SandBox/mmedak.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/mmedak.m (revision 25217) @@ -0,0 +1,195 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %500 km resolution mesh +md.cluster.np=16; + +%parameterize slr solution: +%slr loading: {{{ +nt=100; +md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements+1,nt); +md.solidearth.surfaceload.icethicknesschange(end,:)=2000.5:1:2099.5; +md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; + +posant=find(late <-80); +md.solidearth.surfaceload.icethicknesschange(posant,:)=-50; + +%greenland +posgre=find(late > 70 & late < 80 & longe>-60 & longe<-30); +md.solidearth.surfaceload.icethicknesschange(posgre,:)=-100; + +%alaska : +posala=find(late > 62 & late < 68 & longe>-152 & longe<-147); +md.solidearth.surfaceload.icethicknesschange(posala,:)=-150; + +%himalaya : +poshim=find(late > 25 & late < 32 & longe>81 & longe<86); +md.solidearth.surfaceload.icethicknesschange(poshim,:)=-150; + + +%elastic loading from love numbers: +md.solidearth.lovenumbers=lovenumbers('maxdeg',100); + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(sum(md.solidearth.surfaceload.icethicknesschange(1:end-1,:),2)); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.solidearth.settings.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +%Solution parameters +md.solidearth.settings.reltol=NaN; +md.solidearth.settings.abstol=1e-3; +md.solidearth.settings.computesealevelchange=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.solidearth.settings.maxiter=10; + +%initialize GIA: +md.gia=giamme(); +md.gia.Ngia=ones(md.mesh.numberofvertices,1); +md.gia.Ugia=ones(md.mesh.numberofvertices,1); +md.gia.modelid=1; + +%Uncertainty Quantification%{{{ +md.qmu.variables=struct();; + +ns=size(md.gia.Ngia,2); +ids=(1:(ns+1))'; +p=rand(ns,1); +probs=[p(1:ns); 0]; probs=probs/sum(probs); +md.qmu.variables.giamodelid=histogram_bin_uncertain('GiaModelid',ns+1,ids,probs); + +%partitioning +npart=2; +partition=-ones(md.mesh.numberofelements,1); +partition(posg)=0; +partition(posa)=1; + +%variables: +md.qmu.variables.deltathickness=normal_uncertain('descriptor','scaled_SealevelriseDeltathickness',... + 'mean',ones(2,nt),... + 'stddev',[.1*ones(1,nt);.2*ones(1,nt)],... + 'partition',partition,'nsteps',nt); + +md.qmu.correlation_matrix=zeros(npart*nt,npart*nt); +for i=1:npart, + for j=1:nt, + indi=(i-1)*nt+j; + for k=1:npart, + for l=1:nt, + indj=(k-1)*nt+l; + if i~=k, + md.qmu.correlation_matrix(indi,indj)=0; + else + %same partition: + if j==l, + md.qmu.correlation_matrix(indi,indj)=1; + else + md.qmu.correlation_matrix(indi,indj)=.2; + end + end + end + end + end +end +md.qmu.correlation_matrix=[]; + + + +%responses +md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1'); +md.qmu.responses.sealevel2=response_function('descriptor','Outputdefinition2'); +md.qmu.responses.sealevel3=response_function('descriptor','Outputdefinition3'); +md.qmu.responses.sealevel4=response_function('descriptor','Outputdefinition4'); +md.qmu.responses.sealevel5=response_function('descriptor','Outputdefinition5'); +md.qmu.responses.sealevel6=response_function('descriptor','Outputdefinition6'); +md.qmu.responses.sealevel7=response_function('descriptor','Outputdefinition7'); +md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition8'); +md.qmu.responses.sealevel8=response_function('descriptor','Outputdefinition9'); +md.qmu.responses.sealevel10=response_function('descriptor','Outputdefinition10'); + +%output definitions: +for i=1:10, + if i==1, + md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ... + 'model_string','Sealevel','node',i)}; + else + md.outputdefinition.definitions{end+1}=nodalvalue('name','SNode','definitionstring',['Outputdefinition' num2str(i)], ... + 'model_string','Sealevel','node',i); +end +end +%algorithm: +md.qmu.method =dakota_method('nond_samp'); +md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... + 'seed',1234,... + 'samples',3,... + 'sample_type','lhs'); + +%parameters +md.qmu.params.direct=true; +md.qmu.params.interval_type='forward'; +md.qmu.params.analysis_driver='matlab'; +md.qmu.params.evaluation_scheduling='master'; +md.qmu.params.processors_per_evaluation=1; +md.qmu.params.tabular_graphics_data=true; +md.qmu.output=1; +%}}} + +%transient: +md.timestepping.start_time=2000; +md.timestepping.interp_forcings=0; +md.timestepping.final_time=2002; +md.transient.issmb=0; +md.transient.ismasstransport=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +md.transient.isslr=1; +md.transient.isgia=1; +md.slr.requested_outputs= {'default',... + 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate',... + 'SealevelEustaticMask','SealevelEustaticOceanMask'}; + +%hack: +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.base=-ones(md.mesh.numberofvertices,1); +md.geometry.surface=zeros(md.mesh.numberofvertices,1); + +%eustatic + rigid + elastic + rotation run: +md.verbose=verbose('11111111111'); +%md.verbose.qmu=1; +md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1;md.solidearth.settings.rotation=1; +md.qmu.isdakota=1; +md=solve(md,'slr'); + Index: /issm/trunk-jpl/test/SandBox/res2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/res2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/res2002.m (revision 25217) @@ -0,0 +1,8 @@ +ha=[]; +for i=1:md.qmu.method.params.samples, + md2.results.TransientSolution=md.results.dakota.modelresults{i}.TransientSolution; + h=resultstomatrix(md2,'TransientSolution','SealevelriseDeltathickness'); + ha=[ha;h(posa(1),:)]; +end + +plot(corr(ha(:,1),ha)); Index: /issm/trunk-jpl/test/SandBox/testbayes.m =================================================================== --- /issm/trunk-jpl/test/SandBox/testbayes.m (revision 25216) +++ /issm/trunk-jpl/test/SandBox/testbayes.m (revision 25217) @@ -5,11 +5,10 @@ md=setflowequation(md,'SSA','all'); md.cluster=generic('name',oshostname(),'np',3); -md.materials.rho_ice=10^7; %involved in the mass flux, make it easy -md.geometry.thickness(:)=1; %make it easy -md.geometry.surface=md.geometry.base+md.geometry.thickness; %constrain all velocities to 1 m/yr, in the y-direction -md.stressbalance.spcvx(:)=0; -md.stressbalance.spcvy(:)=1; +pos=find(md.mesh.y==0); +md.stressbalance.spcvx(pos)=0; +md.stressbalance.spcvy(pos)=0; +%md.stressbalance.spcvy(:)=1; md.stressbalance.spcvz(:)=0; @@ -24,5 +23,5 @@ %input/output: -md.qmu.variables.drag_coefficient=normal_uncertain('scaled_FrictionCoefficient',1,0.01); +md.qmu.variables.drag_coefficient=normal_uncertain('scaled_FrictionCoefficient',1,0.10); md.qmu.responses.misfit=calibration_function('MaxVel'); %md.qmu.responses.MaxVel=response_function('MaxVel',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]); @@ -32,5 +31,5 @@ md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... 'seed',1234,... -'samples',20,... +'samples',2000,... 'queso',true,... 'metropolis_hastings',true,... @@ -53,22 +52,2 @@ %Fields and tolerances to track changes md.qmu.results=md.results.dakota; - -error; - -%ok, mass flux of 3 profiles should be -3 Gt/yr -3 Gt/yr and the sum, which is -6 Gt/yr -%we recover those mass fluxes through the mean of the response. -%also, we recover the max velo, which should be 1m/yr. -%we put all that data in the montecarlo field, which we will use to test for success. -%also, check that the stddev are 0. -md.results.dakota.montecarlo=[]; -for i=1:8, - md.results.dakota.montecarlo=[md.results.dakota.montecarlo md.results.dakota.dresp_out(i).mean]; -end -for i=1:8, - md.results.dakota.montecarlo=[md.results.dakota.montecarlo md.results.dakota.dresp_out(i).stddev]; -end -field_names ={'montecarlo'}; -field_tolerances={1e-11}; -field_values={... - md.results.dakota.montecarlo,... - }; Index: /issm/trunk-jpl/test/SandBox/tran2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/tran2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/tran2002.m (revision 25217) @@ -0,0 +1,104 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %500 km resolution mesh +md.cluster.np=16 + +%parameterize slr solution: +%slr loading: {{{ +md.slr.deltathickness=zeros(md.mesh.numberofelements,1); +md.slr.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; +pos=find(late <-80); +md.slr.deltathickness(pos)=-100; +%greenland +pos=find(late > 70 & late < 80 & longe>-60 & longe<-30); +md.slr.deltathickness(pos)=-100; + +%elastic loading from love numbers: +nlov=101; +md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(md.slr.deltathickness); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +pos=find(md.slr.deltathickness); +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.slr.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%New stuff +md.slr.spcthickness = NaN(md.mesh.numberofvertices,1); +md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.slr.maxiter=10; + +%transient settings: +md.timestepping.start_time=0; +md.timestepping.final_time=10; +md.timestepping.time_step=1; +md.transient.isslr=1; +md.transient.issmb=0; +md.transient.isgia=0; +md.transient.ismasstransport=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +dh=md.slr.deltathickness; +deltathickness=zeros(md.mesh.numberofelements+1,10); +for i=1:10, + deltathickness(1:end-1,i)=dh*i; +end +deltathickness(end,:)=0:1:9; +md.slr.deltathickness=deltathickness; + +md.gia=giamme; +md.gia.Ngia=zeros(md.mesh.numberofvertices,1); +md.gia.Ugia=zeros(md.mesh.numberofvertices,1); +md.gia.modelid=1; + + +%dummy: +md.geometry.surface=zeros(md.mesh.numberofvertices,1); +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.base=-ones(md.mesh.numberofvertices,1); +md.geometry.bed=md.geometry.base; + +%eustatic + rigid + elastic + rotation run: +md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1; +md.verbose=verbose('1111111'); +md=solve(md,'slr'); Index: /issm/trunk-jpl/test/SandBox/val2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/val2002.m (revision 25217) +++ /issm/trunk-jpl/test/SandBox/val2002.m (revision 25217) @@ -0,0 +1,86 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %500 km resolution mesh + +%parameterize slr solution: +%slr loading: {{{ +md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); +md.solidearth.surfaceload.waterheightchange=zeros(md.mesh.numberofelements,1); +md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +%US +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; + +pos=find(late < (37+5) & late > (37-5) & longe>(-95-10) & longe<(-95+10)); +md.solidearth.surfaceload.waterheightchange(pos)=-100; + +%elastic loading from love numbers: +md.solidearth.lovenumbers=lovenumbers('maxdeg',100); + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(md.solidearth.surfaceload.waterheightchange); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; +md.mask.ice_levelset(md.mesh.elements(pos,:))=1; + +% }}} + +md.solidearth.settings.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + +%Hydro load: + + +%Miscellaneous +md.miscellaneous.name='test2002'; + +%Solution parameters +md.solidearth.settings.reltol=NaN; +md.solidearth.settings.abstol=1e-3; +md.solidearth.settings.computesealevelchange=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.solidearth.settings.maxiter=10; + +%eustatic run: +md.solidearth.settings.rigid=0; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0; +md=solve(md,'Sealevelrise'); +Seustatic=md.results.SealevelriseSolution.Sealevel; + +%eustatic + rigid run: +md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0; +md=solve(md,'Sealevelrise'); +Srigid=md.results.SealevelriseSolution.Sealevel; + +%eustatic + rigid + elastic run: +md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1;md.solidearth.settings.rotation=0; +md=solve(md,'Sealevelrise'); +Selastic=md.results.SealevelriseSolution.Sealevel; + +%eustatic + rigid + elastic + rotation run: +md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; +md=solve(md,'Sealevelrise'); +Srotation=md.results.SealevelriseSolution.Sealevel; + +%Fields and tolerances to track changes +field_names ={'Eustatic','Rigid','Elastic','Rotation'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13}; +field_values={Seustatic,Srigid,Selastic,Srotation};