Changeset 25751


Ignore:
Timestamp:
11/13/20 15:14:27 (4 years ago)
Author:
Eric.Larour
Message:

CHG: new partition vector for solidearthsettings. New handling of rigid and
elastic table precomputations. New Europa planet radius. New earth default in
model constructor.

Location:
issm/trunk-jpl/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r25655 r25751  
    177177        bool            rigid=false;
    178178        bool            elastic=false;
     179        bool            rotation=false;
    179180
    180181        int     numoutputs;
    … …  
    203204        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum));
    204205        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum));
     206        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum));
    205207        parameters->AddObject(new DoubleParam(CumBslrEnum,0.0));
    206208        parameters->AddObject(new DoubleParam(CumBslrIceEnum,0.0));
    … …  
    212214        planetarea=4*PI*planetradius*planetradius;
    213215        parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
     216
     217        /*Deal with partition of the barystatic contribution:*/
     218        iomodel->FetchData(&npartice,"md.solidearth.npartice");
     219        parameters->AddObject(new IntParam(SolidearthNpartIceEnum,npartice));
     220        if(npartice){
     221                iomodel->FetchData(&partitionice,&nel,NULL,"md.solidearth.partitionice");
     222                parameters->AddObject(new DoubleMatParam(SolidearthPartitionIceEnum,partitionice,nel,1));
     223                parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.npartice",SolidearthNpartIceEnum));
     224                bslrice_partition=xNewZeroInit<IssmDouble>(npartice);
     225                parameters->AddObject(new DoubleMatParam(CumBslrIcePartitionEnum,bslrice_partition,npartice,1));
     226                xDelete<IssmDouble>(partitionice);
     227        }
     228        iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro");
     229        parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro));
     230        if(nparthydro){
     231                iomodel->FetchData(&partitionhydro,&nel,NULL,"md.solidearth.partitionhydro");
     232                parameters->AddObject(new DoubleMatParam(SolidearthPartitionHydroEnum,partitionhydro,nel,1));
     233                parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.nparthydro",SolidearthNpartHydroEnum));
     234                bslrhydro_partition=xNewZeroInit<IssmDouble>(nparthydro);
     235                parameters->AddObject(new DoubleMatParam(CumBslrHydroPartitionEnum,bslrhydro_partition,nparthydro,1));
     236                xDelete<IssmDouble>(partitionhydro);
     237        }
     238
    214239
    215240        /*Deal with dsl multi-model ensembles: {{{*/
    … …  
    234259        iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
    235260        iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
    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
     261        iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
     262
     263        if(elastic | rigid){
     264                /*compute green functions for a range of angles*/
     265                iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
     266                M=reCast<int,IssmDouble>(180./degacc+1.);
     267        }
     268
     269        /*love numbers: */
     270        if(elastic){
     271                iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
     272                iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
     273                iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
     274                iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
     275                iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
     276                iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
     277
     278                parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
     279                parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
     280                parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
     281                parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
     282                parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
     283                parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
     284
     285                // AD performance is sensitive to calls to ensurecontiguous.
     286                // // Providing "t" will cause ensurecontiguous to be called.
     287                #ifdef _HAVE_AD_
     288                G_elastic=xNew<IssmDouble>(M,"t");
     289                U_elastic=xNew<IssmDouble>(M,"t");
     290                H_elastic=xNew<IssmDouble>(M,"t");
     291                #else
     292                G_elastic=xNew<IssmDouble>(M);
     293                U_elastic=xNew<IssmDouble>(M);
     294                H_elastic=xNew<IssmDouble>(M);
     295                #endif
     296        }
     297        if(rigid){
     298                #ifdef _HAVE_AD_
     299                G_rigid=xNew<IssmDouble>(M,"t");
     300                #else
     301                G_rigid=xNew<IssmDouble>(M);
     302                #endif
     303        }
     304       
     305        if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
     306
     307        if(rigid | elastic){
     308
     309                /*compute combined legendre + love number (elastic green function:*/
     310                m=DetermineLocalSize(M,IssmComm::GetComm());
     311                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
     312        }
     313        if(elastic){
     314                #ifdef _HAVE_AD_
     315                G_elastic_local=xNew<IssmDouble>(m,"t");
     316                U_elastic_local=xNew<IssmDouble>(m,"t");
     317                H_elastic_local=xNew<IssmDouble>(m,"t");
     318                #else
     319                G_elastic_local=xNew<IssmDouble>(m);
     320                U_elastic_local=xNew<IssmDouble>(m);
     321                H_elastic_local=xNew<IssmDouble>(m);
     322                #endif
     323        }
     324        if(rigid){
     325                #ifdef _HAVE_AD_
     326                G_rigid_local=xNew<IssmDouble>(m,"t");
     327                #else
     328                G_rigid_local=xNew<IssmDouble>(m);
     329                #endif
     330        }
     331
     332        if(rigid){
     333                for(int i=lower_row;i<upper_row;i++){
     334                        IssmDouble alpha,x;
     335                        alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     336                        G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
     337                }
     338        }
     339        if(elastic){
     340                for(int i=lower_row;i<upper_row;i++){
     341                        IssmDouble alpha,x;
     342                        alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     343
     344                        G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
     345                        U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
     346                        H_elastic_local[i-lower_row]= 0;
     347                        IssmDouble Pn = 0.;
     348                        IssmDouble Pn1 = 0.;
     349                        IssmDouble Pn2 = 0.;
     350                        IssmDouble Pn_p = 0.;
     351                        IssmDouble Pn_p1 = 0.;
     352                        IssmDouble Pn_p2 = 0.;
     353
     354                        for (int n=0;n<nl;n++) {
     355                                IssmDouble deltalove_G;
     356                                IssmDouble deltalove_U;
     357
     358                                deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
     359                                deltalove_U = (love_h[n]-love_h[nl-1]);
     360
     361                                /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
     362                                if(n==0){
     363                                        Pn=1;
     364                                        Pn_p=0;
    330365                                }
     366                                else if(n==1){
     367                                        Pn = cos(alpha);
     368                                        Pn_p = 1;
     369                                }
     370                                else{
     371                                        Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
     372                                        Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
     373                                }
     374                                Pn2=Pn1; Pn1=Pn;
     375                                Pn_p2=Pn_p1; Pn_p1=Pn_p;
     376
     377                                G_elastic_local[i-lower_row] += deltalove_G*Pn;         // gravitational potential
     378                                U_elastic_local[i-lower_row] += deltalove_U*Pn;         // vertical (up) displacement
     379                                H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;              // horizontal displacements
    331380                        }
    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());
     381                }
     382        }
     383        if(rigid){
     384
     385                /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
     386                int* recvcounts=xNew<int>(IssmComm::GetSize());
     387                int* displs=xNew<int>(IssmComm::GetSize());
     388
     389                //recvcounts:
     390                ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     391
     392                /*displs: */
     393                ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     394
     395                /*All gather:*/
     396                ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     397                if(elastic){
    345398                        ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    346399                        ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    347400                        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];
     401                }
     402               
     403                /*free resources: */
     404                xDelete<int>(recvcounts);
     405                xDelete<int>(displs);
     406
     407                /*Avoid singularity at 0: */
     408                G_rigid[0]=G_rigid[1];
     409                G_elastic[0]=G_elastic[1];
     410                U_elastic[0]=U_elastic[1];
     411                H_elastic[0]=H_elastic[1];
     412
     413                /*Save our precomputed tables into parameters*/
     414                parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
     415                if(elastic){
    359416                        parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
    360                         U_elastic[0]=U_elastic[1];
    361417                        parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
    362                         H_elastic[0]=H_elastic[1];
    363418                        parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
    364 
    365                         /*free resources: */
     419                }
     420
     421                /*free resources: */
     422                xDelete<IssmDouble>(G_rigid);
     423                xDelete<IssmDouble>(G_rigid_local);
     424                if(elastic){
    366425                        xDelete<IssmDouble>(love_h);
    367426                        xDelete<IssmDouble>(love_k);
    … …  
    370429                        xDelete<IssmDouble>(love_tk);
    371430                        xDelete<IssmDouble>(love_tl);
    372                         xDelete<IssmDouble>(G_rigid);
    373                         xDelete<IssmDouble>(G_rigid_local);
    374431                        xDelete<IssmDouble>(G_elastic);
    375432                        xDelete<IssmDouble>(G_elastic_local);
    … …  
    378435                        xDelete<IssmDouble>(H_elastic);
    379436                        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         }
    442        
     437                }
     438        }
     439        /*}}}*/
     440
    443441        /*Indicate we have not yet run the Geometry Core module: */
    444442        parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));

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

    r25727 r25751  
    54635463        this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    54645464        this->AddInput(SealevelAreaEnum,&area,P0Enum);
     5465        if(!computerigid)return;
    54655466
    54665467        /*recover precomputed green function kernels:*/
    5467         if(computerigid){
    5468                 DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
    5469                 parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
     5468        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
     5469        parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
     5470
     5471        /*allocate indices:*/
     5472        indices=xNew<IssmDouble>(gsize);
     5473        G=xNewZeroInit<IssmDouble>(gsize);
     5474
     5475        if(computeelastic){
     5476                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
     5477                parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
     5478
     5479                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
     5480                parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
     5481
     5482                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
     5483                parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
    54705484
    54715485                /*allocate indices:*/
    5472                 indices=xNew<IssmDouble>(gsize);
    5473                 G=xNewZeroInit<IssmDouble>(gsize);
    5474 
     5486                GU=xNewZeroInit<IssmDouble>(gsize);
     5487                if(horiz){
     5488                        GN=xNewZeroInit<IssmDouble>(gsize);
     5489                        GE=xNewZeroInit<IssmDouble>(gsize);
     5490                }
     5491        }
     5492
     5493        /* Where is the centroid of this element:{{{*/
     5494
     5495        /*retrieve coordinates: */
     5496        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5497
     5498        /*figure out gravity center of our element (Cartesian): */
     5499        x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
     5500        y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
     5501        z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
     5502
     5503        /*compute gravity center in lat,long: */
     5504        late= asin(z_element/planetradius);
     5505        longe = atan2(y_element,x_element);
     5506        /*}}}*/
     5507
     5508        constant=3/rho_earth*area/planetarea;
     5509        for(int i=0;i<gsize;i++){
     5510                IssmDouble alpha;
     5511                IssmDouble delPhi,delLambda;
     5512
     5513                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     5514                lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     5515                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
     5516                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     5517                indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
     5518                index=reCast<int,IssmDouble>(indices[i]);
     5519
     5520                /*Rigid earth gravitational perturbation: */
     5521                G[i]+=G_rigid_precomputed[index];
     5522                if(elastic) G[i]+=G_elastic_precomputed[index];
     5523                G[i]=G[i]*constant;
     5524
     5525                /*Elastic components:*/
    54755526                if(computeelastic){
    5476                         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    5477                         parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
    5478 
    5479                         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
    5480                         parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
    5481 
    5482                         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
    5483                         parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
    5484 
    5485                         /*allocate indices:*/
    5486                         GU=xNewZeroInit<IssmDouble>(gsize);
     5527                        GU[i] =  constant * U_elastic_precomputed[index];
    54875528                        if(horiz){
    5488                                 GN=xNewZeroInit<IssmDouble>(gsize);
    5489                                 GE=xNewZeroInit<IssmDouble>(gsize);
    5490                         }
    5491 
    5492                         /* Where is the centroid of this element:*/
    5493 
    5494                         /*retrieve coordinates: */
    5495                         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5496 
    5497                         /*figure out gravity center of our element (Cartesian): */
    5498                         x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
    5499                         y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
    5500                         z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
    5501 
    5502                         /*compute gravity center in lat,long: */
    5503                         late= asin(z_element/planetradius);
    5504                         longe = atan2(y_element,x_element);
    5505 
    5506                         constant=3/rho_earth*area/planetarea;
    5507 
    5508                         for(int i=0;i<gsize;i++){
    5509                                 IssmDouble alpha;
    5510                                 IssmDouble delPhi,delLambda;
    5511 
    5512                                 /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    5513                                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    5514                                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
    5515                                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    5516                                 indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
    5517                                 index=reCast<int,IssmDouble>(indices[i]);
    5518 
    5519                                 /*Rigid earth gravitational perturbation: */
    5520                                 G[i]+=G_rigid_precomputed[index];
    5521                                 G[i]+=G_elastic_precomputed[index];
    5522                                 G[i]=G[i]*constant;
    5523 
    5524                                 /*Elastic components:*/
    5525                                 GU[i] =  constant * U_elastic_precomputed[index];
    5526                                 if(horiz){
    5527                                         /*Compute azimuths, both north and east components: */
    5528                                         x = xx[i]; y = yy[i]; z = zz[i];
    5529                                         if(latitude[i]==90){
    5530                                                 x=1e-12; y=1e-12;
    5531                                         }
    5532                                         if(latitude[i]==-90){
    5533                                                 x=1e-12; y=1e-12;
    5534                                         }
    5535                                         dx = x_element-x; dy = y_element-y; dz = z_element-z;
    5536                                         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);
    5537                                         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    5538 
    5539                                         GN[i] = constant*H_elastic_precomputed[index]*N_azim;
    5540                                         GE[i] = constant*H_elastic_precomputed[index]*E_azim;
     5529                                /*Compute azimuths, both north and east components: */
     5530                                x = xx[i]; y = yy[i]; z = zz[i];
     5531                                if(latitude[i]==90){
     5532                                        x=1e-12; y=1e-12;
    55415533                                }
    5542                         }
    5543 
    5544                         /*Add in inputs:*/
    5545                     this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
    5546                     this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
    5547                     this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
    5548                     if(horiz){
    5549                                 this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
    5550                                 this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
    5551                         }
    5552                         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    5553                         this->AddInput(SealevelAreaEnum,&area,P0Enum);
    5554 
    5555                         /*Free allocations:*/
    5556                         #ifdef _HAVE_RESTRICT_
    5557                         delete GU;
    5558                         if(horiz){
    5559                                 delete GN;
    5560                                 delete GE;
    5561                         }
    5562                         #else
    5563                         xDelete(GU);
    5564                         if(horiz){
    5565                                 xDelete(GN);
    5566                                 xDelete(GE);
    5567                         }
    5568                         #endif
    5569                 } else { /*computeelastic==false*/
    5570                         /* Where is the centroid of this element:*/
    5571 
    5572                         /*retrieve coordinates: */
    5573                         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5574 
    5575                         /*figure out gravity center of our element (Cartesian): */
    5576                         x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
    5577                         y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
    5578                         z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
    5579 
    5580                         /*compute gravity center in lat,long: */
    5581                         late= asin(z_element/planetradius);
    5582                         longe = atan2(y_element,x_element);
    5583 
    5584                         constant=3/rho_earth*area/planetarea;
    5585                         for(int i=0;i<gsize;i++){
    5586                                 IssmDouble alpha;
    5587                                 IssmDouble delPhi,delLambda;
    5588 
    5589                                 /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    5590                                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    5591                                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
    5592                                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    5593                                 indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
    5594                                 index=reCast<int,IssmDouble>(indices[i]);
    5595 
    5596                                 /*Rigid earth gravitational perturbation: */
    5597                                 G[i]+=G_rigid_precomputed[index];
    5598                                 G[i]=G[i]*constant;
    5599                         }
    5600 
    5601                         /*Add in inputs:*/
    5602                     this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
    5603                     this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
    5604                         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    5605                         this->AddInput(SealevelAreaEnum,&area,P0Enum);
    5606                 }
    5607 
    5608                 /*Free allocations:*/
    5609                 #ifdef _HAVE_RESTRICT_
    5610                 delete indices;
    5611                 delete G;
    5612                 #else
    5613                 xDelete(indices);
    5614                 xDelete(G);
    5615                 #endif
    5616         } else { /*computerigid==false*/
    5617                 // do nothing
    5618         }
     5534                                if(latitude[i]==-90){
     5535                                        x=1e-12; y=1e-12;
     5536                                }
     5537                                dx = x_element-x; dy = y_element-y; dz = z_element-z;
     5538                                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);
     5539                                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     5540
     5541                                GN[i] = constant*H_elastic_precomputed[index]*N_azim;
     5542                                GE[i] = constant*H_elastic_precomputed[index]*E_azim;
     5543                        }
     5544                }
     5545        }
     5546
     5547        /*Add in inputs:*/
     5548        this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
     5549        this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
     5550        if(computeelastic){
     5551                this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
     5552                if(horiz){
     5553                        this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
     5554                        this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
     5555                }
     5556        }
     5557        this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
     5558        this->AddInput(SealevelAreaEnum,&area,P0Enum);
     5559
     5560        /*Free allocations:*/
     5561        #ifdef _HAVE_RESTRICT_
     5562        delete G;
     5563        if(computeelastic){
     5564                delete GU;
     5565                if(horiz){
     5566                        delete GN;
     5567                        delete GE;
     5568                }
     5569        }
     5570        #else
     5571        xDelete(G);
     5572        if(elastic){
     5573                xDelete(GU);
     5574                if(horiz){
     5575                        xDelete(GN);
     5576                        xDelete(GE);
     5577                }
     5578        }
     5579        #endif
    56195580
    56205581        return;
    … …  
    56315592        bool scaleoceanarea= false;
    56325593        bool computerigid= false;
     5594        int  glfraction=1;
    56335595
    56345596        /*output: */
    … …  
    56905652
    56915653                phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
     5654                if(glfraction==0)phi=1;
     5655                #ifdef _ISSM_DEBUG_
     5656                this->AddInput2(SealevelEustaticMaskEnum,&phi,P0Enum);
     5657                #endif
    56925658        }
    56935659        else phi=1.0;
    … …  
    57375703        }
    57385704
     5705        /*Plug bslrice into barystatic contribution vector:*/
     5706        if(barystatic_contribution){
     5707                id=partition[this->Sid()]+1;
     5708                barystatic_contribution->SetValue(id,bslrice,ADD_VAL);
     5709        }
    57395710        /*return :*/
    57405711        return bslrice;
    … …  
    57625733        IssmDouble bslrhydro = 0;
    57635734
    5764         /*early return if we are on an ice cap:*/
    5765         if(masks->isiceonly[this->lid]){
     5735        /*early return if we are on an ice cap. Nop, we may have hydro
     5736         * and ice on the same masscon:*/
     5737        /*if(masks->isiceonly[this->lid]){
    57665738                bslrhydro=0;
    57675739                return bslrhydro;
    5768         }
     5740        }*/
    57695741
    57705742        /*early return if we are fully floating:*/
    … …  
    58075779        }
    58085780
     5781        /*Plug bslrice into barystatic contribution vector:*/
     5782        if(barystatic_contribution){
     5783                id=partition[this->Sid()]+1;
     5784                barystatic_contribution->SetValue(id,bslrhydro,ADD_VAL);
     5785        }
    58095786        /*output:*/
    58105787        return bslrhydro;

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

    r25741 r25751  
    114114syn keyword cConstant CumBslrIceEnum
    115115syn keyword cConstant CumBslrHydroEnum
     116syn keyword cConstant CumBslrIcePartitionEnum
     117syn keyword cConstant CumBslrHydroPartitionEnum
    116118syn keyword cConstant CumGmtslrEnum
    117119syn keyword cConstant CumGmslrEnum
    … …  
    328330syn keyword cConstant ModelnameEnum
    329331syn keyword cConstant SaveResultsEnum
     332syn keyword cConstant SolidearthPartitionIceEnum
     333syn keyword cConstant SolidearthPartitionHydroEnum
     334syn keyword cConstant SolidearthNpartIceEnum
     335syn keyword cConstant SolidearthNpartHydroEnum
    330336syn keyword cConstant SolidearthPlanetRadiusEnum
    331337syn keyword cConstant SolidearthPlanetAreaEnum
    … …  
    345351syn keyword cConstant SealevelriseGElasticEnum
    346352syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum
     353syn keyword cConstant SolidearthSettingsGlfractionEnum
    347354syn keyword cConstant SolidearthSettingsRunFrequencyEnum
    348355syn keyword cConstant SealevelriseHElasticEnum

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

    r25741 r25751  
    108108        CumBslrIceEnum,
    109109        CumBslrHydroEnum,
     110        CumBslrIcePartitionEnum,
     111        CumBslrHydroPartitionEnum,
    110112        CumGmtslrEnum,
    111113        CumGmslrEnum,
    … …  
    322324        ModelnameEnum,
    323325        SaveResultsEnum,
     326        SolidearthPartitionIceEnum,
     327        SolidearthPartitionHydroEnum,
     328        SolidearthNpartIceEnum,
     329        SolidearthNpartHydroEnum,
    324330        SolidearthPlanetRadiusEnum,
    325331        SolidearthPlanetAreaEnum,
    … …  
    339345        SealevelriseGElasticEnum,
    340346        SolidearthSettingsComputesealevelchangeEnum,
     347        SolidearthSettingsGlfractionEnum,
    341348        SolidearthSettingsRunFrequencyEnum,
    342349        SealevelriseHElasticEnum,

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

    r25741 r25751  
    116116                case CumBslrIceEnum : return "CumBslrIce";
    117117                case CumBslrHydroEnum : return "CumBslrHydro";
     118                case CumBslrIcePartitionEnum : return "CumBslrIcePartition";
     119                case CumBslrHydroPartitionEnum : return "CumBslrHydroPartition";
    118120                case CumGmtslrEnum : return "CumGmtslr";
    119121                case CumGmslrEnum : return "CumGmslr";
    … …  
    330332                case ModelnameEnum : return "Modelname";
    331333                case SaveResultsEnum : return "SaveResults";
     334                case SolidearthPartitionIceEnum : return "SolidearthPartitionIce";
     335                case SolidearthPartitionHydroEnum : return "SolidearthPartitionHydro";
     336                case SolidearthNpartIceEnum : return "SolidearthNpartIce";
     337                case SolidearthNpartHydroEnum : return "SolidearthNpartHydro";
    332338                case SolidearthPlanetRadiusEnum : return "SolidearthPlanetRadius";
    333339                case SolidearthPlanetAreaEnum : return "SolidearthPlanetArea";
    … …  
    347353                case SealevelriseGElasticEnum : return "SealevelriseGElastic";
    348354                case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange";
     355                case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction";
    349356                case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency";
    350357                case SealevelriseHElasticEnum : return "SealevelriseHElastic";

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

    r25741 r25751  
    116116              else if (strcmp(name,"CumBslrIce")==0) return CumBslrIceEnum;
    117117              else if (strcmp(name,"CumBslrHydro")==0) return CumBslrHydroEnum;
     118              else if (strcmp(name,"CumBslrIcePartition")==0) return CumBslrIcePartitionEnum;
     119              else if (strcmp(name,"CumBslrHydroPartition")==0) return CumBslrHydroPartitionEnum;
    118120              else if (strcmp(name,"CumGmtslr")==0) return CumGmtslrEnum;
    119121              else if (strcmp(name,"CumGmslr")==0) return CumGmslrEnum;
    … …  
    135137              else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
    136138              else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
    137               else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
    138               else if (strcmp(name,"DslModel")==0) return DslModelEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"DslModelid")==0) return DslModelidEnum;
     142              if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
     143              else if (strcmp(name,"DslModel")==0) return DslModelEnum;
     144              else if (strcmp(name,"DslModelid")==0) return DslModelidEnum;
    143145              else if (strcmp(name,"DslNummodels")==0) return DslNummodelsEnum;
    144146              else if (strcmp(name,"DslComputeFingerprints")==0) return DslComputeFingerprintsEnum;
    … …  
    258260              else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum;
    259261              else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum;
    260               else if (strcmp(name,"LoveR0")==0) return LoveR0Enum;
    261               else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum;
     265              if (strcmp(name,"LoveR0")==0) return LoveR0Enum;
     266              else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum;
     267              else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum;
    266268              else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
    267269              else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
    … …  
    336338              else if (strcmp(name,"Modelname")==0) return ModelnameEnum;
    337339              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
     340              else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
     341              else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum;
     342              else if (strcmp(name,"SolidearthNpartIce")==0) return SolidearthNpartIceEnum;
     343              else if (strcmp(name,"SolidearthNpartHydro")==0) return SolidearthNpartHydroEnum;
    338344              else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum;
    339345              else if (strcmp(name,"SolidearthPlanetArea")==0) return SolidearthPlanetAreaEnum;
    … …  
    353359              else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
    354360              else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
     361              else if (strcmp(name,"SolidearthSettingsGlfraction")==0) return SolidearthSettingsGlfractionEnum;
    355362              else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    356363              else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum;
    … …  
    376383              else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum;
    377384              else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
    378               else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum;
    379389              else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum;
    380390              else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum;
    … …  
    383393              else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
    384394              else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
     395              else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
    389396              else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum;
    390397              else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum;
    … …  
    499506              else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum;
    500507              else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
    501               else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    502512              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
    503513              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
    … …  
    506516              else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
    507517              else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum;
     518              else if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum;
    512519              else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum;
    513520              else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum;
    … …  
    622629              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    623630              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    624               else if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
    625635              else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum;
    626636              else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum;
    … …  
    629639              else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum;
    630640              else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"NGia")==0) return NGiaEnum;
     641              else if (strcmp(name,"NGia")==0) return NGiaEnum;
    635642              else if (strcmp(name,"NGiaRate")==0) return NGiaRateEnum;
    636643              else if (strcmp(name,"UGia")==0) return UGiaEnum;
    … …  
    745752              else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum;
    746753              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    747               else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    748758              else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum;
    749759              else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum;
    … …  
    752762              else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum;
    753763              else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"SmbA")==0) return SmbAEnum;
     764              else if (strcmp(name,"SmbA")==0) return SmbAEnum;
    758765              else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
    759766              else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
    … …  
    868875              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
    869876              else if (strcmp(name,"Area")==0) return AreaEnum;
    870               else if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum;
    871881              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
    872882              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    … …  
    875885              else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
    876886              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
     887              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    881888              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    882889              else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
    … …  
    991998              else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
    992999              else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
    993               else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
    9941004              else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
    9951005              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
    … …  
    9981008              else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum;
    9991009              else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
     1010              else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
    10041011              else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum;
    10051012              else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum;
    … …  
    11141121              else if (strcmp(name,"FemModel")==0) return FemModelEnum;
    11151122              else if (strcmp(name,"FileParam")==0) return FileParamEnum;
    1116               else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum;
    11171127              else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
    11181128              else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum;
    … …  
    11211131              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
    11221132              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
     1133              else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
    11271134              else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
    11281135              else if (strcmp(name,"Fset")==0) return FsetEnum;
    … …  
    12371244              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    12381245              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    1239               else if (strcmp(name,"MinVx")==0) return MinVxEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"MinVx")==0) return MinVxEnum;
    12401250              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    12411251              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
    … …  
    12441254              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    12451255              else if (strcmp(name,"Mpi")==0) return MpiEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
     1256              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    12501257              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    12511258              else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum;
    … …  
    13601367              else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
    13611368              else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
    1362               else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
    13631373              else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
    13641374              else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum;
    … …  
    13671377              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
    13681378              else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
     1379              else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
    13731380              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
    13741381              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;

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

    r25502 r25751  
    192192        methods
    193193                function md = model(varargin) % {{{
     194               
    194195
    195196                        switch nargin
    196197                                case 0
    197                                         md=setdefaultparameters(md);
    198                                 case 1
    199                                         error('model constructor not supported yet');
    200 
     198                                        md=setdefaultparameters(md,'earth');
    201199                                otherwise
    202                                         error('model constructor error message: 0 of 1 argument only in input.');
     200                                        options=pairoptions(varargin{:});
     201                                        planet=getfieldvalue(options,'planet','earth');
     202                                        md=setdefaultparameters(md,planet);
    203203                                end
    204204                end
    … …  
    12871287                        end
    12881288                end% }}}
    1289                 function md = setdefaultparameters(md) % {{{
     1289                function md = setdefaultparameters(md,planet) % {{{
    12901290
    12911291                        %initialize subclasses
    … …  
    12991299                        md.friction         = friction();
    13001300                        md.rifts            = rifts();
    1301                         md.solidearth       = solidearth();
     1301                        md.solidearth       = solidearth(planet);
    13021302                        md.dsl              = dsl();
    13031303                        md.timestepping     = timestepping();

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

    r25144 r25751  
    1414                requested_outputs      = {};
    1515                transitions            = {};
     16                partitionice              = [];
     17                partitionhydro             = [];
    1618        end
    1719        methods
    … …  
    1921                        switch nargin
    2022                                case 0
    21                                         self=setdefaultparameters(self);
     23                                        self=setdefaultparameters(self,earth);
     24                                case 1
     25                                        self=setdefaultparameters(self,varargin{:});
    2226                                otherwise
    23                                         error('constructor not supported');
     27                                        error('solidearth constructor error message: zero or one argument  only!');
    2428                        end
    2529                end % }}}
    26                 function self = setdefaultparameters(self) % {{{
     30                function self = setdefaultparameters(self,planet) % {{{
    2731               
    2832                %output default:
    … …  
    3236                self.transitions={};
    3337
     38                %no partitions requested for barystatic contribution:
     39                self.partitionice=[];
     40                self.partitionhydro=[];
     41
    3442                %earth radius
    35                 self.planetradius= planetradius('earth');
     43                self.planetradius= planetradius(planet);
    3644               
    3745                end % }}}
    … …  
    6270                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    6371                        fielddisplay(self,'requested_outputs','additional outputs requested');
     72                        fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
     73                        fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
    6474                        self.settings.disp();
    6575                        self.surfaceload.disp();
    … …  
    7383                        WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double');
    7484                        WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
     85               
     86                        if ~isempty(self.partitionice),
     87                                npartice=max(self.partitionice)+2;
     88                        else
     89                                npartice=0;
     90                        end
     91                        if ~isempty(self.partitionhydro),
     92                                nparthydro=max(self.partitionhydro)+2;
     93                        else
     94                                nparthydro=0;
     95                        end
    7596
     97                       
     98                        WriteData(fid,prefix,'object',self,'fieldname','partitionice','mattype',1,'format','DoubleMat');
     99                        WriteData(fid,prefix,'data',npartice,'format','Integer','name','md.solidearth.npartice');
     100                        WriteData(fid,prefix,'object',self,'fieldname','partitionhydro','mattype',1,'format','DoubleMat');
     101                        WriteData(fid,prefix,'data',nparthydro,'format','Integer','name','md.solidearth.nparthydro');
     102                       
    76103                        self.settings.marshall(prefix,md,fid);
    77104                        self.surfaceload.marshall(prefix,md,fid);
    … …  
    98125                        writejscellstring(fid,[modelname '.solidearth.requested_outputs'],self.requested_outputs);
    99126                        writejscellarray(fid,[modelname '.solidearth.transitions'],self.transitions);
     127                        writejscellarray(fid,[modelname '.solidearth.partition'],self.partition);
    100128                end % }}}
    101129                function self = extrude(self,md) % {{{

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

    r25151 r25751  
    6464                        md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]);
    6565
     66                        %checks on computational flags:
     67                        if self.elastic==1 & self.rigid==0,
     68                                error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set');
     69                        end
     70
    6671                        %a coupler to a planet model is provided.
    6772                        if self.computesealevelchange,

  • issm/trunk-jpl/src/m/geometry/planetradius.m

    r24902 r25751  
    1010if strcmpi(planet,'earth'),
    1111        radius=6.371012*10^6;
     12elseif strcmpi(planet,'europa'),
     13        radius=1.5008*10^6;
    1214else
    1315        error(['planet type ' planet ' not supported yet!']);
Note: See TracChangeset for help on using the changeset viewer.