Changeset 26283


Ignore:
Timestamp:
05/19/21 16:03:18 (4 years ago)
Author:
caronlam
Message:

CHG: grdcore viscous mode now computes the delta_deformation between two calls rather than the total deformation, in line with elastic mode. BUG: Fixing issues related to viscous stack & resetting barystatic contribs every call in order to shift sealevel fields by the delta_barycontrib instead of cummulated barycontrib.

Location:
issm/trunk-jpl
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp

    r26278 r26283  
    338338                        for(int i=lower_row;i<upper_row;i++){
    339339                                IssmDouble alpha,x;
    340                                 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     340                                alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
    341341                                G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
    342342                        }
    … …  
    345345                        for(int i=lower_row;i<upper_row;i++){
    346346                                IssmDouble alpha,x;
    347                                 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     347                                alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
    348348
    349349                                for(int t=0;t<ntimesteps;t++){
    … …  
    478478                                                int timeindex=index*nt+t;
    479479                                                int timepreindex= index*ntimesteps+timeindex2;
    480                                                 G_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*G_viscoelastic[timepreindex]+lincoeff*G_viscoelastic[timepreindex+1];
    481                                                 U_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*U_viscoelastic[timepreindex]+lincoeff*U_viscoelastic[timepreindex+1];
    482                                                 if(horiz)H_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*H_viscoelastic[timepreindex]+lincoeff*H_viscoelastic[timepreindex+1];
     480                                                G_viscoelastic_interpolated[timeindex]=(1-lincoeff)*G_viscoelastic[timepreindex]+lincoeff*G_viscoelastic[timepreindex+1];
     481                                                U_viscoelastic_interpolated[timeindex]=(1-lincoeff)*U_viscoelastic[timepreindex]+lincoeff*U_viscoelastic[timepreindex+1];
     482                                                if(horiz)H_viscoelastic_interpolated[timeindex]=(1-lincoeff)*H_viscoelastic[timepreindex]+lincoeff*H_viscoelastic[timepreindex+1];
    483483                                        }
    484484                                }

  • issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp

    r26149 r26283  
    131131
    132132} /*}}}*/
     133void BarystaticContributions::Reset(int eid){ /*{{{*/
     134       
     135        int id;
     136        if(nice){
     137                id=reCast<int>(pice[eid]);
     138                ice->SetValue(id,0.,INS_VAL);
     139        }
     140        else{
     141                ice->SetValue(0,0.,INS_VAL);
     142        }
     143
     144        if(nhydro){
     145                id=reCast<int>(phydro[eid]);
     146                hydro->SetValue(id,0.,INS_VAL);
     147        }
     148        else{
     149                hydro->SetValue(0,0.,INS_VAL);
     150        }
     151
     152        if(nocean){
     153                id=reCast<int>(pocean[eid]);
     154                ocean->SetValue(id,0.,INS_VAL);
     155        }
     156        else{
     157                ocean->SetValue(0,0.,INS_VAL);
     158        }
     159               
     160
     161} /*}}}*/
    133162void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/
    134163

  • issm/trunk-jpl/src/c/classes/BarystaticContributions.h

    r26121 r26283  
    4242                void Save(Results* results, Parameters* parameters, IssmDouble oceanarea);
    4343                void Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue);
     44                void Reset(int eid);
    4445
    4546};

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r26276 r26283  
    62186218        }
    62196219
     6220
    62206221        constant=3/rho_earth/planetarea;
    62216222
    … …  
    62456246                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    62466247                        index=reCast<int,IssmDouble>( alpha/M_PI*reCast<IssmDouble,int>(M-1) );
    6247                         _assert_(index>0 && index<M);
     6248                        _assert_(index>=0 && index<M);
    62486249
    62496250                        if(horiz){
    … …  
    63826383                this->inputs->SetArrayInput(StackRSLEnum,this->lid,stackRSL,3*nt);
    63836384                this->inputs->SetArrayInput(StackUEnum,this->lid,stackU,3*nt);
     6385                this->parameters->SetParam(0,StackIndexEnum);
    63846386                if(horiz){
    63856387                        stackN=xNewZeroInit<IssmDouble>(3*nt);
    … …  
    69556957        IssmDouble* stacktimes=NULL;
    69566958        int         stacknumsteps;
    6957         int         stackindex;
    6958         int         newindex;
     6959        int         stackindex=0;
     6960        int         newindex=0;
    69596961        int         dummy;
    69606962        bool        viscous=false;
    … …  
    69806982                }
    69816983
     6984
     6985                bool foundtime=false;
     6986                int offset=1;
    69826987                lincoeff=0;
    69836988                newindex=stacknumsteps-2;
     6989
    69846990                for(int t=stackindex;t<stacknumsteps;t++){
    69856991                        if (stacktimes[t]>currenttime){
    69866992                                newindex=t-1;
    69876993                                lincoeff=(currenttime-stacktimes[newindex])/(stacktimes[t]-stacktimes[newindex]);
     6994                                foundtime=true;
     6995                                offset=0;
    69886996                                break;
    69896997                        }
    69906998                }
    6991                 if(newindex==(stacknumsteps-2))lincoeff=1;
     6999
     7000                if(!foundtime) lincoeff=1;
    69927001                stacktimes[newindex]=currenttime;
    69937002                for(int i=0;i<NUMVERTICES;i++){
    6994                         stackRSL[i*stacknumsteps+newindex]=(1-lincoeff)*stackRSL[i*stacknumsteps+newindex]+lincoeff*stackRSL[i*stacknumsteps+newindex+1];
    6995                         stackU[i*stacknumsteps+newindex]=(1-lincoeff)*stackU[i*stacknumsteps+newindex]+lincoeff*stackU[i*stacknumsteps+newindex+1];
     7003                        stackRSL[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackRSL[i*stacknumsteps+newindex]+lincoeff*stackRSL[i*stacknumsteps+newindex+1];
     7004                        stackU[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackU[i*stacknumsteps+newindex]+lincoeff*stackU[i*stacknumsteps+newindex+1];
    69967005                        if(horiz){
    6997                                 stackN[i*stacknumsteps+newindex]=(1-lincoeff)*stackN[i*stacknumsteps+newindex]+lincoeff*stackN[i*stacknumsteps+newindex+1];
    6998                                 stackE[i*stacknumsteps+newindex]=(1-lincoeff)*stackE[i*stacknumsteps+newindex]+lincoeff*stackE[i*stacknumsteps+newindex+1];
    6999                         }
    7000                 }
    7001                 stackindex=newindex;
    7002 
     7006                                stackN[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackN[i*stacknumsteps+newindex]+lincoeff*stackN[i*stacknumsteps+newindex+1];
     7007                                stackE[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackE[i*stacknumsteps+newindex]+lincoeff*stackE[i*stacknumsteps+newindex+1];
     7008                        }
     7009                }
     7010                stackindex=newindex+offset;
     7011
     7012               
    70037013                this->parameters->SetParam(stackindex,StackIndexEnum);
    70047014                this->parameters->SetParam(stacktimes,stacknumsteps,StackTimesEnum);
    … …  
    71627172                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    71637173                loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL);
     7174                loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
    71647175        }
    71657176        else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
     7177
    71667178
    71677179        /*add ocean area into a global oceanareas vector:*/
    … …  
    73957407                        sealevelinterp=xNew<IssmDouble>(3*nt);
    73967408                        if(stacktimes[stackindex]<final_time){
    7397                                 lincoeff=1-(stacktimes[stackindex+1]-stacktimes[stackindex])/timeacc;
     7409                                lincoeff=(stacktimes[stackindex+1]-stacktimes[stackindex])/timeacc;
    73987410                                for(int t=stackindex;t<nt;t++){
    73997411                                        if(t==stackindex){
    74007412                                                for(int i=0;i<NUMVERTICES;i++){
     7413                                                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    74017414                                                        sealevelinterp[i*nt+t]=  sealevel[i*nt+0];
    74027415                                                }
    … …  
    74047417                                        else{
    74057418                                                for(int i=0;i<NUMVERTICES;i++){
     7419                                                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    74067420                                                        sealevelinterp[i*nt+t]=  (1-lincoeff)*sealevel[i*nt+(t-stackindex-1)]+lincoeff*sealevel[i*nt+(t-stackindex)];
    74077421                                                }
    … …  
    74137427                /*update sealevel at present time using stack at present time: */
    74147428                for(int i=0;i<NUMVERTICES;i++){
     7429                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    74157430                        sealevel[i*nt+0]+=stack[i*stacknumsteps+stackindex];
    74167431                }
    74177432
    74187433                if(computefuture){ /*update stack with future deformation from present load: */
    7419                         for(int t=stackindex;t<nt;t++){
     7434
     7435                        for(int t=nt-1;t>=stackindex;t--){
    74207436                                for(int i=0;i<NUMVERTICES;i++){
    7421                                         stack[i*stacknumsteps+t]+=sealevelinterp[i*nt+t];
     7437                                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7438                                        stack[i*stacknumsteps+t]+=sealevelinterp[i*nt+t]-sealevelinterp[i*nt+stackindex]-stack[i*stacknumsteps+stackindex];
    74227439                                }
    74237440                        }
    … …  
    74337450        if(percpu){
    74347451                for(i=0;i<NUMVERTICES;i++){
    7435                         if(slgeom->lids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0];
     7452                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     7453                                sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0];
     7454                        }
    74367455                }
    74377456        }
    … …  
    74397458                for(i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i*nt+0];
    74407459        }
     7460
    74417461
    74427462} /*}}}*/

  • issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

    r26272 r26283  
    388388                barycontrib->Cumulate(femmodel->parameters);
    389389                barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
     390
     391                //Avoid double counting barystatic contributions in future time steps: reset all non-cumulative bscl to 0
     392                for(Object* & object : femmodel->elements->objects){
     393                        Element* element = xDynamicCast<Element*>(object);
     394                        barycontrib->Reset(element->Sid());
     395                }
    390396        }
    391397
    … …  
    714720        delete vsubsealevelloadsvolume;
    715721       
    716         //return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea;
    717         return (sealevelloadsaverage)/oceanarea;
     722        return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea;
    718723} /*}}}*/
    719724void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/

  • issm/trunk-jpl/test/NightlyRun/test2090.m

    r26272 r26283  
    33%mesh earth:
    44md=model;
    5 md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh
     5md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh
    66
    77%Geometry for the bed, arbitrary thickness of 1000:
    88md.geometry.bed=-ones(md.mesh.numberofvertices,1);
    99md.geometry.base=md.geometry.bed;
    10 md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
     10md.geometry.thickness=100*ones(md.mesh.numberofvertices,1);
    1111md.geometry.surface=md.geometry.bed+md.geometry.thickness;
    1212
    … …  
    1515%solidearth loading:  {{{
    1616md.masstransport.spcthickness=[md.geometry.thickness;0];
    17 md.dsl.global_average_thermosteric_sea_level=[0;0];
    18 md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1);
    19 md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
     17
    2018md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
    2119%antarctica
    … …  
    2826longe=atan2d(ye,xe);
    2927pos=find(late < -80);
    30 md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
     28%pos=10;
     29md.masstransport.spcthickness(md.mesh.elements(pos,:))=md.masstransport.spcthickness(md.mesh.elements(pos,:)) -100;
    3130posant=pos;
    3231%greenland
    3332pos=find(late>60 & late<90 & longe>-75 & longe<-15);
    34 md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
     33md.masstransport.spcthickness(md.mesh.elements(pos,:))=md.masstransport.spcthickness(md.mesh.elements(pos,:)) -100;
    3534posgre=pos;
    3635
    … …  
    3938load ../Data/lnb_temporal.mat
    4039md.solidearth.lovenumbers=lovenumbers('maxdeg',1000);
    41 md.solidearth.lovenumbers.h=[zeros(1,100); h1];
    42 md.solidearth.lovenumbers.k=[zeros(1,100); k1];
    43 md.solidearth.lovenumbers.l=[zeros(1,100); l1];
    44 md.solidearth.lovenumbers.th=repmat(md.solidearth.lovenumbers.th,1,100);
    45 md.solidearth.lovenumbers.tk=repmat(md.solidearth.lovenumbers.tk,1,100);
    46 md.solidearth.lovenumbers.tl=repmat(md.solidearth.lovenumbers.tl,1,100);
     40md.solidearth.lovenumbers.timefreq=[0];
     41
     42md.solidearth.lovenumbers.h=[h1'];
     43md.solidearth.lovenumbers.k=[k1'];
     44md.solidearth.lovenumbers.l=[l1'];
     45%md.solidearth.lovenumbers.h=repmat(md.solidearth.lovenumbers.h,1,100);
     46%md.solidearth.lovenumbers.k=repmat(md.solidearth.lovenumbers.k,1,100);
     47%md.solidearth.lovenumbers.l=repmat(md.solidearth.lovenumbers.l,1,100);
     48md.solidearth.lovenumbers.th=repmat(md.solidearth.lovenumbers.th,1,101);
     49md.solidearth.lovenumbers.tk=repmat(md.solidearth.lovenumbers.tk,1,101);
     50md.solidearth.lovenumbers.tl=repmat(md.solidearth.lovenumbers.tl,1,101);
    4751md.solidearth.lovenumbers.timefreq=time;
    4852
    … …  
    5660icemask(md.mesh.elements(posant))=-1;
    5761icemask(md.mesh.elements(posgre))=-1;
     62%oceanmask(md.mesh.elements(posant))=1;
    5863
    5964md.mask.ice_levelset=icemask;
    … …  
    6368%time stepping:
    6469md.timestepping.start_time=0;
    65 md.timestepping.time_step=1;
    66 md.timestepping.final_time=10;
     70md.timestepping.time_step=100;
     71md.timestepping.final_time=1000;
     72
     73time1=md.timestepping.start_time:md.timestepping.time_step:md.timestepping.final_time;
     74md.masstransport.spcthickness=repmat(md.masstransport.spcthickness, [1 length(time1)]);
     75md.masstransport.spcthickness(1:end-1,2:end)=0;
     76md.masstransport.spcthickness(end,:)=time1;
    6777
    6878md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
    … …  
    8292
    8393%Solution parameters
    84 md.cluster.np=3;
     94md.cluster.np=1;
    8595md.solidearth.settings.reltol=NaN;
    8696md.solidearth.settings.abstol=1e-3;
    … …  
    8999md.solidearth.settings.ocean_area_scaling=0;
    90100md.solidearth.settings.grdmodel=1;
     101md.solidearth.settings.timeacc=md.timestepping.time_step;
     102md.solidearth.settings.degacc=.01;
    91103
    92 %Physics: 
     104%Physics:
    93105md.transient.issmb=0;
    94106md.transient.isstressbalance=0;
    … …  
    96108md.transient.ismasstransport=1;
    97109md.transient.isslc=1;
    98 md.solidearth.requested_outputs={'Sealevel','Bed','BedGRD'};
    99 
     110md.solidearth.requested_outputs={'Sealevel','Bed'};
    100111
    101112% max number of iteration reverted back to 10 (i.e., the original default value)
    … …  
    105116md.solidearth.settings.rigid=1;
    106117md.solidearth.settings.elastic=1;
     118md.solidearth.settings.viscous=1;
    107119md.solidearth.settings.rotation=0;
    108120md=solve(md,'tr');
    109 S=md.results.TransientSolution.Sealevel;
     121
     122
     123clear S B
     124for i=length(time1)-1
     125S(:,i)=md.results.TransientSolution(i).Sealevel;
     126B(:,i)=md.results.TransientSolution(i).Bed;
     127end
    110128
    111129%Fields and tolerances to track changes
    112 field_names={'Sealevel'};
    113 field_tolerances={1e-13};
    114 field_values={S};
     130field_names={'Sealevel','Bed'};
     131field_tolerances={1e-13,1e-13};
     132field_values={S,B};
Note: See TracChangeset for help on using the changeset viewer.