Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26282) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26283) @@ -338,5 +338,5 @@ for(int i=lower_row;i(i)*degacc * PI / 180.0; + alpha= reCast(i)*degacc * M_PI / 180.0; G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); } @@ -345,5 +345,5 @@ for(int i=lower_row;i(i)*degacc * PI / 180.0; + alpha= reCast(i)*degacc * M_PI / 180.0; for(int t=0;t(pice[eid]); + ice->SetValue(id,0.,INS_VAL); + } + else{ + ice->SetValue(0,0.,INS_VAL); + } + + if(nhydro){ + id=reCast(phydro[eid]); + hydro->SetValue(id,0.,INS_VAL); + } + else{ + hydro->SetValue(0,0.,INS_VAL); + } + + if(nocean){ + id=reCast(pocean[eid]); + ocean->SetValue(id,0.,INS_VAL); + } + else{ + ocean->SetValue(0,0.,INS_VAL); + } + + +} /*}}}*/ void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/BarystaticContributions.h =================================================================== --- /issm/trunk-jpl/src/c/classes/BarystaticContributions.h (revision 26282) +++ /issm/trunk-jpl/src/c/classes/BarystaticContributions.h (revision 26283) @@ -42,4 +42,5 @@ void Save(Results* results, Parameters* parameters, IssmDouble oceanarea); void Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue); + void Reset(int eid); }; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26282) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26283) @@ -6218,4 +6218,5 @@ } + constant=3/rho_earth/planetarea; @@ -6245,5 +6246,5 @@ alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); index=reCast( alpha/M_PI*reCast(M-1) ); - _assert_(index>0 && index=0 && indexinputs->SetArrayInput(StackRSLEnum,this->lid,stackRSL,3*nt); this->inputs->SetArrayInput(StackUEnum,this->lid,stackU,3*nt); + this->parameters->SetParam(0,StackIndexEnum); if(horiz){ stackN=xNewZeroInit(3*nt); @@ -6955,6 +6957,6 @@ IssmDouble* stacktimes=NULL; int stacknumsteps; - int stackindex; - int newindex; + int stackindex=0; + int newindex=0; int dummy; bool viscous=false; @@ -6980,25 +6982,33 @@ } + + bool foundtime=false; + int offset=1; lincoeff=0; newindex=stacknumsteps-2; + for(int t=stackindex;tcurrenttime){ newindex=t-1; lincoeff=(currenttime-stacktimes[newindex])/(stacktimes[t]-stacktimes[newindex]); + foundtime=true; + offset=0; break; } } - if(newindex==(stacknumsteps-2))lincoeff=1; + + if(!foundtime) lincoeff=1; stacktimes[newindex]=currenttime; for(int i=0;iparameters->SetParam(stackindex,StackIndexEnum); this->parameters->SetParam(stacktimes,stacknumsteps,StackTimesEnum); @@ -7162,6 +7172,8 @@ int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL); + loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL); } else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); + /*add ocean area into a global oceanareas vector:*/ @@ -7395,8 +7407,9 @@ sealevelinterp=xNew(3*nt); if(stacktimes[stackindex]lids[this->vertices[i]->lid]!=this->lid)continue; sealevelinterp[i*nt+t]= sealevel[i*nt+0]; } @@ -7404,4 +7417,5 @@ else{ for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; sealevelinterp[i*nt+t]= (1-lincoeff)*sealevel[i*nt+(t-stackindex-1)]+lincoeff*sealevel[i*nt+(t-stackindex)]; } @@ -7413,11 +7427,14 @@ /*update sealevel at present time using stack at present time: */ for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; sealevel[i*nt+0]+=stack[i*stacknumsteps+stackindex]; } if(computefuture){ /*update stack with future deformation from present load: */ - for(int t=stackindex;t=stackindex;t--){ for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; + stack[i*stacknumsteps+t]+=sealevelinterp[i*nt+t]-sealevelinterp[i*nt+stackindex]-stack[i*stacknumsteps+stackindex]; } } @@ -7433,5 +7450,7 @@ if(percpu){ for(i=0;ilids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0]; + if(slgeom->lids[this->vertices[i]->lid]==this->lid){ + sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0]; + } } } @@ -7439,4 +7458,5 @@ for(i=0;iCumulate(femmodel->parameters); barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); + + //Avoid double counting barystatic contributions in future time steps: reset all non-cumulative bscl to 0 + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + barycontrib->Reset(element->Sid()); + } } @@ -714,6 +720,5 @@ delete vsubsealevelloadsvolume; - //return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea; - return (sealevelloadsaverage)/oceanarea; + return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea; } /*}}}*/ void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/ Index: /issm/trunk-jpl/test/NightlyRun/test2090.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2090.m (revision 26282) +++ /issm/trunk-jpl/test/NightlyRun/test2090.m (revision 26283) @@ -3,10 +3,10 @@ %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh %Geometry for the bed, arbitrary thickness of 1000: md.geometry.bed=-ones(md.mesh.numberofvertices,1); md.geometry.base=md.geometry.bed; -md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1); +md.geometry.thickness=100*ones(md.mesh.numberofvertices,1); md.geometry.surface=md.geometry.bed+md.geometry.thickness; @@ -15,7 +15,5 @@ %solidearth loading: {{{ md.masstransport.spcthickness=[md.geometry.thickness;0]; -md.dsl.global_average_thermosteric_sea_level=[0;0]; -md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1); -md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); %antarctica @@ -28,9 +26,10 @@ longe=atan2d(ye,xe); pos=find(late < -80); -md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +%pos=10; +md.masstransport.spcthickness(md.mesh.elements(pos,:))=md.masstransport.spcthickness(md.mesh.elements(pos,:)) -100; posant=pos; %greenland pos=find(late>60 & late<90 & longe>-75 & longe<-15); -md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +md.masstransport.spcthickness(md.mesh.elements(pos,:))=md.masstransport.spcthickness(md.mesh.elements(pos,:)) -100; posgre=pos; @@ -39,10 +38,15 @@ load ../Data/lnb_temporal.mat md.solidearth.lovenumbers=lovenumbers('maxdeg',1000); -md.solidearth.lovenumbers.h=[zeros(1,100); h1]; -md.solidearth.lovenumbers.k=[zeros(1,100); k1]; -md.solidearth.lovenumbers.l=[zeros(1,100); l1]; -md.solidearth.lovenumbers.th=repmat(md.solidearth.lovenumbers.th,1,100); -md.solidearth.lovenumbers.tk=repmat(md.solidearth.lovenumbers.tk,1,100); -md.solidearth.lovenumbers.tl=repmat(md.solidearth.lovenumbers.tl,1,100); +md.solidearth.lovenumbers.timefreq=[0]; + +md.solidearth.lovenumbers.h=[h1']; +md.solidearth.lovenumbers.k=[k1']; +md.solidearth.lovenumbers.l=[l1']; +%md.solidearth.lovenumbers.h=repmat(md.solidearth.lovenumbers.h,1,100); +%md.solidearth.lovenumbers.k=repmat(md.solidearth.lovenumbers.k,1,100); +%md.solidearth.lovenumbers.l=repmat(md.solidearth.lovenumbers.l,1,100); +md.solidearth.lovenumbers.th=repmat(md.solidearth.lovenumbers.th,1,101); +md.solidearth.lovenumbers.tk=repmat(md.solidearth.lovenumbers.tk,1,101); +md.solidearth.lovenumbers.tl=repmat(md.solidearth.lovenumbers.tl,1,101); md.solidearth.lovenumbers.timefreq=time; @@ -56,4 +60,5 @@ icemask(md.mesh.elements(posant))=-1; icemask(md.mesh.elements(posgre))=-1; +%oceanmask(md.mesh.elements(posant))=1; md.mask.ice_levelset=icemask; @@ -63,6 +68,11 @@ %time stepping: md.timestepping.start_time=0; -md.timestepping.time_step=1; -md.timestepping.final_time=10; +md.timestepping.time_step=100; +md.timestepping.final_time=1000; + +time1=md.timestepping.start_time:md.timestepping.time_step:md.timestepping.final_time; +md.masstransport.spcthickness=repmat(md.masstransport.spcthickness, [1 length(time1)]); +md.masstransport.spcthickness(1:end-1,2:end)=0; +md.masstransport.spcthickness(end,:)=time1; md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); @@ -82,5 +92,5 @@ %Solution parameters -md.cluster.np=3; +md.cluster.np=1; md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; @@ -89,6 +99,8 @@ md.solidearth.settings.ocean_area_scaling=0; md.solidearth.settings.grdmodel=1; +md.solidearth.settings.timeacc=md.timestepping.time_step; +md.solidearth.settings.degacc=.01; -%Physics: +%Physics: md.transient.issmb=0; md.transient.isstressbalance=0; @@ -96,6 +108,5 @@ md.transient.ismasstransport=1; md.transient.isslc=1; -md.solidearth.requested_outputs={'Sealevel','Bed','BedGRD'}; - +md.solidearth.requested_outputs={'Sealevel','Bed'}; % max number of iteration reverted back to 10 (i.e., the original default value) @@ -105,10 +116,17 @@ md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; +md.solidearth.settings.viscous=1; md.solidearth.settings.rotation=0; md=solve(md,'tr'); -S=md.results.TransientSolution.Sealevel; + + +clear S B +for i=length(time1)-1 +S(:,i)=md.results.TransientSolution(i).Sealevel; +B(:,i)=md.results.TransientSolution(i).Bed; +end %Fields and tolerances to track changes -field_names={'Sealevel'}; -field_tolerances={1e-13}; -field_values={S}; +field_names={'Sealevel','Bed'}; +field_tolerances={1e-13,1e-13}; +field_values={S,B};