Changeset 26272


Ignore:
Timestamp:
05/14/21 17:27:17 (4 years ago)
Author:
Eric.Larour
Message:

CHG: visco elastic solid-Earth/sea-level solver now activated.

Location:
issm/trunk-jpl
Files:
1 added
15 edited

Legend:

Unmodified
Added
Removed

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

    r26242 r26272  
    8080
    8181        int         nl;
     82        int         ntimesteps;
    8283        IssmDouble* love_h=NULL;
    8384        IssmDouble* love_k=NULL;
    … …  
    8687        IssmDouble* love_tk=NULL;
    8788        IssmDouble* love_tl=NULL;
     89        IssmDouble* love_timefreq=NULL;
     90        bool        love_istime=true;
    8891        int         externalnature=0;
    8992        int         isexternal=0;
    … …  
    9194        IssmDouble* G_rigid = NULL;
    9295        IssmDouble* G_rigid_local = NULL;
    93         IssmDouble* G_elastic = NULL;
    94         IssmDouble* G_elastic_local = NULL;
    95         IssmDouble* U_elastic = NULL;
    96         IssmDouble* U_elastic_local = NULL;
    97         IssmDouble* H_elastic = NULL;
    98         IssmDouble* H_elastic_local = NULL;
     96        IssmDouble* G_viscoelastic = NULL;
     97        IssmDouble* G_viscoelastic_interpolated = NULL;
     98        IssmDouble* G_viscoelastic_local = NULL;
     99        IssmDouble* U_viscoelastic = NULL;
     100        IssmDouble* U_viscoelastic_interpolated = NULL;
     101        IssmDouble* U_viscoelastic_local = NULL;
     102        IssmDouble* H_viscoelastic = NULL;
     103        IssmDouble* H_viscoelastic_interpolated= NULL;
     104        IssmDouble* H_viscoelastic_local = NULL;
    99105        int         M,m,lower_row,upper_row;
    100106        IssmDouble  degacc=.01;
     107        IssmDouble  timeacc;
    101108        IssmDouble  planetradius=0;
    102109        IssmDouble  planetarea=0;
    103110        bool            rigid=false;
    104111        bool            elastic=false;
     112        bool            viscous=false;
    105113        bool            rotation=false;
     114        int         ndeg;
     115        int         horiz;
     116
     117        bool istime=true;
     118        IssmDouble start_time,final_time;
     119        int  nt,precomputednt;
    106120
    107121        int     numoutputs;
    … …  
    120134        IssmDouble*  bslcocean_partition=NULL;
    121135        int          npartice,nparthydro,npartocean,nel;
    122 
     136        int          grdmodel;
    123137
    124138        /*some constant parameters: */
    … …  
    134148        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
    135149        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));
     150        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.viscous",SolidearthSettingsViscousEnum));
    136151        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum));
    137152        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
    … …  
    211226        }
    212227
    213         /*Deal with elasticity {{{*/
    214         iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
    215         iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
    216         iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
    217 
    218         if(elastic | rigid){
    219                 /*compute green functions for a range of angles*/
    220                 iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
    221                 M=reCast<int,IssmDouble>(180./degacc+1.);
    222         }
    223 
    224         /*love numbers: */
    225         if(elastic){
    226                 iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
    227                 iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
    228                 iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
    229                 iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
    230                 iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
    231                 iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
    232 
    233                 parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
    234                 parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
    235                 parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
    236                 parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
    237                 parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
    238                 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
    239 
    240                 // AD performance is sensitive to calls to ensurecontiguous.
    241                 // // Providing "t" will cause ensurecontiguous to be called.
    242                 #ifdef _HAVE_AD_
    243                 G_elastic=xNew<IssmDouble>(M,"t");
    244                 U_elastic=xNew<IssmDouble>(M,"t");
    245                 H_elastic=xNew<IssmDouble>(M,"t");
    246                 #else
    247                 G_elastic=xNew<IssmDouble>(M);
    248                 U_elastic=xNew<IssmDouble>(M);
    249                 H_elastic=xNew<IssmDouble>(M);
    250                 #endif
    251         }
    252         if(rigid){
    253                 #ifdef _HAVE_AD_
    254                 G_rigid=xNew<IssmDouble>(M,"t");
    255                 #else
    256                 G_rigid=xNew<IssmDouble>(M);
    257                 #endif
    258         }
    259        
    260         if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
    261 
    262         if(rigid | elastic){
    263 
    264                 /*compute combined legendre + love number (elastic green function:*/
    265                 m=DetermineLocalSize(M,IssmComm::GetComm());
    266                 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    267         }
    268         if(elastic){
    269                 #ifdef _HAVE_AD_
    270                 G_elastic_local=xNew<IssmDouble>(m,"t");
    271                 U_elastic_local=xNew<IssmDouble>(m,"t");
    272                 H_elastic_local=xNew<IssmDouble>(m,"t");
    273                 #else
    274                 G_elastic_local=xNew<IssmDouble>(m);
    275                 U_elastic_local=xNew<IssmDouble>(m);
    276                 H_elastic_local=xNew<IssmDouble>(m);
    277                 #endif
    278         }
    279         if(rigid){
    280                 #ifdef _HAVE_AD_
    281                 G_rigid_local=xNew<IssmDouble>(m,"t");
    282                 #else
    283                 G_rigid_local=xNew<IssmDouble>(m);
    284                 #endif
    285         }
    286 
    287         if(rigid){
    288                 for(int i=lower_row;i<upper_row;i++){
    289                         IssmDouble alpha,x;
    290                         alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
    291                         G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
    292                 }
    293         }
    294         if(elastic){
    295                 for(int i=lower_row;i<upper_row;i++){
    296                         IssmDouble alpha,x;
    297                         alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
    298 
    299                         G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
    300                         U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
    301                         H_elastic_local[i-lower_row]= 0;
    302                         IssmDouble Pn = 0.;
    303                         IssmDouble Pn1 = 0.;
    304                         IssmDouble Pn2 = 0.;
    305                         IssmDouble Pn_p = 0.;
    306                         IssmDouble Pn_p1 = 0.;
    307                         IssmDouble Pn_p2 = 0.;
    308 
    309                         for (int n=0;n<nl;n++) {
    310                                 IssmDouble deltalove_G;
    311                                 IssmDouble deltalove_U;
    312 
    313                                 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
    314                                 deltalove_U = (love_h[n]-love_h[nl-1]);
    315 
    316                                 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
    317                                 if(n==0){
    318                                         Pn=1;
    319                                         Pn_p=0;
     228        parameters->FindParam(&grdmodel,GrdModelEnum);
     229        if(grdmodel!=IvinsEnum){
     230                /*Deal with elasticity {{{*/
     231                iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
     232                iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
     233                iomodel->FetchData(&viscous,"md.solidearth.settings.viscous");
     234                iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
     235                iomodel->FetchData(&horiz,"md.solidearth.settings.horiz");
     236
     237                if(rigid){
     238                        /*compute green functions for a range of angles*/
     239                        iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
     240                        M=reCast<int,IssmDouble>(180./degacc+1.);
     241                }
     242
     243                /*love numbers: */
     244                if(viscous || elastic){
     245                        int dummy;
     246
     247                        iomodel->FetchData(&timeacc,"md.solidearth.settings.timeacc");
     248                        iomodel->FetchData(&start_time,"md.timestepping.start_time");
     249                        iomodel->FetchData(&final_time,"md.timestepping.final_time");
     250                        iomodel->FetchData(&love_istime,"md.solidearth.lovenumbers.istime");
     251                        iomodel->FetchData(&love_timefreq,&precomputednt,&dummy,"md.solidearth.lovenumbers.timefreq");
     252                        iomodel->FetchData(&love_h,&ndeg,&precomputednt,"md.solidearth.lovenumbers.h");
     253                        iomodel->FetchData(&love_k,&ndeg,&precomputednt,"md.solidearth.lovenumbers.k");
     254                        iomodel->FetchData(&love_l,&ndeg,&precomputednt,"md.solidearth.lovenumbers.l");
     255                        iomodel->FetchData(&love_th,&ndeg,&precomputednt,"md.solidearth.lovenumbers.th");
     256                        iomodel->FetchData(&love_tk,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tk");
     257                        iomodel->FetchData(&love_tl,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tl");
     258
     259                        parameters->AddObject(new DoubleParam(SolidearthSettingsTimeAccEnum,timeacc));
     260                        parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,ndeg,precomputednt));
     261                        parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,ndeg,precomputednt));
     262                        parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,ndeg,precomputednt));
     263                        parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,ndeg,precomputednt));
     264                        parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,ndeg,precomputednt));
     265                        parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,ndeg,precomputednt));
     266                        parameters->AddObject(new DoubleMatParam(LoveTimeFreqEnum,love_timefreq,precomputednt,1));
     267                        parameters->AddObject(new BoolParam(LoveIsTimeEnum,love_istime));
     268
     269                        // AD performance is sensitive to calls to ensurecontiguous.
     270                        // // Providing "t" will cause ensurecontiguous to be called.
     271                        if(viscous){
     272                                IssmDouble* stacktimes=NULL;
     273                                ntimesteps=precomputednt;
     274                                nt=reCast<int>((final_time-start_time)/timeacc)+1;
     275
     276                                parameters->AddObject(new IntParam(StackNumStepsEnum,nt));
     277                                /*Initialize stack times:*/
     278                                stacktimes=xNew<IssmDouble>(nt);
     279                                for(int t=0;t<nt;t++){
     280                                        stacktimes[t]=start_time+timeacc*t;
    320281                                }
    321                                 else if(n==1){
    322                                         Pn = cos(alpha);
    323                                         Pn_p = 1;
     282                                parameters->AddObject(new DoubleVecParam(StackTimesEnum,stacktimes,nt));
     283                                parameters->AddObject(new IntParam(StackIndexEnum,0));
     284                                xDelete<IssmDouble>(stacktimes);
     285                        }
     286                        else {
     287                                ntimesteps=1;
     288                                nt=1;
     289                        }
     290
     291#ifdef _HAVE_AD_
     292                        G_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
     293                        U_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
     294                        if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
     295#else
     296                        G_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
     297                        U_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
     298                        if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
     299#endif
     300                }
     301                if(rigid){
     302#ifdef _HAVE_AD_
     303                        G_rigid=xNew<IssmDouble>(M,"t");
     304#else
     305                        G_rigid=xNew<IssmDouble>(M);
     306#endif
     307                }
     308
     309                if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
     310
     311                if(rigid){
     312
     313                        /*compute combined legendre + love number (elastic green function:*/
     314                        m=DetermineLocalSize(M,IssmComm::GetComm());
     315                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
     316                }
     317                if(viscous | elastic){
     318#ifdef _HAVE_AD_
     319                        G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
     320                        U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
     321                        if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
     322#else
     323                        G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
     324                        U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
     325                        if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
     326#endif
     327                }
     328                if(rigid){
     329#ifdef _HAVE_AD_
     330                        G_rigid_local=xNew<IssmDouble>(m,"t");
     331#else
     332                        G_rigid_local=xNew<IssmDouble>(m);
     333#endif
     334                }
     335
     336                if(rigid){
     337                        for(int i=lower_row;i<upper_row;i++){
     338                                IssmDouble alpha,x;
     339                                alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     340                                G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
     341                        }
     342                }
     343                if(viscous | elastic){
     344                        for(int i=lower_row;i<upper_row;i++){
     345                                IssmDouble alpha,x;
     346                                alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     347
     348                                for(int t=0;t<ntimesteps;t++){
     349                                        G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_k[(ndeg-1)*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t])*G_rigid_local[i-lower_row];
     350                                        U_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_h[(ndeg-1)*precomputednt+t])*G_rigid_local[i-lower_row];
     351                                        if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t]= 0;
    324352                                }
    325                                 else{
    326                                         Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
    327                                         Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
     353
     354                                IssmDouble Pn = 0.;
     355                                IssmDouble Pn1 = 0.;
     356                                IssmDouble Pn2 = 0.;
     357                                IssmDouble Pn_p = 0.;
     358                                IssmDouble Pn_p1 = 0.;
     359                                IssmDouble Pn_p2 = 0.;
     360
     361                                for (int n=0;n<ndeg;n++) {
     362
     363                                        /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
     364                                        if(n==0){
     365                                                Pn=1;
     366                                                Pn_p=0;
     367                                        }
     368                                        else if(n==1){
     369                                                Pn = cos(alpha);
     370                                                Pn_p = 1;
     371                                        }
     372                                        else{
     373                                                Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
     374                                                Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
     375                                        }
     376                                        Pn2=Pn1; Pn1=Pn;
     377                                        Pn_p2=Pn_p1; Pn_p1=Pn_p;
     378
     379                                        for(int t=0;t<ntimesteps;t++){
     380                                                IssmDouble deltalove_G;
     381                                                IssmDouble deltalove_U;
     382
     383                                                deltalove_G = (love_k[n*precomputednt+t]-love_k[(ndeg-1)*precomputednt+t]-love_h[n*precomputednt+t]+love_h[(ndeg-1)*precomputednt+t]);
     384                                                deltalove_U = (love_h[n*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t]);
     385
     386                                                G_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_G*Pn;                             // gravitational potential
     387                                                U_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_U*Pn;                             // vertical (up) displacement
     388                                                if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t] += sin(alpha)*love_l[n*precomputednt+t]*Pn_p;         // horizontal displacements
     389                                        }
    328390                                }
    329                                 Pn2=Pn1; Pn1=Pn;
    330                                 Pn_p2=Pn_p1; Pn_p1=Pn_p;
    331 
    332                                 G_elastic_local[i-lower_row] += deltalove_G*Pn;         // gravitational potential
    333                                 U_elastic_local[i-lower_row] += deltalove_U*Pn;         // vertical (up) displacement
    334                                 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;              // horizontal displacements
    335                         }
    336                 }
    337         }
    338         if(rigid){
    339 
    340                 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
    341                 int* recvcounts=xNew<int>(IssmComm::GetSize());
    342                 int* displs=xNew<int>(IssmComm::GetSize());
    343 
    344                 //recvcounts:
    345                 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
    346 
    347                 /*displs: */
    348                 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
    349 
    350                 /*All gather:*/
    351                 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    352                 if(elastic){
    353                         ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    354                         ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    355                         ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    356                 }
    357                
    358                 /*free resources: */
    359                 xDelete<int>(recvcounts);
    360                 xDelete<int>(displs);
    361 
    362                 /*Avoid singularity at 0: */
    363                 G_rigid[0]=G_rigid[1];
    364                 if(elastic){
    365                         G_elastic[0]=G_elastic[1];
    366                         U_elastic[0]=U_elastic[1];
    367                         H_elastic[0]=H_elastic[1];
    368                 }
    369                
    370 
    371                 /*Save our precomputed tables into parameters*/
    372                 parameters->AddObject(new DoubleVecParam(SealevelchangeGRigidEnum,G_rigid,M));
    373                 if(elastic){
    374                         parameters->AddObject(new DoubleVecParam(SealevelchangeGElasticEnum,G_elastic,M));
    375                         parameters->AddObject(new DoubleVecParam(SealevelchangeUElasticEnum,U_elastic,M));
    376                         parameters->AddObject(new DoubleVecParam(SealevelchangeHElasticEnum,H_elastic,M));
    377                 }
    378 
    379                 /*free resources: */
    380                 xDelete<IssmDouble>(G_rigid);
    381                 xDelete<IssmDouble>(G_rigid_local);
    382                 if(elastic){
    383                         xDelete<IssmDouble>(love_h);
    384                         xDelete<IssmDouble>(love_k);
    385                         xDelete<IssmDouble>(love_l);
    386                         xDelete<IssmDouble>(love_th);
    387                         xDelete<IssmDouble>(love_tk);
    388                         xDelete<IssmDouble>(love_tl);
    389                         xDelete<IssmDouble>(G_elastic);
    390                         xDelete<IssmDouble>(G_elastic_local);
    391                         xDelete<IssmDouble>(U_elastic);
    392                         xDelete<IssmDouble>(U_elastic_local);
    393                         xDelete<IssmDouble>(H_elastic);
    394                         xDelete<IssmDouble>(H_elastic_local);
    395                 }
    396         } /*}}}*/
    397 
    398         /*Indicate we have not yet run the Geometry Core module: */
    399         parameters->AddObject(new BoolParam(SealevelchangeGeometryDoneEnum,false));
    400         /*}}}*/
     391                        }
     392                }
     393                if(rigid){
     394
     395                        /*merge G_viscoelastic_local into G_viscoelastic; U_viscoelastic_local into U_viscoelastic; H_viscoelastic_local to H_viscoelastic:{{{*/
     396                        int* recvcounts=xNew<int>(IssmComm::GetSize());
     397                        int* displs=xNew<int>(IssmComm::GetSize());
     398                        int  rc;
     399                        int  offset;
     400
     401                        //deal with rigid first:
     402                        ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     403
     404                        /*displs: */
     405                        ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     406
     407                        /*All gather:*/
     408                        ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     409
     410                        if(elastic){
     411                                rc=m*ntimesteps;
     412                                offset=lower_row*ntimesteps;
     413                                ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     414                                ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     415
     416                                ISSM_MPI_Allgatherv(G_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, G_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     417                                ISSM_MPI_Allgatherv(U_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, U_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     418                                if(horiz)ISSM_MPI_Allgatherv(H_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, H_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     419                        }
     420
     421                        /*free resources: */
     422                        xDelete<int>(recvcounts);
     423                        xDelete<int>(displs);
     424
     425                       
     426
     427                        /*Avoid singularity at 0: */
     428                        G_rigid[0]=G_rigid[1];
     429                        if(elastic){
     430                                for(int t=0;t<ntimesteps;t++){
     431                                        G_viscoelastic[t]=G_viscoelastic[ntimesteps+t];
     432                                        U_viscoelastic[t]=U_viscoelastic[ntimesteps+t];
     433                                        if(horiz)H_viscoelastic[t]=H_viscoelastic[ntimesteps+t];
     434                                }
     435                        }
     436
     437                       
     438
     439                        /*Reinterpolate viscoelastic green kernels onto a regular gridded time
     440                         *with steps equal to timeacc:*/
     441                        if(viscous){
     442                                nt=reCast<int>((final_time-start_time)/timeacc)+1;
     443#ifdef _HAVE_AD_
     444                                G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
     445                                U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
     446                                if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
     447#else
     448                                G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
     449                                U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
     450                                if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
     451#endif
     452
     453                                for(int t=0;t<nt;t++){
     454                                        IssmDouble lincoeff;
     455                                        IssmDouble viscoelastic_time=t*timeacc;
     456                                        int        timeindex2=-1;
     457
     458                                        /*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
     459                                        if(t!=0){
     460                                                for(int t2=0;t2<ntimesteps;t2++){
     461                                                        if (viscoelastic_time<love_timefreq[t2]){
     462                                                                timeindex2=t2-1;
     463                                                                if(timeindex2<0)_error_("Temporal Love numbers are computed  with a time accuracy superior to the requested solution time step!");
     464                                                                lincoeff=(viscoelastic_time-love_timefreq[t2-1])/(love_timefreq[t2]-love_timefreq[t2-1]);
     465                                                                break;
     466                                                        }
     467                                                }
     468                                                if(timeindex2==-1)_error_("Temporal love numbers should be extended in time to encompass the requested solution time interval!");
     469                                        }
     470                                        else{
     471                                                timeindex2=0;
     472                                                lincoeff=0;
     473                                        }
     474
     475                                        for(int index=0;index<M;index++){
     476
     477                                                int timeindex=index*nt+t;
     478                                                int timepreindex= index*ntimesteps+timeindex2;
     479                                                G_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*G_viscoelastic[timepreindex]+lincoeff*G_viscoelastic[timepreindex+1];
     480                                                U_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*U_viscoelastic[timepreindex]+lincoeff*U_viscoelastic[timepreindex+1];
     481                                                if(horiz)H_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*H_viscoelastic[timepreindex]+lincoeff*H_viscoelastic[timepreindex+1];
     482                                        }
     483                                }
     484                                               
     485                        }
     486                        else if(elastic){
     487                                nt=1; //in elastic, or if we run only rigid, we need only one step
     488#ifdef _HAVE_AD_
     489                                G_viscoelastic_interpolated=G_viscoelastic;
     490                                U_viscoelastic_interpolated=U_viscoelastic;
     491                                if(horiz)H_viscoelastic_interpolated=H_viscoelastic;
     492#else
     493                                G_viscoelastic_interpolated=G_viscoelastic;
     494                                U_viscoelastic_interpolated=U_viscoelastic;
     495                                if(horiz)H_viscoelastic_interpolated=H_viscoelastic;
     496#endif
     497                        }       
     498
     499                        /*Save our precomputed tables into parameters*/
     500                        parameters->AddObject(new DoubleVecParam(SealevelchangeGRigidEnum,G_rigid,M));
     501                        if(viscous || elastic){
     502                                parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt));
     503                                parameters->AddObject(new DoubleVecParam(SealevelchangeUViscoElasticEnum,U_viscoelastic_interpolated,M*nt));
     504                                if(horiz)parameters->AddObject(new DoubleVecParam(SealevelchangeHViscoElasticEnum,H_viscoelastic_interpolated,M*nt));
     505                        }
     506
     507                        /*free resources: */
     508                        xDelete<IssmDouble>(G_rigid);
     509                        xDelete<IssmDouble>(G_rigid_local);
     510                        if(elastic){
     511                                xDelete<IssmDouble>(love_h);
     512                                xDelete<IssmDouble>(love_k);
     513                                xDelete<IssmDouble>(love_l);
     514                                xDelete<IssmDouble>(love_th);
     515                                xDelete<IssmDouble>(love_tk);
     516                                xDelete<IssmDouble>(love_tl);
     517                                xDelete<IssmDouble>(G_viscoelastic);
     518                                xDelete<IssmDouble>(G_viscoelastic_local);
     519                                xDelete<IssmDouble>(U_viscoelastic);
     520                                xDelete<IssmDouble>(U_viscoelastic_local);
     521                                if(horiz){
     522                                        xDelete<IssmDouble>(H_viscoelastic);
     523                                        xDelete<IssmDouble>(H_viscoelastic_local);
     524                                }
     525                        }
     526                } /*}}}*/
     527
     528                /*Indicate we have not yet run the Geometry Core module: */
     529                parameters->AddObject(new BoolParam(SealevelchangeGeometryDoneEnum,false));
     530                /*}}}*/
     531        }
    401532
    402533        /*Transitions:{{{ */

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r26230 r26272  
    401401                virtual void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
    402402                virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
     403                virtual void       SealevelchangeUpdateViscousFields()=0;
    403404                #endif
    404405

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r26230 r26272  
    234234                void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    235235                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     236                void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
    236237                #endif
    237238

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r26230 r26272  
    189189                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    190190                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     191                void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
    191192#endif
    192193

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r26230 r26272  
    196196                void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    197197                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     198                void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
    198199#endif
    199200

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

    r26242 r26272  
    61006100        IssmDouble xyz_list[NUMVERTICES][3];
    61016101
     6102        /*stacks:*/
     6103        IssmDouble* stackRSL = NULL;
     6104        IssmDouble* stackU = NULL;
     6105        IssmDouble* stackN = NULL;
     6106        IssmDouble* stackE = NULL;
     6107
    61026108        #ifdef _HAVE_RESTRICT_
    61036109        IssmDouble* __restrict__ G=NULL;
    … …  
    61056111        IssmDouble* __restrict__ GN=NULL;
    61066112        IssmDouble* __restrict__ GE=NULL;
    6107         IssmDouble* __restrict__ G_elastic_precomputed=NULL;
     6113        IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
    61086114        IssmDouble* __restrict__ G_rigid_precomputed=NULL;
    6109         IssmDouble* __restrict__ U_elastic_precomputed=NULL;
    6110         IssmDouble* __restrict__ H_elastic_precomputed=NULL;
     6115        IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
     6116        IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
    61116117        #else
    61126118        IssmDouble* G=NULL;
    … …  
    61146120        IssmDouble* GN=NULL;
    61156121        IssmDouble* GE=NULL;
    6116         IssmDouble* G_elastic_precomputed=NULL;
     6122        IssmDouble* G_viscoelastic_precomputed=NULL;
    61176123        IssmDouble* G_rigid_precomputed=NULL;
    6118         IssmDouble* U_elastic_precomputed=NULL;
    6119         IssmDouble* H_elastic_precomputed=NULL;
     6124        IssmDouble* U_viscoelastic_precomputed=NULL;
     6125        IssmDouble* H_viscoelastic_precomputed=NULL;
    61206126        #endif
    61216127       
    6122         /*elastic green function:*/
     6128        /*viscoelastic green function:*/
    61236129        int index;
    61246130        int         M;
    … …  
    61286134        bool computeelastic = false;
    61296135        bool computerotation = false;
     6136        bool computeviscous = false;
    61306137        int  horiz;
     6138        bool istime=true;
     6139        IssmDouble timeacc=0;
     6140        IssmDouble start_time,final_time;
     6141        int  nt,precomputednt;
    61316142
    61326143        /*Rotational:*/
    … …  
    61496160        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    61506161        this->parameters->FindParam(&computerotation,SolidearthSettingsRotationEnum);
     6162        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    61516163        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    61526164        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    … …  
    61736185
    61746186        if(computeelastic){
    6175                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter);
    6176                 parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
    6177 
    6178                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter);
    6179                 parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
    6180 
    6181                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter);
    6182                 parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
    6183 
     6187                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
     6188                parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
    61846189               
     6190                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
     6191                parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
     6192
     6193                if(horiz){
     6194                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
     6195                        parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
     6196                }
    61856197        }
    61866198        /*}}}*/
    … …  
    61946206        /*}}}*/
    61956207        /*Compute green functions:{{{ */
     6208        if(computeviscous){
     6209                this->parameters->FindParam(&istime,LoveIsTimeEnum);
     6210                if(!istime)_error_("Frequency love numbers not supported yet!");
     6211                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     6212                this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
     6213                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     6214                nt=reCast<int>((final_time-start_time)/timeacc)+1;
     6215        }
     6216        else{
     6217                nt=1; //in elastic, or if we run only rigid, we need only one step
     6218        }
     6219
    61966220        constant=3/rho_earth/planetarea;
    6197        
    6198         //Allocate:
    6199         G=xNewZeroInit<IssmDouble>(3*nel);
     6221
     6222        G=xNewZeroInit<IssmDouble>(3*nel*nt);
    62006223        if(computeelastic){
    6201                 GU=xNewZeroInit<IssmDouble>(3*nel);
     6224                GU=xNewZeroInit<IssmDouble>(3*nel*nt);
    62026225                if(horiz){
    6203                         GN=xNewZeroInit<IssmDouble>(3*nel);
    6204                         GE=xNewZeroInit<IssmDouble>(3*nel);
     6226                        GN=xNewZeroInit<IssmDouble>(3*nel*nt);
     6227                        GE=xNewZeroInit<IssmDouble>(3*nel*nt);
    62056228                }
    62066229        }
    … …  
    62246247                        _assert_(index>0 && index<M);
    62256248
    6226                         /*Rigid earth gravitational perturbation: */
    6227                         G[i*nel+e]+=G_rigid_precomputed[index];
    6228                        
    6229                         if(computeelastic){
    6230                                 G[i*nel+e]+=G_elastic_precomputed[index];
    6231                         }
    6232                         G[i*nel+e]=G[i*nel+e]*ratioe;
    6233 
    6234                         /*Elastic components:*/
    6235                         if(computeelastic){
    6236                                 GU[i*nel+e] =  ratioe * U_elastic_precomputed[index];
    6237                                 if(horiz){
    6238                                         /*Compute azimuths, both north and east components: */
    6239                                         x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
    6240                                         if(lati==M_PI/2){
    6241                                                 x=1e-12; y=1e-12;
     6249                        if(horiz){
     6250                                /*Compute azimuths, both north and east components: */
     6251                                x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
     6252                                if(lati==M_PI/2){
     6253                                        x=1e-12; y=1e-12;
     6254                                }
     6255                                if(lati==-M_PI/2){
     6256                                        x=1e-12; y=1e-12;
     6257                                }
     6258                                dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z;
     6259                                N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     6260                                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     6261                        }
     6262
     6263                        for(int t=0;t<nt;t++){
     6264                                int timeindex=i*nel*nt+e*nt+t;
     6265
     6266                                /*Rigid earth gravitational perturbation: */
     6267                                G[timeindex]+=G_rigid_precomputed[index];
     6268
     6269                                /*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
     6270                                if(computeelastic){
     6271                                        G[timeindex]+=G_viscoelastic_precomputed[index*nt+t];
     6272                                }
     6273                                G[timeindex]=G[timeindex]*ratioe;
     6274
     6275                                /*Elastic components:*/
     6276                                if(computeelastic){
     6277                                        GU[timeindex] =  ratioe * U_viscoelastic_precomputed[index*nt+t];
     6278                                        if(horiz){
     6279                                                GN[timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*N_azim;
     6280                                                GE[timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*E_azim;
    62426281                                        }
    6243                                         if(lati==-M_PI/2){
    6244                                                 x=1e-12; y=1e-12;
    6245                                         }
    6246                                         dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z;
    6247                                         N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6248                                         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6249 
    6250                                         GN[i*nel+e] = ratioe*H_elastic_precomputed[index]*N_azim;
    6251                                         GE[i*nel+e] = ratioe*H_elastic_precomputed[index]*E_azim;
    62526282                                }
    6253                         }
    6254                 }
    6255         }
    6256 
     6283                        } //for(int t=0;t<nt;t++)
     6284                } //for (int i=0;i<3;i++)
     6285        } //for(int e=0;e<nel;e++)
    62576286
    62586287        /*Add in inputs:*/
    6259         this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3);
     6288        this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3*nt);
    62606289        if(computeelastic){
    6261                 this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3);
     6290                this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3*nt);
    62626291                if(horiz){
    6263                         this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3);
    6264                         this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3);
     6292                        this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3*nt);
     6293                        this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3*nt);
    62656294                }
    62666295        }
    … …  
    63466375        }
    63476376        /*}}}*/
     6377        /*Initialize stacks: {{{*/
     6378        if(computeviscous){
     6379                stackRSL=xNewZeroInit<IssmDouble>(3*nt);
     6380                stackU=xNewZeroInit<IssmDouble>(3*nt);
     6381
     6382                this->inputs->SetArrayInput(StackRSLEnum,this->lid,stackRSL,3*nt);
     6383                this->inputs->SetArrayInput(StackUEnum,this->lid,stackU,3*nt);
     6384                if(horiz){
     6385                        stackN=xNewZeroInit<IssmDouble>(3*nt);
     6386                        stackE=xNewZeroInit<IssmDouble>(3*nt);
     6387                        this->inputs->SetArrayInput(StackNEnum,this->lid,stackRSL,3*nt);
     6388                        this->inputs->SetArrayInput(StackEEnum,this->lid,stackU,3*nt);
     6389                }
     6390        }
     6391        /*}}}*/
     6392
    63486393        /*Free allocations:{{{*/
    63496394        #ifdef _HAVE_RESTRICT_
    … …  
    66926737
    66936738        #ifdef _HAVE_RESTRICT_
    6694         IssmDouble* __restrict__ G_elastic_precomputed=NULL;
     6739        IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
    66956740        IssmDouble* __restrict__ G_rigid_precomputed=NULL;
    6696         IssmDouble* __restrict__ U_elastic_precomputed=NULL;
    6697         IssmDouble* __restrict__ H_elastic_precomputed=NULL;
     6741        IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
     6742        IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
    66986743        IssmDouble** __restrict__ Gsubel=NULL;
    66996744        IssmDouble** __restrict__ GUsubel=NULL;
    … …  
    67026747
    67036748        #else
    6704         IssmDouble* G_elastic_precomputed=NULL;
     6749        IssmDouble* G_viscoelastic_precomputed=NULL;
    67056750        IssmDouble* G_rigid_precomputed=NULL;
    6706         IssmDouble* U_elastic_precomputed=NULL;
    6707         IssmDouble* H_elastic_precomputed=NULL;
     6751        IssmDouble* U_viscoelastic_precomputed=NULL;
     6752        IssmDouble* H_viscoelastic_precomputed=NULL;
    67086753        IssmDouble** Gsubel=NULL;
    67096754        IssmDouble** GUsubel=NULL;
    … …  
    67196764        bool computerigid = false;
    67206765        bool computeelastic = false;
     6766        bool computeviscous = false;
    67216767        int  horiz;
     6768
     6769        bool istime=true;
     6770        IssmDouble timeacc=0;
     6771        IssmDouble start_time,final_time;
     6772        int  nt,precomputednt;
    67226773
    67236774        /*}}}*/
    … …  
    67266777        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    67276778        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
     6779        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    67286780        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    67296781        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    … …  
    67406792
    67416793        if(computeelastic){
    6742                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter);
    6743                 parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
    6744 
    6745                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter);
    6746                 parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
    6747 
    6748                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter);
    6749                 parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
     6794                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
     6795                parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
     6796
     6797                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
     6798                parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
     6799
     6800                if(horiz){
     6801                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
     6802                        parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
     6803
     6804                }
    67506805        }
    67516806        /*}}}*/
    … …  
    67606815        constant=3/rho_earth/planetarea;
    67616816
     6817        if(computeviscous){
     6818                this->parameters->FindParam(&istime,LoveIsTimeEnum);
     6819                if(!istime)_error_("Frequency love numbers not supported yet!");
     6820                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     6821                this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
     6822                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     6823                nt=reCast<int>((final_time-start_time)/timeacc)+1;
     6824        }
     6825        else{
     6826                nt=1; //in elastic, or if we run only rigid, we need only one step
     6827        }
    67626828        Gsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    67636829        if(computeelastic){
    … …  
    67726838        for(int l=0;l<SLGEOM_NUMLOADS;l++){
    67736839                int nbar=slgeom->nbar[l];
    6774                 Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
     6840                Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    67756841                if(computeelastic){
    6776                         GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
     6842                        GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    67776843                        if(horiz){
    6778                                 GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
    6779                                 GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
     6844                                GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
     6845                                GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    67806846                        }
    67816847                }
    … …  
    67956861                                longi=longitude[i];
    67966862
     6863                                if(horiz){
     6864                                        /*Compute azimuths, both north and east components: */
     6865                                        x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
     6866                                        if(lati==90){
     6867                                                x=1e-12; y=1e-12;
     6868                                        }
     6869                                        if(lati==-90){
     6870                                                x=1e-12; y=1e-12;
     6871                                        }
     6872                                        IssmDouble xbar=planetradius*cos(late)*cos(longe);
     6873                                        IssmDouble ybar=planetradius*cos(late)*sin(longe);
     6874                                        IssmDouble zbar=planetradius*sin(late);
     6875
     6876                                        dx = xbar-x; dy = ybar-y; dz = zbar-z;
     6877                                        N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     6878                                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     6879                                }
     6880
    67976881                                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    67986882                                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    … …  
    68006884                                index=reCast<int,IssmDouble>( alpha/M_PI*reCast<IssmDouble,int>(M-1) );
    68016885
    6802                                 /*Rigid earth gravitational perturbation: */
    6803                                 Gsubel[l][i*nbar+e]+=G_rigid_precomputed[index];
    6804 
    6805                                 if(computeelastic){
    6806                                         Gsubel[l][i*nbar+e]+=G_elastic_precomputed[index];
    6807                                 }
    6808                                 Gsubel[l][i*nbar+e]*=ratioe;
    6809 
    6810                                 /*Elastic components:*/
    6811                                 if(computeelastic){
    6812                                         GUsubel[l][i*nbar+e] =  ratioe * U_elastic_precomputed[index];
    6813                                         if(horiz){
    6814                                                 /*Compute azimuths, both north and east components: */
    6815                                                 x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
    6816                                                 if(lati==90){
    6817                                                         x=1e-12; y=1e-12;
     6886
     6887                                for(int t=0;t<nt;t++){
     6888                                        int timeindex=i*nbar*nt+e*nt+t;
     6889
     6890                                        /*Rigid earth gravitational perturbation: */
     6891                                        Gsubel[l][timeindex]+=G_rigid_precomputed[index];
     6892
     6893                                        if(computeelastic){
     6894                                                Gsubel[l][timeindex]+=G_viscoelastic_precomputed[index*nt+t];
     6895                                        }
     6896                                        Gsubel[l][timeindex]*=ratioe;
     6897
     6898                                        /*Elastic components:*/
     6899                                        if(computeelastic){
     6900                                                GUsubel[l][timeindex] =  ratioe * U_viscoelastic_precomputed[index*nt+t];
     6901
     6902                                                if(horiz){
     6903                                                        GNsubel[l][timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*N_azim;
     6904                                                        GEsubel[l][timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*E_azim;
    68186905                                                }
    6819                                                 if(lati==-90){
    6820                                                         x=1e-12; y=1e-12;
    6821                                                 }
    6822                                                 IssmDouble xbar=planetradius*cos(late)*cos(longe);
    6823                                                 IssmDouble ybar=planetradius*cos(late)*sin(longe);
    6824                                                 IssmDouble zbar=planetradius*sin(late);
    6825 
    6826                                                 dx = xbar-x; dy = ybar-y; dz = zbar-z;
    6827                                                 N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6828                                                 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6829 
    6830                                                 GNsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*N_azim;
    6831                                                 GEsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*E_azim;
    68326906                                        }
    68336907                                }
    … …  
    68386912        /*Save all these arrayins for each element:*/
    68396913        for (int l=0;l<SLGEOM_NUMLOADS;l++){
    6840                 this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3);
     6914                this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3*nt);
    68416915                if(computeelastic){
    6842                         this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3);
     6916                        this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3*nt);
    68436917                        if(horiz){
    6844                                 this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3);
    6845                                 this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3);
     6918                                this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3*nt);
     6919                                this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3*nt);
    68466920                        }
    68476921                }
    … …  
    68696943        /*}}}*/
    68706944        return;
     6945
     6946}
     6947/*}}}*/
     6948void       Tria::SealevelchangeUpdateViscousFields(){ /*{{{*/
     6949       
     6950        /*Inputs:*/
     6951        IssmDouble* stackRSL=NULL;
     6952        IssmDouble* stackU=NULL;
     6953        IssmDouble* stackN=NULL;
     6954        IssmDouble* stackE=NULL;
     6955        IssmDouble* stacktimes=NULL;
     6956        int         stacknumsteps;
     6957        int         stackindex;
     6958        int         newindex;
     6959        int         dummy;
     6960        bool        viscous=false;
     6961        IssmDouble  currenttime;
     6962        IssmDouble  lincoeff=0;
     6963        int horiz;
     6964               
     6965        this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
     6966       
     6967        if(viscous){
     6968                this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     6969               
     6970                this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum);
     6971                this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum);
     6972                this->parameters->FindParam(&stackindex,StackIndexEnum);
     6973                this->parameters->FindParam(&currenttime,TimeEnum);
     6974
     6975                this->inputs->GetArrayPtr(StackRSLEnum,this->lid,&stackRSL,&dummy);
     6976                this->inputs->GetArrayPtr(StackUEnum,this->lid,&stackU,&dummy);
     6977                if(horiz){
     6978                        this->inputs->GetArrayPtr(StackNEnum,this->lid,&stackN,&dummy);
     6979                        this->inputs->GetArrayPtr(StackEEnum,this->lid,&stackE,&dummy);
     6980                }
     6981
     6982                lincoeff=0;
     6983                newindex=stacknumsteps-2;
     6984                for(int t=stackindex;t<stacknumsteps;t++){
     6985                        if (stacktimes[t]>currenttime){
     6986                                newindex=t-1;
     6987                                lincoeff=(currenttime-stacktimes[newindex])/(stacktimes[t]-stacktimes[newindex]);
     6988                                break;
     6989                        }
     6990                }
     6991                if(newindex==(stacknumsteps-2))lincoeff=1;
     6992                stacktimes[newindex]=currenttime;
     6993                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];
     6996                        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
     7003                this->parameters->SetParam(stackindex,StackIndexEnum);
     7004                this->parameters->SetParam(stacktimes,stacknumsteps,StackTimesEnum);
     7005
     7006                /*free allocations:*/
     7007                xDelete<IssmDouble>(stacktimes);
     7008        }
     7009
    68717010
    68727011}
    … …  
    69597098        IssmDouble oceanarea=0;
    69607099        IssmDouble rho_water;
     7100        bool computefuture=false;
    69617101       
    69627102        bool sal = false;
     7103        bool viscous = false;
    69637104        bool rotation= false;
    69647105        bool percpu= false;
    … …  
    69687109        IssmDouble Grotm2[3];
    69697110        IssmDouble Grotm3[3];
    6970        
     7111               
    69717112        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
     7113        this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
    69727114        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    69737115        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    … …  
    69807122                this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
    69817123
    6982                 this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true);
     7124                this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true,StackRSLEnum,computefuture=false);
    69837125        }
    69847126
    … …  
    70517193        IssmDouble* GNsub[SLGEOM_NUMLOADS];
    70527194        IssmDouble* GEsub[SLGEOM_NUMLOADS];
     7195        bool computefuture=false;
    70537196
    70547197        int horiz;
    … …  
    71027245                        }
    71037246                }
    7104 
    7105                 this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false);
     7247                this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false,StackRSLEnum,computefuture=true);
    71067248
    71077249                if(elastic){
    7108                         this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false);
     7250                        this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false,StackUEnum,computefuture=true);
    71097251                        if(horiz ){
    7110                                 this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false);
    7111                                 this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false);
     7252                                this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false,StackNEnum,computefuture=true);
     7253                                this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false,StackEEnum,computefuture=true);
    71127254                        }
    71137255                }
    … …  
    71667308
    71677309} /*}}}*/
    7168 void       Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu) { /*{{{*/
    7169 
    7170         IssmDouble sealevel[3]={0};
    7171         int i,e,l,nbar;
     7310void       Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu, int stackenum, bool computefuture) { /*{{{*/
     7311
     7312        IssmDouble* sealevel=NULL;
     7313        int i,e,l,t,nbar;
     7314        bool viscous=false;
     7315        IssmDouble* stack=NULL;
     7316        int nt=1; //important
     7317        int stackindex=0; //important
     7318        int stacknumsteps=1; //important
     7319
     7320        this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
     7321        if(viscous){
     7322                this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum);
     7323                if(computefuture) nt=stacknumsteps;
     7324                else nt=1;
     7325
     7326                //allocate
     7327                sealevel=xNewZeroInit<IssmDouble>(3*nt);
     7328        }
     7329        else sealevel=xNewZeroInit<IssmDouble>(3*nt);
    71727330
    71737331        if(loads->sealevelloads){
     7332
    71747333                for(i=0;i<NUMVERTICES;i++) {
    71757334                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    71767335                        for (e=0;e<nel;e++){
    7177                                 sealevel[i]+=G[i*nel+e]*(loads->sealevelloads[e]+loads->loads[e]);
     7336                                for(t=0;t<nt;t++){
     7337                                        int index=i*nel*stacknumsteps+e*stacknumsteps+t;
     7338                                        sealevel[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]);
     7339                                }
    71787340                        }
    71797341                        for(l=0;l<SLGEOM_NUMLOADS;l++){
    71807342                                nbar=slgeom->nbar[l];
    71817343                                for (e=0;e<nbar;e++){
    7182                                         sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subloads[l][e]);
     7344                                        for(t=0;t<nt;t++){
     7345                                                int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
     7346                                                sealevel[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
     7347                                        }
    71837348                                }
    71847349                                if(l==SLGEOM_OCEAN){
    71857350                                        for (e=0;e<nbar;e++){
    7186                                                 sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subsealevelloads[e]);
     7351                                                for(t=0;t<nt;t++){
     7352                                                        int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
     7353                                                        sealevel[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]);
     7354                                                }
    71877355                                        }
    71887356                                }
    … …  
    71917359        }
    71927360        else{  //this is the initial convolution where only loads are provided
     7361
    71937362                for(i=0;i<NUMVERTICES;i++) {
    71947363                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    71957364                        for (e=0;e<nel;e++){
    7196                                 sealevel[i]+=G[i*nel+e]*(loads->loads[e]);
     7365                                for(t=0;t<nt;t++){
     7366                                        int index=i*nel*stacknumsteps+e*stacknumsteps+t;
     7367                                        sealevel[i*nt+t]+=G[index]*(loads->loads[e]);
     7368                                }
    71977369                        }
    71987370                        for(l=0;l<SLGEOM_NUMLOADS;l++){
    71997371                                nbar=slgeom->nbar[l];
    72007372                                for (e=0;e<nbar;e++){
    7201                                         sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subloads[l][e]);
     7373                                        for(t=0;t<nt;t++){
     7374                                                int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
     7375                                                sealevel[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
     7376                                        }
    72027377                                }
    72037378                        }
    72047379                }
     7380        }
     7381       
     7382        if(viscous){
     7383                IssmDouble* sealevelinterp=NULL;
     7384                IssmDouble* stacktimes=NULL;
     7385                IssmDouble  final_time;
     7386                IssmDouble  lincoeff;
     7387                IssmDouble  timeacc;
     7388
     7389                this->parameters->FindParam(&stackindex,StackIndexEnum);
     7390                this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum);
     7391                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     7392                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     7393                this->inputs->GetArrayPtr(stackenum,this->lid,&stack,NULL);
     7394                if(computefuture){
     7395                        sealevelinterp=xNew<IssmDouble>(3*nt);
     7396                        if(stacktimes[stackindex]<final_time){
     7397                                lincoeff=1-(stacktimes[stackindex+1]-stacktimes[stackindex])/timeacc;
     7398                                for(int t=stackindex;t<nt;t++){
     7399                                        if(t==stackindex){
     7400                                                for(int i=0;i<NUMVERTICES;i++){
     7401                                                        sealevelinterp[i*nt+t]=  sealevel[i*nt+0];
     7402                                                }
     7403                                        }
     7404                                        else{
     7405                                                for(int i=0;i<NUMVERTICES;i++){
     7406                                                        sealevelinterp[i*nt+t]=  (1-lincoeff)*sealevel[i*nt+(t-stackindex-1)]+lincoeff*sealevel[i*nt+(t-stackindex)];
     7407                                                }
     7408                                        }
     7409                                }
     7410                        }
     7411                }
     7412
     7413                /*update sealevel at present time using stack at present time: */
     7414                for(int i=0;i<NUMVERTICES;i++){
     7415                        sealevel[i*nt+0]+=stack[i*stacknumsteps+stackindex];
     7416                }
     7417
     7418                if(computefuture){ /*update stack with future deformation from present load: */
     7419                        for(int t=stackindex;t<nt;t++){
     7420                                for(int i=0;i<NUMVERTICES;i++){
     7421                                        stack[i*stacknumsteps+t]+=sealevelinterp[i*nt+t];
     7422                                }
     7423                        }
     7424                        /*Re-add stack now that we updated:*/
     7425                        this->inputs->SetArrayInput(stackenum,this->lid,stack,3*stacknumsteps);
     7426                }
     7427
     7428                /*Free allocatoins:*/
     7429                xDelete<IssmDouble>(stacktimes);
    72057430        }
    72067431
    … …  
    72087433        if(percpu){
    72097434                for(i=0;i<NUMVERTICES;i++){
    7210                         if(slgeom->lids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i];
     7435                        if(slgeom->lids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0];
    72117436                }
    72127437        }
    72137438        else{
    7214                 for(i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i];
     7439                for(i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i*nt+0];
    72157440        }
    72167441

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

    r26230 r26272  
    180180                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom);
    181181                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom);
     182                void       SealevelchangeUpdateViscousFields();
    182183                #endif
    183184                /*}}}*/
    … …  
    243244                void           UpdateConstraintsExtrudeFromBase(void);
    244245                void           UpdateConstraintsExtrudeFromTop(void);
    245                 void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool optimize);
     246                void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu,int stackenum,bool computefuture);
    246247                /*}}}*/
    247248

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

    r26249 r26272  
    10191019        doubletype*  LoveHf = NULL;
    10201020        doubletype*  LoveLf = NULL;
     1021       
    10211022        doubletype*  LoveKernels= NULL;
    10221023        LoveVariables* vars=NULL;
    … …  
    11041105        //Temporal love numbers
    11051106        if (istemporal && !complex_computation){
     1107
     1108                IssmDouble*  LoveHtDouble=NULL;
     1109                IssmDouble*  LoveKtDouble=NULL;
     1110                IssmDouble*  LoveLtDouble=NULL;
     1111                doubletype*  LoveHt=NULL;
     1112                doubletype*  LoveLt=NULL;
     1113                doubletype*  LoveKt=NULL;
     1114
    11061115                int NTit;
    11071116                femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
    11081117                int nt = nfreq/2/NTit;
    11091118
    1110                 doubletype* LoveHt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
    1111                 doubletype* LoveLt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
    1112                 doubletype* LoveKt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
     1119                LoveHt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
     1120                LoveLt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
     1121                LoveKt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
    11131122
    11141123                love_freq_to_temporal<doubletype>(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel);
    11151124
     1125                /*Downcast and add into parameters:*/
     1126                LoveHtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
     1127                LoveLtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
     1128                LoveKtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
     1129                for(int i=0;i<(sh_nmax+1)*nt;i++){
     1130                        LoveHtDouble[i]=std::real(LoveHt[i]);
     1131                        LoveLtDouble[i]=std::real(LoveLt[i]);
     1132                        LoveKtDouble[i]=std::real(LoveKt[i]);
     1133                }
     1134                if(forcing_type==9){ //tidal loading
     1135                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LoveHtDouble,(sh_nmax+1)*nt,1));
     1136                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LoveKtDouble,(sh_nmax+1)*nt,1));
     1137                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LoveLtDouble,(sh_nmax+1)*nt,1));
     1138                }
     1139                else if(forcing_type==11){ //surface loading
     1140                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LoveHtDouble,(sh_nmax+1)*nt,1));
     1141                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LoveKtDouble,(sh_nmax+1)*nt,1));
     1142                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LoveLtDouble,(sh_nmax+1)*nt,1));
     1143                }
     1144                xDelete<IssmDouble>(LoveHtDouble);
     1145                xDelete<IssmDouble>(LoveKtDouble);
     1146                xDelete<IssmDouble>(LoveLtDouble);
     1147       
     1148                /*Add into external results, no need to downcast, we can handle complexes/quad precision: */
    11161149                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0));
    11171150                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0));
    … …  
    11221155                xDelete<doubletype>(LoveKt);
    11231156        }
    1124 
     1157        else{
     1158
     1159                IssmDouble*  LoveHfDouble=NULL;
     1160                IssmDouble*  LoveKfDouble=NULL;
     1161                IssmDouble*  LoveLfDouble=NULL;
     1162
     1163                /*Add into parameters:*/
     1164                LoveHfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
     1165                LoveLfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
     1166                LoveKfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
     1167                for(int i=0;i<(sh_nmax+1)*nfreq;i++){
     1168                        LoveHfDouble[i]=std::real(LoveHf[i]);
     1169                        LoveLfDouble[i]=std::real(LoveLf[i]);
     1170                        LoveKfDouble[i]=std::real(LoveKf[i]);
     1171                }
     1172                if(forcing_type==9){ //tidal loading
     1173                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LoveHfDouble,(sh_nmax+1)*nfreq,1));
     1174                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LoveKfDouble,(sh_nmax+1)*nfreq,1));
     1175                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LoveLfDouble,(sh_nmax+1)*nfreq,1));
     1176                }
     1177                else if(forcing_type==11){ //surface loading
     1178                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LoveHfDouble,(sh_nmax+1)*nfreq,1));
     1179                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LoveKfDouble,(sh_nmax+1)*nfreq,1));
     1180                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LoveLfDouble,(sh_nmax+1)*nfreq,1));
     1181                }
     1182                xDelete<IssmDouble>(LoveHfDouble);
     1183                xDelete<IssmDouble>(LoveKfDouble);
     1184                xDelete<IssmDouble>(LoveLfDouble);
     1185        }
     1186
     1187        /*Add into external results:*/
    11251188        femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0));
    11261189        femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0));

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

    r26230 r26272  
    245245        int  grdmodel;
    246246        int  computesealevel=0;
     247        bool viscous=false;
    247248        IssmDouble*           sealevelpercpu=NULL;
    248249
    … …  
    259260        femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
    260261        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
     262        femmodel->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
    261263        /*}}}*/
    262264
    … …  
    296298
    297299        if(VerboseSolution()) _printf0_("         starting  GRD convolutions\n");
     300       
     301        /*update viscous RSL:*/
     302        for(Object* & object : femmodel->elements->objects){
     303                Element* element = xDynamicCast<Element*>(object);
     304                element->SealevelchangeUpdateViscousFields();
     305        }
    298306
    299307        /*buildup loads: */

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r26271 r26272  
    364364syn keyword cConstant RotationalAngularVelocityEnum
    365365syn keyword cConstant SolidearthSettingsElasticEnum
     366syn keyword cConstant SolidearthSettingsViscousEnum
    366367syn keyword cConstant SealevelchangeGeometryDoneEnum
    367368syn keyword cConstant RotationalEquatorialMoiEnum
    … …  
    373374syn keyword cConstant LoadLoveKEnum
    374375syn keyword cConstant LoadLoveLEnum
     376syn keyword cConstant LoveTimeFreqEnum
     377syn keyword cConstant LoveIsTimeEnum
    375378syn keyword cConstant SealevelchangeGRigidEnum
    376 syn keyword cConstant SealevelchangeGElasticEnum
     379syn keyword cConstant SealevelchangeGViscoElasticEnum
    377380syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum
    378381syn keyword cConstant SolidearthSettingsGRDEnum
    379382syn keyword cConstant SolidearthSettingsRunFrequencyEnum
    380 syn keyword cConstant SealevelchangeHElasticEnum
     383syn keyword cConstant SolidearthSettingsTimeAccEnum
     384syn keyword cConstant SealevelchangeHViscoElasticEnum
    381385syn keyword cConstant SolidearthSettingsHorizEnum
    382386syn keyword cConstant SolidearthSettingsMaxiterEnum
    … …  
    389393syn keyword cConstant SealevelchangeRunCountEnum
    390394syn keyword cConstant SealevelchangeTransitionsEnum
    391 syn keyword cConstant SealevelchangeUElasticEnum
     395syn keyword cConstant SealevelchangeUViscoElasticEnum
    392396syn keyword cConstant SettingsIoGatherEnum
    393397syn keyword cConstant SettingsNumResultsOnNodesEnum
    … …  
    398402syn keyword cConstant SettingsSolverResidueThresholdEnum
    399403syn keyword cConstant SettingsWaitonlockEnum
     404syn keyword cConstant StackNumStepsEnum
     405syn keyword cConstant StackTimesEnum
     406syn keyword cConstant StackIndexEnum
    400407syn keyword cConstant SmbAIceEnum
    401408syn keyword cConstant SmbAIdxEnum
    … …  
    949956syn keyword cConstant SolidearthExternalGeoidRateEnum
    950957syn keyword cConstant SolidearthExternalBarystaticSeaLevelRateEnum
     958syn keyword cConstant StackRSLEnum
     959syn keyword cConstant StackUEnum
     960syn keyword cConstant StackNEnum
     961syn keyword cConstant StackEEnum
    951962syn keyword cConstant StrainRateeffectiveEnum
    952963syn keyword cConstant StrainRateparallelEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r26271 r26272  
    358358        RotationalAngularVelocityEnum,
    359359        SolidearthSettingsElasticEnum,
     360        SolidearthSettingsViscousEnum,
    360361        SealevelchangeGeometryDoneEnum,
    361362        RotationalEquatorialMoiEnum,
    … …  
    367368        LoadLoveKEnum,
    368369        LoadLoveLEnum,
     370        LoveTimeFreqEnum,
     371        LoveIsTimeEnum,
    369372        SealevelchangeGRigidEnum,
    370         SealevelchangeGElasticEnum,
     373        SealevelchangeGViscoElasticEnum,
    371374        SolidearthSettingsComputesealevelchangeEnum,
    372375        SolidearthSettingsGRDEnum,
    373376        SolidearthSettingsRunFrequencyEnum,
    374         SealevelchangeHElasticEnum,
     377        SolidearthSettingsTimeAccEnum,
     378        SealevelchangeHViscoElasticEnum,
    375379        SolidearthSettingsHorizEnum,
    376380        SolidearthSettingsMaxiterEnum,
    … …  
    383387        SealevelchangeRunCountEnum,
    384388        SealevelchangeTransitionsEnum,
    385         SealevelchangeUElasticEnum,
     389        SealevelchangeUViscoElasticEnum,
    386390        SettingsIoGatherEnum,
    387391        SettingsNumResultsOnNodesEnum,
    … …  
    392396        SettingsSolverResidueThresholdEnum,
    393397        SettingsWaitonlockEnum,
     398        StackNumStepsEnum,
     399        StackTimesEnum,
     400        StackIndexEnum,
    394401        SmbAIceEnum,
    395402        SmbAIdxEnum,
    … …  
    946953        SolidearthExternalGeoidRateEnum,
    947954        SolidearthExternalBarystaticSeaLevelRateEnum,
     955        StackRSLEnum,
     956        StackUEnum,
     957        StackNEnum,
     958        StackEEnum,
    948959        StrainRateeffectiveEnum,
    949960        StrainRateparallelEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r26271 r26272  
    366366                case RotationalAngularVelocityEnum : return "RotationalAngularVelocity";
    367367                case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic";
     368                case SolidearthSettingsViscousEnum : return "SolidearthSettingsViscous";
    368369                case SealevelchangeGeometryDoneEnum : return "SealevelchangeGeometryDone";
    369370                case RotationalEquatorialMoiEnum : return "RotationalEquatorialMoi";
    … …  
    375376                case LoadLoveKEnum : return "LoadLoveK";
    376377                case LoadLoveLEnum : return "LoadLoveL";
     378                case LoveTimeFreqEnum : return "LoveTimeFreq";
     379                case LoveIsTimeEnum : return "LoveIsTime";
    377380                case SealevelchangeGRigidEnum : return "SealevelchangeGRigid";
    378                 case SealevelchangeGElasticEnum : return "SealevelchangeGElastic";
     381                case SealevelchangeGViscoElasticEnum : return "SealevelchangeGViscoElastic";
    379382                case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange";
    380383                case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD";
    381384                case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency";
    382                 case SealevelchangeHElasticEnum : return "SealevelchangeHElastic";
     385                case SolidearthSettingsTimeAccEnum : return "SolidearthSettingsTimeAcc";
     386                case SealevelchangeHViscoElasticEnum : return "SealevelchangeHViscoElastic";
    383387                case SolidearthSettingsHorizEnum : return "SolidearthSettingsHoriz";
    384388                case SolidearthSettingsMaxiterEnum : return "SolidearthSettingsMaxiter";
    … …  
    391395                case SealevelchangeRunCountEnum : return "SealevelchangeRunCount";
    392396                case SealevelchangeTransitionsEnum : return "SealevelchangeTransitions";
    393                 case SealevelchangeUElasticEnum : return "SealevelchangeUElastic";
     397                case SealevelchangeUViscoElasticEnum : return "SealevelchangeUViscoElastic";
    394398                case SettingsIoGatherEnum : return "SettingsIoGather";
    395399                case SettingsNumResultsOnNodesEnum : return "SettingsNumResultsOnNodes";
    … …  
    400404                case SettingsSolverResidueThresholdEnum : return "SettingsSolverResidueThreshold";
    401405                case SettingsWaitonlockEnum : return "SettingsWaitonlock";
     406                case StackNumStepsEnum : return "StackNumSteps";
     407                case StackTimesEnum : return "StackTimes";
     408                case StackIndexEnum : return "StackIndex";
    402409                case SmbAIceEnum : return "SmbAIce";
    403410                case SmbAIdxEnum : return "SmbAIdx";
    … …  
    951958                case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
    952959                case SolidearthExternalBarystaticSeaLevelRateEnum : return "SolidearthExternalBarystaticSeaLevelRate";
     960                case StackRSLEnum : return "StackRSL";
     961                case StackUEnum : return "StackU";
     962                case StackNEnum : return "StackN";
     963                case StackEEnum : return "StackE";
    953964                case StrainRateeffectiveEnum : return "StrainRateeffective";
    954965                case StrainRateparallelEnum : return "StrainRateparallel";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r26271 r26272  
    372372              else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum;
    373373              else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum;
     374              else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum;
    374375              else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum;
    375376              else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
    … …  
    381382              else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum;
    382383              else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum;
    383               else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum;
    384               else if (strcmp(name,"SealevelchangeGElastic")==0) return SealevelchangeGElasticEnum;
     384              else if (strcmp(name,"LoveTimeFreq")==0) return LoveTimeFreqEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
     388              if (strcmp(name,"LoveIsTime")==0) return LoveIsTimeEnum;
     389              else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum;
     390              else if (strcmp(name,"SealevelchangeGViscoElastic")==0) return SealevelchangeGViscoElasticEnum;
     391              else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
    389392              else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
    390393              else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    391               else if (strcmp(name,"SealevelchangeHElastic")==0) return SealevelchangeHElasticEnum;
     394              else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum;
     395              else if (strcmp(name,"SealevelchangeHViscoElastic")==0) return SealevelchangeHViscoElasticEnum;
    392396              else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
    393397              else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
    … …  
    400404              else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
    401405              else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum;
    402               else if (strcmp(name,"SealevelchangeUElastic")==0) return SealevelchangeUElasticEnum;
     406              else if (strcmp(name,"SealevelchangeUViscoElastic")==0) return SealevelchangeUViscoElasticEnum;
    403407              else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
    404408              else if (strcmp(name,"SettingsNumResultsOnNodes")==0) return SettingsNumResultsOnNodesEnum;
    … …  
    409413              else if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum;
    410414              else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
     415              else if (strcmp(name,"StackNumSteps")==0) return StackNumStepsEnum;
     416              else if (strcmp(name,"StackTimes")==0) return StackTimesEnum;
     417              else if (strcmp(name,"StackIndex")==0) return StackIndexEnum;
    411418              else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum;
    412419              else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
    … …  
    499506              else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
    500507              else if (strcmp(name,"TimesteppingInterpForcing")==0) return TimesteppingInterpForcingEnum;
    501               else if (strcmp(name,"TimesteppingCycleForcing")==0) return TimesteppingCycleForcingEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"TimesteppingCycleForcing")==0) return TimesteppingCycleForcingEnum;
    502512              else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
    503513              else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
    … …  
    506516              else if (strcmp(name,"TimesteppingType")==0) return TimesteppingTypeEnum;
    507517              else if (strcmp(name,"ToMITgcmComm")==0) return ToMITgcmCommEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
     518              else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
    512519              else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
    513520              else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
    … …  
    622629              else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
    623630              else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
    624               else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
    625635              else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
    626636              else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
    … …  
    629639              else if (strcmp(name,"DistanceToCalvingfront")==0) return DistanceToCalvingfrontEnum;
    630640              else if (strcmp(name,"DistanceToGroundingline")==0) return DistanceToGroundinglineEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
     641              else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
    635642              else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum;
    636643              else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
    … …  
    745752              else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
    746753              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
    747               else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
    748758              else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
    749759              else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
    … …  
    752762              else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
    753763              else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
     764              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    758765              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    759766              else if (strcmp(name,"Node")==0) return NodeEnum;
    … …  
    868875              else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
    869876              else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
    870               else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
    871881              else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
    872882              else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    … …  
    875885              else if (strcmp(name,"SmbC")==0) return SmbCEnum;
    876886              else if (strcmp(name,"SmbCcsnowValue")==0) return SmbCcsnowValueEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
     887              else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
    881888              else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
    882889              else if (strcmp(name,"SmbD")==0) return SmbDEnum;
    … …  
    972979              else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
    973980              else if (strcmp(name,"SolidearthExternalBarystaticSeaLevelRate")==0) return SolidearthExternalBarystaticSeaLevelRateEnum;
     981              else if (strcmp(name,"StackRSL")==0) return StackRSLEnum;
     982              else if (strcmp(name,"StackU")==0) return StackUEnum;
     983              else if (strcmp(name,"StackN")==0) return StackNEnum;
     984              else if (strcmp(name,"StackE")==0) return StackEEnum;
    974985              else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
    975986              else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
    … …  
    987998              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    988999              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    989               else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
    9901004              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    9911005              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
    … …  
    9981012              else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
    9991013              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
     1014              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    10041015              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    10051016              else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
    … …  
    11101121              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    11111122              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
    1112               else if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
    11131127              else if (strcmp(name,"Outputdefinition65")==0) return Outputdefinition65Enum;
    11141128              else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum;
    … …  
    11211135              else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
    11221136              else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
     1137              else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
    11271138              else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
    11281139              else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
    … …  
    12331244              else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    12341245              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    1235               else if (strcmp(name,"Element")==0) return ElementEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Element")==0) return ElementEnum;
    12361250              else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
    12371251              else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
    … …  
    12441258              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    12451259              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
     1260              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    12501261              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
    12511262              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
    … …  
    13561367              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    13571368              else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
    1358               else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    13591373              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    13601374              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    … …  
    13671381              else if (strcmp(name,"Melange")==0) return MelangeEnum;
    13681382              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
     1383              else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
    13731384              else if (strcmp(name,"MeshX")==0) return MeshXEnum;
    13741385              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    … …  
    14791490              else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
    14801491              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
    1481               else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
    14821496              else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
    14831497              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    … …  
    14901504              else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
    14911505              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
     1506              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    14961507              else if (strcmp(name,"Tetra")==0) return TetraEnum;
    14971508              else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;

  • issm/trunk-jpl/src/m/classes/lovenumbers.m

    r26233 r26272  
    99classdef lovenumbers
    1010        properties (SetAccess=public)
    11                
     11
    1212                %regular love numbers:
    1313                h           = []; %provided by PREM model
    … …  
    2020                tl          = [];
    2121                tk2secular  = 0;  %deg 2 secular number.
     22
     23                %time/frequency for visco-elastic love numbers
     24                timefreq    = [];
     25                istime      = 1;
    2226
    2327        end
    … …  
    4145                        %secular fluid love number:
    4246                        self.tk2secular=0.942;
    43 
     47                       
     48                        %time:
     49                        self.istime=1; %temporal love numbers by default
     50                        self.timefreq=0; %elastic case by default.
    4451
    4552                end % }}}
    … …  
    5966                        md = checkfield(md,'fieldname','solidearth.lovenumbers.tl','NaN',1,'Inf',1);
    6067                        md = checkfield(md,'fieldname','solidearth.lovenumbers.tk2secular','NaN',1,'Inf',1);
     68                        md = checkfield(md,'fieldname','solidearth.lovenumbers.timefreq','NaN',1,'Inf',1);
     69                        md = checkfield(md,'fieldname','solidearth.lovenumbers.istime','NaN',1,'Inf',1,'values',[0 1]);
    6170
    6271                        %check that love numbers are provided at the same level of accuracy:
    … …  
    6574                        end
    6675                       
     76                        ntf=length(self.timefreq);
     77                        if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf ),
     78                                error('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector');
     79                        end
    6780
    6881                end % }}}
    … …  
    8194                        fielddisplay(self,'tl','tidal load Love number (deg 2)');
    8295                        fielddisplay(self,'tk2secular','secular fluid Love number');
     96                       
     97                        fielddisplay(self,'istime','time (default=1) or frequency love numbers (0)');
     98                        fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
    8399
    84100                end % }}}
    … …  
    94110                        WriteData(fid,prefix,'object',self,'data',self.tk2secular,'fieldname','lovenumbers.tk2secular','format','Double');
    95111
     112                        if self.istime,
     113                                scale=md.constants.yts;
     114                        else
     115                                scale=1.0/md.constants.yts;
     116                        end
     117                        WriteData(fid,prefix,'object',self,'fieldname','istime','name','md.solidearth.lovenumbers.istime','format','Boolean');
     118                        WriteData(fid,prefix,'object',self,'fieldname','timefreq','name','md.solidearth.lovenumbers.timefreq','format','DoubleMat','mattype',1,'scale',scale);
     119
    96120                end % }}}
    97121                function savemodeljs(self,fid,modelname) % {{{
    … …  
    99123                        writejs1Darray(fid,[modelname '.lovenumbers.k'],self.k);
    100124                        writejs1Darray(fid,[modelname '.lovenumbers.l'],self.l);
     125                        writejs1Darray(fid,[modelname '.lovenumbers.istime'],self.istime);
     126                        writejs1Darray(fid,[modelname '.lovenumbers.time'],self.time);
    101127                end % }}}
    102128                function self = extrude(self,md) % {{{

  • issm/trunk-jpl/src/m/classes/solidearthsettings.m

    r26229 r26272  
    1111                rigid                  = 0;
    1212                elastic                = 0;
     13                viscous                = 0;
    1314                rotation               = 0;
    1415                ocean_area_scaling     = 0;
    … …  
    1718                isgrd                  = 0; %will GRD patterns be computed?
    1819                compute_bp_grd         = 0; %will GRD patterns for bottom pressures be computed?
    19                 degacc                 = 0; %degree increment for resolution of Green tables
     20                degacc                 = 0; %degree increment for resolution of Green tables.
     21                timeacc                = 0; %time step accurary required to compute Green tables
    2022                horiz                  = 0; %compute horizontal deformation
    2123                grdmodel               = 0; %grd model (0 by default, 1 for elastic, 2 for Ivins)
    … …  
    4345                self.rigid=1;
    4446                self.elastic=1;
     47                self.viscous=1;
    4548                self.rotation=1;
    4649                self.ocean_area_scaling=0;
    … …  
    5154                %numerical discretization accuracy
    5255                self.degacc=.01;
     56                self.timeacc=1;
    5357
    5458                %how many time steps we skip before we run solidearthsettings solver during transient
    … …  
    7579                        md = checkfield(md,'fieldname','solidearth.settings.runfrequency','size',[1 1],'>=',1);
    7680                        md = checkfield(md,'fieldname','solidearth.settings.degacc','size',[1 1],'>=',1e-10);
     81                        md = checkfield(md,'fieldname','solidearth.settings.timeacc','size',[1 1],'>',0);
    7782                        md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]);
    7883                        md = checkfield(md,'fieldname','solidearth.settings.grdmodel','>=',0,'<=',2);
    … …  
    8287                        if self.elastic==1 & self.rigid==0,
    8388                                error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set');
     89                        end
     90                        if self.viscous==1 & self.elastic==0,
     91                                error('solidearthsettings checkconsistency error message: need elastic on if viscous flag is set');
    8492                        end
    8593
    … …  
    101109                        end
    102110
    103 
    104111                end % }}}
    105112                function disp(self) % {{{
    … …  
    116123                        fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
    117124                        fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
     125                        fielddisplay(self,'viscous','viscous earth graviational potential perturbation');
    118126                        fielddisplay(self,'rotation','earth rotational potential perturbation');
    119127                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
     128                        fielddisplay(self,'timeacc','time accuracy (default 1 yr) for numerical discretization of the Green''s functions');
    120129                        fielddisplay(self,'grdmodel','type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins');
    121130                        fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
    … …  
    127136                        WriteData(fid,prefix,'object',self,'fieldname','rigid','name','md.solidearth.settings.rigid','format','Boolean');
    128137                        WriteData(fid,prefix,'object',self,'fieldname','elastic','name','md.solidearth.settings.elastic','format','Boolean');
     138                        WriteData(fid,prefix,'object',self,'fieldname','viscous','name','md.solidearth.settings.viscous','format','Boolean');
    129139                        WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.rotation','format','Boolean');
    130140                        WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','name','md.solidearth.settings.ocean_area_scaling','format','Boolean');
    131141                        WriteData(fid,prefix,'object',self,'fieldname','runfrequency','name','md.solidearth.settings.runfrequency','format','Integer');
    132142                        WriteData(fid,prefix,'object',self,'fieldname','degacc','name','md.solidearth.settings.degacc','format','Double');
     143                        WriteData(fid,prefix,'object',self,'fieldname','timeacc','name','md.solidearth.settings.timeacc','format','Double','scale',md.constants.yts);
    133144                        WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer');
    134145                        WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer');
    … …  
    145156                        writejsdouble(fid,[modelname '.solidearth.settings.rigid'],self.rigid);
    146157                        writejsdouble(fid,[modelname '.solidearth.settings.elastic'],self.elastic);
     158                        writejsdouble(fid,[modelname '.solidearth.settings.viscous'],self.viscous);
    147159                        writejsdouble(fid,[modelname '.solidearth.settings.rotation'],self.rotation);
    148160                        writejsdouble(fid,[modelname '.solidearth.settings.ocean_area_scaling'],self.ocean_area_scaling);
    149161                        writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency);
    150162                        writejsdouble(fid,[modelname '.solidearth.settings.degacc'],self.degacc);
     163                        writejsdouble(fid,[modelname '.solidearth.settings.timeacc'],self.timeacc);
    151164                        writejsdouble(fid,[modelname '.solidearth.settings.cross_section_shape'],self.cross_section_shape);
    152165                end % }}}
Note: See TracChangeset for help on using the changeset viewer.