Changeset 25655


Ignore:
Timestamp:
10/08/20 00:06:51 (4 years ago)
Author:
jdquinn
Message:

BUG: SLR rigid/elastic fix

Location:
issm/trunk-jpl/src/c
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r25627 r25655  
    175175        IssmDouble  planetradius=0;
    176176        IssmDouble  planetarea=0;
    177         bool         elastic=false;
     177        bool            rigid=false;
     178        bool            elastic=false;
    178179
    179180        int     numoutputs;
    … …  
    231232        } /*}}}*/
    232233        /*Deal with elasticity {{{*/
     234        iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
    233235        iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
    234 // if(elastic){
    235                 /*love numbers: */
    236                 iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
    237                 iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
    238                 iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
    239                 iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
    240                 iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
    241                 iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
    242                 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
    243 
    244                 parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
    245                 parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
    246                 parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
    247                 parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
    248                 parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
    249                 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
    250 
    251                 /*compute elastic green function for a range of angles*/
    252                 iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
    253                 M=reCast<int,IssmDouble>(180./degacc+1.);
    254 
    255                 // AD performance is sensitive to calls to ensurecontiguous.
    256                 // // Providing "t" will cause ensurecontiguous to be called.
    257                 #ifdef _HAVE_AD_
    258                 G_rigid=xNew<IssmDouble>(M,"t");
    259                 G_elastic=xNew<IssmDouble>(M,"t");
    260                 U_elastic=xNew<IssmDouble>(M,"t");
    261                 H_elastic=xNew<IssmDouble>(M,"t");
    262                 #else
    263                 G_rigid=xNew<IssmDouble>(M);
    264                 G_elastic=xNew<IssmDouble>(M);
    265                 U_elastic=xNew<IssmDouble>(M);
    266                 H_elastic=xNew<IssmDouble>(M);
    267                 #endif
    268 
    269                 /*compute combined legendre + love number (elastic green function:*/
    270                 m=DetermineLocalSize(M,IssmComm::GetComm());
    271                 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    272                 // AD performance is sensitive to calls to ensurecontiguous.
    273                 // // Providing "t" will cause ensurecontiguous to be called.
    274                 #ifdef _HAVE_AD_
    275                 G_elastic_local=xNew<IssmDouble>(m,"t");
    276                 G_rigid_local=xNew<IssmDouble>(m,"t");
    277                 U_elastic_local=xNew<IssmDouble>(m,"t");
    278                 H_elastic_local=xNew<IssmDouble>(m,"t");
    279                 #else
    280                 G_elastic_local=xNew<IssmDouble>(m);
    281                 G_rigid_local=xNew<IssmDouble>(m);
    282                 U_elastic_local=xNew<IssmDouble>(m);
    283                 H_elastic_local=xNew<IssmDouble>(m);
    284                 #endif
    285 
    286                 for(int i=lower_row;i<upper_row;i++){
    287                         IssmDouble alpha,x;
    288                         alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
    289 
    290                         G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
    291                         G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
    292                         U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
    293                         H_elastic_local[i-lower_row]= 0;
    294                         IssmDouble Pn = 0.;
    295                         IssmDouble Pn1 = 0.;
    296                         IssmDouble Pn2 = 0.;
    297                         IssmDouble Pn_p = 0.;
    298                         IssmDouble Pn_p1 = 0.;
    299                         IssmDouble Pn_p2 = 0.;
    300 
    301                         for (int n=0;n<nl;n++) {
    302                                 IssmDouble deltalove_G;
    303                                 IssmDouble deltalove_U;
    304 
    305                                 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
    306                                 deltalove_U = (love_h[n]-love_h[nl-1]);
    307 
    308                                 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
    309                                 if(n==0){
    310                                         Pn=1;
    311                                         Pn_p=0;
     236        if (rigid) {
     237                if (elastic) {
     238                        /*love numbers: */
     239                        iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
     240                        iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
     241                        iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
     242                        iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
     243                        iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
     244                        iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
     245                        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
     246
     247                        parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
     248                        parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
     249                        parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
     250                        parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
     251                        parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
     252                        parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
     253
     254                        /*compute elastic green function for a range of angles*/
     255                        iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
     256                        M=reCast<int,IssmDouble>(180./degacc+1.);
     257
     258                        // AD performance is sensitive to calls to ensurecontiguous.
     259                        // // Providing "t" will cause ensurecontiguous to be called.
     260                        #ifdef _HAVE_AD_
     261                        G_rigid=xNew<IssmDouble>(M,"t");
     262                        G_elastic=xNew<IssmDouble>(M,"t");
     263                        U_elastic=xNew<IssmDouble>(M,"t");
     264                        H_elastic=xNew<IssmDouble>(M,"t");
     265                        #else
     266                        G_rigid=xNew<IssmDouble>(M);
     267                        G_elastic=xNew<IssmDouble>(M);
     268                        U_elastic=xNew<IssmDouble>(M);
     269                        H_elastic=xNew<IssmDouble>(M);
     270                        #endif
     271
     272                        /*compute combined legendre + love number (elastic green function:*/
     273                        m=DetermineLocalSize(M,IssmComm::GetComm());
     274                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
     275                        // AD performance is sensitive to calls to ensurecontiguous.
     276                        // // Providing "t" will cause ensurecontiguous to be called.
     277                        #ifdef _HAVE_AD_
     278                        G_elastic_local=xNew<IssmDouble>(m,"t");
     279                        G_rigid_local=xNew<IssmDouble>(m,"t");
     280                        U_elastic_local=xNew<IssmDouble>(m,"t");
     281                        H_elastic_local=xNew<IssmDouble>(m,"t");
     282                        #else
     283                        G_elastic_local=xNew<IssmDouble>(m);
     284                        G_rigid_local=xNew<IssmDouble>(m);
     285                        U_elastic_local=xNew<IssmDouble>(m);
     286                        H_elastic_local=xNew<IssmDouble>(m);
     287                        #endif
     288
     289                        for(int i=lower_row;i<upper_row;i++){
     290                                IssmDouble alpha,x;
     291                                alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     292
     293                                G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
     294                                G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
     295                                U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
     296                                H_elastic_local[i-lower_row]= 0;
     297                                IssmDouble Pn = 0.;
     298                                IssmDouble Pn1 = 0.;
     299                                IssmDouble Pn2 = 0.;
     300                                IssmDouble Pn_p = 0.;
     301                                IssmDouble Pn_p1 = 0.;
     302                                IssmDouble Pn_p2 = 0.;
     303
     304                                for (int n=0;n<nl;n++) {
     305                                        IssmDouble deltalove_G;
     306                                        IssmDouble deltalove_U;
     307
     308                                        deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
     309                                        deltalove_U = (love_h[n]-love_h[nl-1]);
     310
     311                                        /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
     312                                        if(n==0){
     313                                                Pn=1;
     314                                                Pn_p=0;
     315                                        }
     316                                        else if(n==1){
     317                                                Pn = cos(alpha);
     318                                                Pn_p = 1;
     319                                        }
     320                                        else{
     321                                                Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
     322                                                Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
     323                                        }
     324                                        Pn2=Pn1; Pn1=Pn;
     325                                        Pn_p2=Pn_p1; Pn_p1=Pn_p;
     326
     327                                        G_elastic_local[i-lower_row] += deltalove_G*Pn;         // gravitational potential
     328                                        U_elastic_local[i-lower_row] += deltalove_U*Pn;         // vertical (up) displacement
     329                                        H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;              // horizontal displacements
    312330                                }
    313                                 else if(n==1){
    314                                         Pn = cos(alpha);
    315                                         Pn_p = 1;
    316                                 }
    317                                 else{
    318                                         Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
    319                                         Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
    320                                 }
    321                                 Pn2=Pn1; Pn1=Pn;
    322                                 Pn_p2=Pn_p1; Pn_p1=Pn_p;
    323 
    324                                 G_elastic_local[i-lower_row] += deltalove_G*Pn;         // gravitational potential
    325                                 U_elastic_local[i-lower_row] += deltalove_U*Pn;         // vertical (up) displacement
    326                                 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;              // horizontal displacements
    327331                        }
    328                 }
    329 
    330                 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
    331                 int* recvcounts=xNew<int>(IssmComm::GetSize());
    332                 int* displs=xNew<int>(IssmComm::GetSize());
    333 
    334                 //recvcounts:
    335                 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
    336 
    337                 /*displs: */
    338                 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
    339 
    340                 /*All gather:*/
    341                 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    342                 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    343                 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    344                 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    345                
    346                 /*free resources: */
    347                 xDelete<int>(recvcounts);
    348                 xDelete<int>(displs);
    349 
    350                 /*}}}*/
    351 
    352                 /*Avoid singularity at 0: */
    353                 G_rigid[0]=G_rigid[1];
    354                 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
    355                 G_elastic[0]=G_elastic[1];
    356                 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
    357                 U_elastic[0]=U_elastic[1];
    358                 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
    359                 H_elastic[0]=H_elastic[1];
    360                 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
    361 
    362                 /*free resources: */
    363                 xDelete<IssmDouble>(love_h);
    364                 xDelete<IssmDouble>(love_k);
    365                 xDelete<IssmDouble>(love_l);
    366                 xDelete<IssmDouble>(love_th);
    367                 xDelete<IssmDouble>(love_tk);
    368                 xDelete<IssmDouble>(love_tl);
    369                 xDelete<IssmDouble>(G_rigid);
    370                 xDelete<IssmDouble>(G_rigid_local);
    371                 xDelete<IssmDouble>(G_elastic);
    372                 xDelete<IssmDouble>(G_elastic_local);
    373                 xDelete<IssmDouble>(U_elastic);
    374                 xDelete<IssmDouble>(U_elastic_local);
    375                 xDelete<IssmDouble>(H_elastic);
    376                 xDelete<IssmDouble>(H_elastic_local);
    377         // }
     332
     333                        /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
     334                        int* recvcounts=xNew<int>(IssmComm::GetSize());
     335                        int* displs=xNew<int>(IssmComm::GetSize());
     336
     337                        //recvcounts:
     338                        ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     339
     340                        /*displs: */
     341                        ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     342
     343                        /*All gather:*/
     344                        ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     345                        ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     346                        ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     347                        ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     348                       
     349                        /*free resources: */
     350                        xDelete<int>(recvcounts);
     351                        xDelete<int>(displs);
     352
     353                        /*}}}*/
     354
     355                        /*Avoid singularity at 0: */
     356                        G_rigid[0]=G_rigid[1];
     357                        parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
     358                        G_elastic[0]=G_elastic[1];
     359                        parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
     360                        U_elastic[0]=U_elastic[1];
     361                        parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
     362                        H_elastic[0]=H_elastic[1];
     363                        parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
     364
     365                        /*free resources: */
     366                        xDelete<IssmDouble>(love_h);
     367                        xDelete<IssmDouble>(love_k);
     368                        xDelete<IssmDouble>(love_l);
     369                        xDelete<IssmDouble>(love_th);
     370                        xDelete<IssmDouble>(love_tk);
     371                        xDelete<IssmDouble>(love_tl);
     372                        xDelete<IssmDouble>(G_rigid);
     373                        xDelete<IssmDouble>(G_rigid_local);
     374                        xDelete<IssmDouble>(G_elastic);
     375                        xDelete<IssmDouble>(G_elastic_local);
     376                        xDelete<IssmDouble>(U_elastic);
     377                        xDelete<IssmDouble>(U_elastic_local);
     378                        xDelete<IssmDouble>(H_elastic);
     379                        xDelete<IssmDouble>(H_elastic_local);
     380                } else { /*elastic==false*/
     381                        /*compute elastic green function for a range of angles*/
     382                        iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
     383                        M=reCast<int,IssmDouble>(180./degacc+1.);
     384                       
     385                        // AD performance is sensitive to calls to ensurecontiguous.
     386                        // // Providing "t" will cause ensurecontiguous to be called.
     387                        #ifdef _HAVE_AD_
     388                        G_rigid=xNew<IssmDouble>(M,"t");
     389                        #else
     390                        G_rigid=xNew<IssmDouble>(M);
     391                        #endif
     392
     393                        /*compute combined legendre + love number (elastic green function:*/
     394                        m=DetermineLocalSize(M,IssmComm::GetComm());
     395                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
     396                        // AD performance is sensitive to calls to ensurecontiguous.
     397                        // // Providing "t" will cause ensurecontiguous to be called.
     398                        #ifdef _HAVE_AD_
     399                        G_rigid_local=xNew<IssmDouble>(m,"t");
     400                        #else
     401                        G_rigid_local=xNew<IssmDouble>(m);
     402                        #endif
     403
     404                        for(int i=lower_row;i<upper_row;i++){
     405                                IssmDouble alpha,x;
     406                                alpha=reCast<IssmDouble>(i)*degacc*PI/180.0;
     407                                G_rigid_local[i-lower_row]=.5/sin(alpha/2.0);
     408                        }
     409
     410                        /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
     411                        int* recvcounts=xNew<int>(IssmComm::GetSize());
     412                        int* displs=xNew<int>(IssmComm::GetSize());
     413
     414                        //recvcounts:
     415                        ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     416
     417                        /*displs: */
     418                        ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     419
     420                        /*All gather:*/
     421                        ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     422                       
     423                        /*free resources: */
     424                        xDelete<int>(recvcounts);
     425                        xDelete<int>(displs);
     426
     427                        /*}}}*/
     428
     429                        /*Avoid singularity at 0: */
     430                        G_rigid[0]=G_rigid[1];
     431                        parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
     432
     433                        /*free resources: */
     434                        xDelete<IssmDouble>(G_rigid);
     435                        xDelete<IssmDouble>(G_rigid_local);
     436                }
     437        } else { /*rigid==false*/
     438                if (elastic) {
     439                        _error_("Must set md.solidearth.settings.rigid=1 when setting md.solidearth.settings.elastic=1");
     440                }
     441        }
    378442       
    379         /*Indicate we have not yet ran the Geometry Core module: */
     443        /*Indicate we have not yet run the Geometry Core module: */
    380444        parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));
    381445
    382 /*}}}*/
    383446        /*Transitions:{{{ */
    384447        iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions");

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

    r25627 r25655  
    54755475
    54765476        /*recover precomputed green function kernels:*/
    5477         DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
    5478         parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
    5479 
    5480         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    5481         parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
    5482 
    5483         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
    5484         parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
    5485 
    5486         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
    5487         parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
    5488 
    5489         /*allocate indices:*/
    5490         indices=xNew<IssmDouble>(gsize);
    5491         G=xNewZeroInit<IssmDouble>(gsize);
    5492         GU=xNewZeroInit<IssmDouble>(gsize);
    5493         if(horiz){
    5494                 GN=xNewZeroInit<IssmDouble>(gsize);
    5495                 GE=xNewZeroInit<IssmDouble>(gsize);
    5496         }
    5497 
    5498        
    5499         /* Where is the centroid of this element:*/
    5500 
    5501         /*retrieve coordinates: */
    5502         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5503 
    5504         /*figure out gravity center of our element (Cartesian): */
    5505         x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
    5506         y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
    5507         z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
    5508 
    5509         /*compute gravity center in lat,long: */
    5510         late= asin(z_element/planetradius);
    5511         longe = atan2(y_element,x_element);
    5512 
    5513         constant=3/rho_earth*area/planetarea;
    5514 
    5515         for(int i=0;i<gsize;i++){
    5516                 IssmDouble alpha;
    5517                 IssmDouble delPhi,delLambda;
    5518 
    5519                 /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    5520                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    5521                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
    5522                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    5523                 indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
    5524                 index=reCast<int,IssmDouble>(indices[i]);
    5525 
    5526                 /*Rigid earth gravitational perturbation: */
    5527                 if(computerigid){
    5528                         G[i]+=G_rigid_precomputed[index];
    5529                 }
     5477        if(computerigid){
     5478                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
     5479                parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
     5480
     5481                /*allocate indices:*/
     5482                indices=xNew<IssmDouble>(gsize);
     5483                G=xNewZeroInit<IssmDouble>(gsize);
     5484
    55305485                if(computeelastic){
    5531                         G[i]+=G_elastic_precomputed[index];
    5532                 }
    5533                 G[i]=G[i]*constant;
    5534 
    5535                 /*Elastic components:*/
    5536                 GU[i] =  constant * U_elastic_precomputed[index];
    5537                 if(horiz){
    5538                         /*Compute azimuths, both north and east components: */
    5539                         x = xx[i]; y = yy[i]; z = zz[i];
    5540                         if(latitude[i]==90){
    5541                                 x=1e-12; y=1e-12;
    5542                         }
    5543                         if(latitude[i]==-90){
    5544                                 x=1e-12; y=1e-12;
    5545                         }
    5546                         dx = x_element-x; dy = y_element-y; dz = z_element-z;
    5547                         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);
    5548                         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    5549 
    5550                         GN[i] = constant*H_elastic_precomputed[index]*N_azim;
    5551                         GE[i] = constant*H_elastic_precomputed[index]*E_azim;
    5552                 }
    5553         }
    5554 
    5555         /*Add in inputs:*/
    5556     this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
    5557     this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
    5558     this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
    5559     if(horiz){
    5560                 this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
    5561                 this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
    5562         }
    5563 
    5564         /* NOTE: May need to remove these two lines (from merge) */
    5565         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    5566         this->AddInput(SealevelAreaEnum,&area,P0Enum);
    5567 
    5568         /*Free allocations:*/
    5569         #ifdef _HAVE_RESTRICT_
    5570         delete indices;
    5571         delete G;
    5572         delete GU;
    5573         if(horiz){
    5574                 delete GN;
    5575                 delete GE;
    5576         }
    5577         #else
    5578         xDelete(indices);
    5579         xDelete(G);
    5580         xDelete(GU);
    5581         if(horiz){
    5582                 xDelete(GN);
    5583                 xDelete(GE);
    5584         }
    5585         #endif
     5486                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
     5487                        parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
     5488
     5489                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
     5490                        parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
     5491
     5492                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
     5493                        parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
     5494
     5495                        /*allocate indices:*/
     5496                        GU=xNewZeroInit<IssmDouble>(gsize);
     5497                        if(horiz){
     5498                                GN=xNewZeroInit<IssmDouble>(gsize);
     5499                                GE=xNewZeroInit<IssmDouble>(gsize);
     5500                        }
     5501
     5502                        /* Where is the centroid of this element:*/
     5503
     5504                        /*retrieve coordinates: */
     5505                        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5506
     5507                        /*figure out gravity center of our element (Cartesian): */
     5508                        x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
     5509                        y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
     5510                        z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
     5511
     5512                        /*compute gravity center in lat,long: */
     5513                        late= asin(z_element/planetradius);
     5514                        longe = atan2(y_element,x_element);
     5515
     5516                        constant=3/rho_earth*area/planetarea;
     5517
     5518                        for(int i=0;i<gsize;i++){
     5519                                IssmDouble alpha;
     5520                                IssmDouble delPhi,delLambda;
     5521
     5522                                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     5523                                lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     5524                                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
     5525                                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     5526                                indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
     5527                                index=reCast<int,IssmDouble>(indices[i]);
     5528
     5529                                /*Rigid earth gravitational perturbation: */
     5530                                if(computerigid){
     5531                                        G[i]+=G_rigid_precomputed[index];
     5532                                }
     5533                                if(computeelastic){
     5534                                        G[i]+=G_elastic_precomputed[index];
     5535                                }
     5536                                G[i]=G[i]*constant;
     5537
     5538                                /*Elastic components:*/
     5539                                GU[i] =  constant * U_elastic_precomputed[index];
     5540                                if(horiz){
     5541                                        /*Compute azimuths, both north and east components: */
     5542                                        x = xx[i]; y = yy[i]; z = zz[i];
     5543                                        if(latitude[i]==90){
     5544                                                x=1e-12; y=1e-12;
     5545                                        }
     5546                                        if(latitude[i]==-90){
     5547                                                x=1e-12; y=1e-12;
     5548                                        }
     5549                                        dx = x_element-x; dy = y_element-y; dz = z_element-z;
     5550                                        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);
     5551                                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     5552
     5553                                        GN[i] = constant*H_elastic_precomputed[index]*N_azim;
     5554                                        GE[i] = constant*H_elastic_precomputed[index]*E_azim;
     5555                                }
     5556                        }
     5557
     5558                        /*Add in inputs:*/
     5559                    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
     5560                    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
     5561                    this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
     5562                    if(horiz){
     5563                                this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
     5564                                this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
     5565                        }
     5566                        this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
     5567                        this->AddInput(SealevelAreaEnum,&area,P0Enum);
     5568
     5569                        /*Free allocations:*/
     5570                        #ifdef _HAVE_RESTRICT_
     5571                        delete GU;
     5572                        if(horiz){
     5573                                delete GN;
     5574                                delete GE;
     5575                        }
     5576                        #else
     5577                        xDelete(GU);
     5578                        if(horiz){
     5579                                xDelete(GN);
     5580                                xDelete(GE);
     5581                        }
     5582                        #endif
     5583                } else { /*computeelastic==false*/
     5584                        /* Where is the centroid of this element:*/
     5585
     5586                        /*retrieve coordinates: */
     5587                        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5588
     5589                        /*figure out gravity center of our element (Cartesian): */
     5590                        x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
     5591                        y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
     5592                        z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
     5593
     5594                        /*compute gravity center in lat,long: */
     5595                        late= asin(z_element/planetradius);
     5596                        longe = atan2(y_element,x_element);
     5597
     5598                        constant=3/rho_earth*area/planetarea;
     5599                        for(int i=0;i<gsize;i++){
     5600                                IssmDouble alpha;
     5601                                IssmDouble delPhi,delLambda;
     5602
     5603                                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     5604                                lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     5605                                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
     5606                                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     5607                                indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
     5608                                index=reCast<int,IssmDouble>(indices[i]);
     5609
     5610                                /*Rigid earth gravitational perturbation: */
     5611                                G[i]+=G_rigid_precomputed[index];
     5612                                G[i]=G[i]*constant;
     5613                        }
     5614
     5615                        /*Add in inputs:*/
     5616                    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
     5617                    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
     5618                        this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
     5619                        this->AddInput(SealevelAreaEnum,&area,P0Enum);
     5620
     5621                        /*Free allocations:*/
     5622                        #ifdef _HAVE_RESTRICT_
     5623                        delete indices;
     5624                        delete G;
     5625                        #else
     5626                        xDelete(indices);
     5627                        xDelete(G);
     5628                        #endif
     5629                }
     5630        } else { /*computerigid==false*/
     5631                // do nothing
     5632        }
    55865633
    55875634        return;
    … …  
    56465693
    56475694        /*retrieve precomputed G:*/
    5648 
    5649         /* NOTE: May need to remove condition here (from merge) */
    56505695        if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
    56515696
    … …  
    57525797
    57535798        /*retrieve precomputed G:*/
    5754 
    5755         /* NOTE: May need to remove condition here (from merge) */
    57565799        if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
    57575800
Note: See TracChangeset for help on using the changeset viewer.