Index: /issm/trunk-jpl/test/SandBox/Makefile =================================================================== --- /issm/trunk-jpl/test/SandBox/Makefile (revision 25267) +++ /issm/trunk-jpl/test/SandBox/Makefile (revision 25268) @@ -1,2 +1,2 @@ clean: - rm -rf *.tar.gz *.errlog *.outlog *.bin *.qmu.in *.qmu.out *qmu.err *.queue *.toolkits *.out + rm -rf *.tar.gz *.errlog *.outlog *.bin *.qmu.in *.qmu.out *qmu.err *.queue *.toolkits *.out *.outbin.* Index: /issm/trunk-jpl/test/SandBox/hier2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/hier2002.m (revision 25267) +++ /issm/trunk-jpl/test/SandBox/hier2002.m (revision 25268) @@ -3,11 +3,11 @@ %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %500 km resolution mesh -md.cluster.np=16; +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.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements+1,nt); +nsteps=100; +md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements+1,nsteps); md.solidearth.surfaceload.icethicknesschange(end,:)=2000.5:1:2099.5; md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); @@ -24,15 +24,15 @@ %greenland -posgre=find(late > 70 & late < 80 & longe>-60 & longe<-30); +posgre=find(late > 70 & late < 80 & longe>-50 & longe<-30); md.solidearth.surfaceload.icethicknesschange(posgre,:)=-100; %alaska : -posala=find(late > 62 & late < 68 & longe>-152 & longe<-147); +posala=find(late > 55 & late < 68 & longe>-162 & longe<-140); md.solidearth.surfaceload.icethicknesschange(posala,:)=-150; -%himalaya : -poshim=find(late > 25 & late < 32 & longe>81 & longe<86); -md.solidearth.surfaceload.icethicknesschange(poshim,:)=-200; - + +%hawaii : +poshaw=find(late > 15 & late < 25 & longe>-170 & longe<-130); +md.solidearth.surfaceload.icethicknesschange(poshaw,:)=10; %elastic loading from love numbers: @@ -78,67 +78,99 @@ %initialize GIA: md.gia=giamme(); -md.gia.Ngia=ones(md.mesh.numberofvertices,1); -md.gia.Ugia=ones(md.mesh.numberofvertices,1); +md.gia.Ngia=zeros(md.mesh.numberofvertices,1); +md.gia.Ugia=zeros(md.mesh.numberofvertices,1); md.gia.modelid=1; -%Uncertainty Quantification%{{{ +%Uncertainty Quantification md.qmu.variables=struct();; -%partition vector: -npart=4; +%partitioning scaling +npart=1; partition=-ones(md.mesh.numberofelements,1); -partition(posgre)=0; -partition(posant)=1; +partition(posala)=0; + +%variable scaling +md.qmu.variables.surfaceload0=normal_uncertain(... + 'descriptor','scaled_SurfaceloadIceThicknessChange',... + 'mean',ones(npart,nsteps),... + 'stddev',.1*ones(npart,nsteps),... + 'partition',partition,'nsteps',nsteps); + +%partitioning distributed +npart=3; +partition=-ones(md.mesh.numberofelements,1); + +partition(posant)=0; +partition(posgre)=1; partition(posala)=2; -partition(poshim)=3; + +vposant=unique(md.mesh.elements(posant,:)); +vposgre=unique(md.mesh.elements(posgre,:)); +vposala=unique(md.mesh.elements(posala,:)); + %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); +nant=10; ngre=20; nala=30; nmax=max([nant,ngre,nala]); + +%make ocean levelset temporal: +md.mask.ocean_levelset=[md.mask.ocean_levelset;2000]; + +%background arrays: +bkgd=md.solidearth.surfaceload.icethicknesschange; +bkgdnan=bkgd; bkgdnan(1:end-1,:)=NaN; + +bkgdocean=md.mask.ocean_levelset; +bkgdoceannan=bkgdocean; bkgdoceannan(1:end-1,1)=NaN; + +%mme arrays: +md.solidearth.surfaceload.icethicknesschange=cell(30,1); +md.mask.ocean_levelset=cell(30,1); +for i=1:nmax, + if i==1, + md.solidearth.surfaceload.icethicknesschange{i}=bkgd; + md.mask.ocean_levelset{i}=bkgdocean; + else + bi=bkgdnan; + oi=bkgdoceannan; + if i<=10, + bi(posant,:)=(-50+i)*ones(length(posant),nsteps); + %oi(vposant,:)=bkgdocean(vposant,:)-50+i; + oi(vposant,:)=i; + end + if i<=20, + bi(posgre,:)=(-100+i)*ones(length(posgre),nsteps); + %oi(vposgre,:)=bkgdocean(vposgre,:)-100+i; + oi(vposgre,:)=2*i; + end + if i<=30, + bi(posala,:)=(-150+i)*ones(length(posala),nsteps); + %oi(vposala,:)=bkgdocean(vposala,:)-150+i; + oi(vposala,:)=3*i; + end + md.solidearth.surfaceload.icethicknesschange{i}=bi; + ois=repmat(oi,1,nsteps); ois(end,:)=bi(end,:); + md.mask.ocean_levelset{i}=ois; + end + md.mask.ocean_levelset{i}=bkgdocean; +end +%md.mask.ice_levelset=md.mask.ocean_levelset; +%md.solidearth.surfaceload.icethicknesschange=bkgd; +%md.mask.ocean_levelset=bkgdocean; + +%create distributed histograms for each partition: +pairs_per_variable=zeros(npart,1); +abscissas=cell(npart,1); +counts=cell(npart,1); +[abscissas{1} counts{1} pairs_per_variable(1)]=equiprobable_histogram_uncertain(10); +[abscissas{2} counts{2} pairs_per_variable(2)]=equiprobable_histogram_uncertain(20); +[abscissas{3} counts{3} pairs_per_variable(3)]=equiprobable_histogram_uncertain(30); %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.variables.surfaceload1=histogram_bin_uncertain(... + 'descriptor','distributed_SurfaceloadModelid',... + 'pairs_per_variable',pairs_per_variable,'abscissas',abscissas,'counts',counts,'partition',partition); + + md.qmu.correlation_matrix=[]; - - %responses @@ -179,9 +211,8 @@ md.qmu.params.tabular_graphics_data=true; md.qmu.output=1; -%}}} %transient: md.timestepping.start_time=2000; -md.timestepping.interp_forcings=0; +md.timestepping.interp_forcings=1; md.timestepping.final_time=2002; md.transient.issmb=0; @@ -191,8 +222,8 @@ md.transient.isslr=1; md.transient.isgia=1; -md.slr.requested_outputs= {'default',... - 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... +md.solidearth.requested_outputs= {'default',... + 'SurfaceloadIceThicknessChange','Sealevel','SealevelRSLRate',... 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate',... - 'SealevelEustaticMask','SealevelEustaticOceanMask'}; + 'SealevelEustaticMask','SealevelEustaticOceanMask','MaskOceanLevelset','MaskIceLevelset'}; %hack: @@ -203,7 +234,8 @@ %eustatic + rigid + elastic + rotation run: md.verbose=verbose('11111111111'); -%md.verbose.qmu=1; +%md.verbose=verbose(0); +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'); - +md=solve(md,'tr'); + Index: /issm/trunk-jpl/test/SandBox/partmme2002.m =================================================================== --- /issm/trunk-jpl/test/SandBox/partmme2002.m (revision 25267) +++ /issm/trunk-jpl/test/SandBox/partmme2002.m (revision 25268) @@ -3,6 +3,6 @@ %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',400.); %500 km resolution mesh -md.cluster.np=16; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %500 km resolution mesh +md.cluster.np=2; %parameterize slr solution: @@ -140,5 +140,9 @@ md.mask.ocean_levelset{i}=ois; end + md.mask.ocean_levelset{i}=bkgdocean; end +%md.mask.ice_levelset=md.mask.ocean_levelset; +%md.solidearth.surfaceload.icethicknesschange=bkgd; +%md.mask.ocean_levelset=bkgdocean; %create distributed histograms for each partition: @@ -190,5 +194,5 @@ md.qmu.params.analysis_driver='matlab'; md.qmu.params.evaluation_scheduling='master'; -md.qmu.params.processors_per_evaluation=5; +md.qmu.params.processors_per_evaluation=1; md.qmu.params.tabular_graphics_data=true; md.qmu.output=1; @@ -216,6 +220,6 @@ %eustatic + rigid + elastic + rotation run: -%md.verbose=verbose('11111111111'); -md.verbose=verbose(0); +md.verbose=verbose('11111111111'); +%md.verbose=verbose(0); md.verbose.qmu=1; md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1;md.solidearth.settings.rotation=1; Index: /issm/trunk-jpl/test/SandBox/scale2002_element.m =================================================================== --- /issm/trunk-jpl/test/SandBox/scale2002_element.m (revision 25268) +++ /issm/trunk-jpl/test/SandBox/scale2002_element.m (revision 25268) @@ -0,0 +1,174 @@ +%Test Name: Dakota EarthSlr Scaling. + +%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: +nsteps=100; +md.solidearth.surfaceload.icethicknesschange=ones(md.mesh.numberofelements+1,nsteps); +md.solidearth.surfaceload.icethicknesschange(end,:)=2000.5:1:2099.5; +md.solidearth.sealevel=zeros(md.mesh.numberofvertices,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>-50 & longe<-30); +%md.solidearth.surfaceload.icethicknesschange(posgre,:)=-100; + +%alaska : +posala=find(late > 55 & late < 68 & longe>-162 & longe<-140); +%md.solidearth.surfaceload.icethicknesschange(posala,:)=-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; + +%dsl +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); + +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='scale2002'; + +%Uncertainty Quantification% +md.qmu.variables=struct();; + +%partitioning scaling +npart=3; +partition=-ones(md.mesh.numberofelements,1); +partition(posgre)=0; +partition(posant)=1; +partition(posala)=2; + +%variable scaling +md.qmu.variables.surfaceload0=normal_uncertain(... + 'descriptor','scaled_SurfaceloadIceThicknessChange',... + 'mean',ones(npart,nsteps),... + 'stddev',.1*ones(npart,nsteps),... + 'partition',partition,'nsteps',nsteps); + +%responses +md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1'); + +%output definitions: +md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ... + 'model_string','Sealevel','node',1)}; + +%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','random'); + +%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.solidearth.settings.computesealevelchange=1; +md.timestepping.start_time=1999.5; +md.timestepping.interp_forcings=0; +md.timestepping.time_step=1; +md.timestepping.final_time=2099.5 +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.solidearth.requested_outputs= {'default',... + 'SurfaceloadIceThicknessChange','Sealevel','SealevelRSLRate',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate',... + 'SealevelEustaticMask','SealevelEustaticOceanMask','MaskOceanLevelset','MaskIceLevelset'}; + +%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); + +%settings: +md.verbose.qmu=1; +md.verbose=verbose(0); +md.qmu.isdakota=1; + +md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1;md.solidearth.settings.rotation=1; + +md=solve(md,'tr'); + +%check scaling worked ok: +n=md.qmu.method.params.samples; +field_values={}; + +for i=1:n, + md2=model(); md2.results=md.results.dakota.modelresults{i}; + uq=md2.results.TransientSolution(1).uq_variables; + part1=uq(1:100); + part2=uq(101:200); + part3=uq(201:300); + + h2=ones(md.mesh.numberofelements+1,100); + + for j=1:length(posgre), + for k=1:100, + h2(posgre(j),k)= h2(posgre(j),k)*part1(k); + end + end + for j=1:length(posant), + for k=1:100, + h2(posant(j),k)= h2(posant(j),k)*part2(k); + end + end + for j=1:length(posala), + for k=1:100, + h2(posala(j),k)= h2(posala(j),k)*part3(k); + end + end + + h=resultstomatrix(md2,'TransientSolution','SurfaceloadIceThicknessChange'); + h2(end,:)=h(end,:); + + field_values{end+1}=size(find(h2-h)); +end + + +%Fields and tolerances to track changes +field_names ={'dH1','dH2','dH3'}; +field_tolerances={1e-5,1e-5,1e-5}; + Index: /issm/trunk-jpl/test/SandBox/scale2002_vertices.m =================================================================== --- /issm/trunk-jpl/test/SandBox/scale2002_vertices.m (revision 25268) +++ /issm/trunk-jpl/test/SandBox/scale2002_vertices.m (revision 25268) @@ -0,0 +1,81 @@ +%Test Name: SquareSheetShelfDiadSSA3dDakota +md=triangle(model(),'../Exp/Square.exp',300000.); +md=setmask(md,'../Exp/SquareShelf.exp',''); +md=parameterize(md,'../Par/SquareSheetShelf.par'); +md=setflowequation(md,'SSA','all'); +md.cluster=generic('name',oshostname(),'np',2); + +%partitioning +partition=partitioner(md,'package','linear','npart',md.mesh.numberofvertices)-1; +md.qmu.isdakota=1; + +nsteps=100; +md.friction.coefficient=ones(md.mesh.numberofvertices+1,nsteps); +md.friction.coefficient(end,:)=2000.5:1:2099.5; + +%variables +md.qmu.variables.drag_coefficient=normal_uncertain('descriptor','scaled_FrictionCoefficient',... + 'mean',ones(md.mesh.numberofvertices,nsteps),... + 'stddev',.01*ones(md.mesh.numberofvertices,nsteps),... + 'partition',partition,'nsteps',nsteps); + +%responses +md.qmu.responses.MaxVel=response_function('descriptor','MaxVel'); + +%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','random'); + +%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; + +md.verbose=verbose(0); + +md.timestepping.start_time=1999.5; +md.timestepping.interp_forcings=0; +md.timestepping.time_step=1; +md.timestepping.final_time=2099.5 +md.transient.issmb=0; +md.transient.ismasstransport=0; +md.transient.isthermal=0; + +md.stressbalance.requested_outputs{end+1}='FrictionCoefficient'; + +%imperative! +md.stressbalance.reltol=10^-5; %tighten for qmu analysese + +%solve +md=solve(md,'tr','overwrite','y'); + +%check scaling worked ok: +n=md.qmu.method.params.samples; +field_values={}; + +for i=1:n, + md2=model(); md2.results=md.results.dakota.modelresults{i}; + uq=md2.results.TransientSolution(1).uq_variables; + + h2=ones(md.mesh.numberofvertices+1,100); + + for j=1:md.mesh.numberofvertices, + partj=uq((j-1)*100+1 : j*100) + h2(j,:)=h2(j,:)'.*partj; + end + + h=resultstomatrix(md2,'TransientSolution','FrictionCoefficient'); + h2(end,:)=h(end,:); + size(find(h2-h)) + + field_values{end+1}=size(find(h2-h)); +end + +