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.

File:
1 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:{{{ */
Note: See TracChangeset for help on using the changeset viewer.