Changeset 25763


Ignore:
Timestamp:
11/16/20 16:00:17 (4 years ago)
Author:
Eric.Larour
Message:

CHG: fixing issues with elastic vs rigid vs rotation.

Location:
issm/trunk-jpl
Files:
40 edited

Legend:

Unmodified
Added
Removed

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

    r25439 r25763  
    3232        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    3333        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    34         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     34        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    3535        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    3636        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);

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

    r25439 r25763  
    9999        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    100100        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    101         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     101        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    102102        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    103103        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);

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

    r25439 r25763  
    4141        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    4242        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    43         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     43        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    4444        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    4545        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);

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

    r25539 r25763  
    136136        iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum);
    137137        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    138         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     138        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    139139        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    140140        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);

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

    r25539 r25763  
    7474
    7575        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    76         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     76        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    7777        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    7878        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum,0.);

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

    r25439 r25763  
    7777
    7878        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    79         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     79        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    8080        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    8181        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);

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

    r25697 r25763  
    129129        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    130130        iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
    131         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); /*Needed for friction*/
     131        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); /*Needed for friction*/
    132132        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
    133133        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);

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

    r25439 r25763  
    5858        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    5959        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    60         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     60        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    6161        if(iomodel->domaintype!=Domain2DhorizontalEnum){
    6262                iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);

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

    r25439 r25763  
    4242        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    4343   iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    44         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     44        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    4545        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    4646        if(iomodel->domaintype!=Domain2DhorizontalEnum & iomodel->domaintype!=Domain3DsurfaceEnum){

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

    r25539 r25763  
    148148        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    149149        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    150         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     150        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    151151        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    152152        iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);

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

    r25439 r25763  
    5656        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    5757        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    58         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     58        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    5959        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    6060        if(iomodel->domaintype!=Domain2DhorizontalEnum){

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

    r25759 r25763  
    4040        iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange");
    4141        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum);
    42         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     42        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    4343        iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
    4444        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightChangeEnum);
    … …  
    175175        IssmDouble  planetradius=0;
    176176        IssmDouble  planetarea=0;
     177        bool            rigid=false;
    177178        bool            elastic=false;
     179        bool            rotation=false;
    178180
    179181        int     numoutputs;
    … …  
    207209        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum));
    208210        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum));
     211        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.isgrd",SolidearthSettingsGRDEnum));
    209212        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum));
    210213        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum));
    … …  
    260263        } /*}}}*/
    261264        /*Deal with elasticity {{{*/
     265        iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
    262266        iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
     267        iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
     268
     269        if(elastic | rigid){
     270                /*compute green functions for a range of angles*/
     271                iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
     272                M=reCast<int,IssmDouble>(180./degacc+1.);
     273        }
     274
     275        /*love numbers: */
    263276        if(elastic){
    264                 /*love numbers: */
    265277                iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
    266278                iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
    … …  
    269281                iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
    270282                iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
    271                 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
    272283
    273284                parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
    … …  
    278289                parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
    279290
    280                 /*compute elastic green function for a range of angles*/
    281                 iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
    282                 M=reCast<int,IssmDouble>(180./degacc+1.);
    283 
    284291                // AD performance is sensitive to calls to ensurecontiguous.
    285292                // // Providing "t" will cause ensurecontiguous to be called.
    286293                #ifdef _HAVE_AD_
    287                 G_rigid=xNew<IssmDouble>(M,"t");
    288294                G_elastic=xNew<IssmDouble>(M,"t");
    289295                U_elastic=xNew<IssmDouble>(M,"t");
    290296                H_elastic=xNew<IssmDouble>(M,"t");
    291297                #else
    292                 G_rigid=xNew<IssmDouble>(M);
    293298                G_elastic=xNew<IssmDouble>(M);
    294299                U_elastic=xNew<IssmDouble>(M);
    295300                H_elastic=xNew<IssmDouble>(M);
    296301                #endif
     302        }
     303        if(rigid){
     304                #ifdef _HAVE_AD_
     305                G_rigid=xNew<IssmDouble>(M,"t");
     306                #else
     307                G_rigid=xNew<IssmDouble>(M);
     308                #endif
     309        }
     310       
     311        if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
     312
     313        if(rigid | elastic){
    297314
    298315                /*compute combined legendre + love number (elastic green function:*/
    299316                m=DetermineLocalSize(M,IssmComm::GetComm());
    300317                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    301                 // AD performance is sensitive to calls to ensurecontiguous.
    302                 // // Providing "t" will cause ensurecontiguous to be called.
     318        }
     319        if(elastic){
    303320                #ifdef _HAVE_AD_
    304321                G_elastic_local=xNew<IssmDouble>(m,"t");
    305                 G_rigid_local=xNew<IssmDouble>(m,"t");
    306322                U_elastic_local=xNew<IssmDouble>(m,"t");
    307323                H_elastic_local=xNew<IssmDouble>(m,"t");
    308324                #else
    309325                G_elastic_local=xNew<IssmDouble>(m);
    310                 G_rigid_local=xNew<IssmDouble>(m);
    311326                U_elastic_local=xNew<IssmDouble>(m);
    312327                H_elastic_local=xNew<IssmDouble>(m);
    313328                #endif
    314 
     329        }
     330        if(rigid){
     331                #ifdef _HAVE_AD_
     332                G_rigid_local=xNew<IssmDouble>(m,"t");
     333                #else
     334                G_rigid_local=xNew<IssmDouble>(m);
     335                #endif
     336        }
     337
     338        if(rigid){
    315339                for(int i=lower_row;i<upper_row;i++){
    316340                        IssmDouble alpha,x;
    317341                        alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
    318 
    319342                        G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
     343                }
     344        }
     345        if(elastic){
     346                for(int i=lower_row;i<upper_row;i++){
     347                        IssmDouble alpha,x;
     348                        alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
     349
    320350                        G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
    321351                        U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
    … …  
    356386                        }
    357387                }
     388        }
     389        if(rigid){
    358390
    359391                /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
    … …  
    369401                /*All gather:*/
    370402                ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    371                 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    372                 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    373                 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    374 
     403                if(elastic){
     404                        ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     405                        ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     406                        ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     407                }
     408               
    375409                /*free resources: */
    376410                xDelete<int>(recvcounts);
    377411                xDelete<int>(displs);
    378                 /*}}}*/
    379412
    380413                /*Avoid singularity at 0: */
    381414                G_rigid[0]=G_rigid[1];
    382                 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
    383                 G_elastic[0]=G_elastic[1];
    384                 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
    385                 U_elastic[0]=U_elastic[1];
    386                 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
    387                 H_elastic[0]=H_elastic[1];
    388                 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
     415                if(elastic){
     416                        G_elastic[0]=G_elastic[1];
     417                        U_elastic[0]=U_elastic[1];
     418                        H_elastic[0]=H_elastic[1];
     419                }
     420               
    389421
    390422                /*Save our precomputed tables into parameters*/
    … …  
    397429
    398430                /*free resources: */
    399                 xDelete<IssmDouble>(love_h);
    400                 xDelete<IssmDouble>(love_k);
    401                 xDelete<IssmDouble>(love_l);
    402                 xDelete<IssmDouble>(love_th);
    403                 xDelete<IssmDouble>(love_tk);
    404                 xDelete<IssmDouble>(love_tl);
    405431                xDelete<IssmDouble>(G_rigid);
    406432                xDelete<IssmDouble>(G_rigid_local);
    407                 xDelete<IssmDouble>(G_elastic);
    408                 xDelete<IssmDouble>(G_elastic_local);
    409                 xDelete<IssmDouble>(U_elastic);
    410                 xDelete<IssmDouble>(U_elastic_local);
    411                 xDelete<IssmDouble>(H_elastic);
    412                 xDelete<IssmDouble>(H_elastic_local);
    413         }
     433                if(elastic){
     434                        xDelete<IssmDouble>(love_h);
     435                        xDelete<IssmDouble>(love_k);
     436                        xDelete<IssmDouble>(love_l);
     437                        xDelete<IssmDouble>(love_th);
     438                        xDelete<IssmDouble>(love_tk);
     439                        xDelete<IssmDouble>(love_tl);
     440                        xDelete<IssmDouble>(G_elastic);
     441                        xDelete<IssmDouble>(G_elastic_local);
     442                        xDelete<IssmDouble>(U_elastic);
     443                        xDelete<IssmDouble>(U_elastic_local);
     444                        xDelete<IssmDouble>(H_elastic);
     445                        xDelete<IssmDouble>(H_elastic_local);
     446                }
     447        } /*}}}*/
    414448
    415449        /*Indicate we have not yet run the Geometry Core module: */
    … …  
    436470        /*}}}*/
    437471
    438         /*}}}*/
    439472}/*}}}*/
    440473

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

    r25680 r25763  
    768768        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    769769        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    770         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     770        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    771771        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    772772        iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);

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

    r25610 r25763  
    105105        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    106106        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    107         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     107        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    108108        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    109109        if(iomodel->domaintype!=Domain2DhorizontalEnum){

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

    r25439 r25763  
    133133        iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
    134134        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    135         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
     135        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    136136        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    137137        iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);

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

    r25758 r25763  
    54635463        this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    54645464        this->AddInput(SealevelAreaEnum,&area,P0Enum);
    5465         if(!computerigid & !computeelastic)return;
     5465        if(!computerigid)return;
    54665466
    54675467        /*recover precomputed green function kernels:*/
    … …  
    54695469        parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
    54705470
    5471         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    5472         parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
    5473 
    5474         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
    5475         parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
    5476 
    5477         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
    5478         parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
     5471        if(computeelastic){
     5472                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
     5473                parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
     5474
     5475                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
     5476                parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
     5477
     5478                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
     5479                parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
     5480        }
    54795481
    54805482        /*allocate indices:*/
    54815483        indices=xNew<IssmDouble>(gsize);
    54825484        G=xNewZeroInit<IssmDouble>(gsize);
    5483         GU=xNewZeroInit<IssmDouble>(gsize);
    5484         if(horiz){
    5485                 GN=xNewZeroInit<IssmDouble>(gsize);
    5486                 GE=xNewZeroInit<IssmDouble>(gsize);
    5487         }
    5488 
    5489         /* Where is the centroid of this element:{{{*/
     5485        if(computeelastic){
     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:
    54905494
    54915495        /*retrieve coordinates: */
    … …  
    55155519
    55165520                /*Rigid earth gravitational perturbation: */
    5517                 if(computerigid){
    5518                         G[i]+=G_rigid_precomputed[index];
    5519                 }
     5521                G[i]+=G_rigid_precomputed[index];
     5522               
    55205523                if(computeelastic){
    55215524                        G[i]+=G_elastic_precomputed[index];
    … …  
    55245527
    55255528                /*Elastic components:*/
    5526                 GU[i] =  constant * U_elastic_precomputed[index];
    5527                 if(horiz){
    5528                         /*Compute azimuths, both north and east components: */
    5529                         x = xx[i]; y = yy[i]; z = zz[i];
    5530                         if(latitude[i]==90){
    5531                                 x=1e-12; y=1e-12;
    5532                         }
    5533                         if(latitude[i]==-90){
    5534                                 x=1e-12; y=1e-12;
    5535                         }
    5536                         dx = x_element-x; dy = y_element-y; dz = z_element-z;
    5537                         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);
    5538                         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    5539 
    5540                         GN[i] = constant*H_elastic_precomputed[index]*N_azim;
    5541                         GE[i] = constant*H_elastic_precomputed[index]*E_azim;
     5529                if(computeelastic){
     5530                        GU[i] =  constant * U_elastic_precomputed[index];
     5531                        if(horiz){
     5532                                /*Compute azimuths, both north and east components: */
     5533                                x = xx[i]; y = yy[i]; z = zz[i];
     5534                                if(latitude[i]==90){
     5535                                        x=1e-12; y=1e-12;
     5536                                }
     5537                                if(latitude[i]==-90){
     5538                                        x=1e-12; y=1e-12;
     5539                                }
     5540                                dx = x_element-x; dy = y_element-y; dz = z_element-z;
     5541                                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);
     5542                                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     5543
     5544                                GN[i] = constant*H_elastic_precomputed[index]*N_azim;
     5545                                GE[i] = constant*H_elastic_precomputed[index]*E_azim;
     5546                        }
    55425547                }
    55435548        }
    … …  
    55465551    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
    55475552    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
    5548     this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
    5549     if(horiz){
    5550                 this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
    5551                 this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
     5553        if(computeelastic){
     5554                this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
     5555                if(horiz){
     5556                        this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
     5557                        this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
     5558                }
    55525559        }
    55535560
    … …  
    55565563        delete indices;
    55575564        delete G;
    5558         delete GU;
    5559         if(horiz){
    5560                 delete GN;
    5561                 delete GE;
     5565        if(computeelastic){
     5566                delete GU;
     5567                if(horiz){
     5568                        delete GN;
     5569                        delete GE;
     5570                }
    55625571        }
    55635572        #else
    55645573        xDelete(indices);
    55655574        xDelete(G);
    5566         xDelete(GU);
    5567         if(horiz){
    5568                 xDelete(GN);
    5569                 xDelete(GE);
     5575        if(computeelastic){
     5576                xDelete(GU);
     5577                if(horiz){
     5578                        xDelete(GN);
     5579                        xDelete(GE);
     5580                }
    55705581        }
    55715582        #endif

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

    r25757 r25763  
    2424        bool isslr=0;
    2525        bool isgia=0;
     26        int  grd=0;
    2627        int solution_type;
    2728
    … …  
    3031        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    3132        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     33        femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
    3234
    3335        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
    … …  
    5961
    6062        /*Run geodetic:*/
    61         grd_core(femmodel);
     63        if(grd)grd_core(femmodel);
    6264
    6365        /*Run steric core for sure:*/
    … …  
    109111        int  frequency,count;
    110112        int  horiz;
    111         int  geodetic=0;
    112113        IssmDouble dt;
    113114        IssmDouble oceanarea;
    114115        int        bp_compute_fingerprints=0;
    115116
    116         /*Should we even be here?:*/
    117         femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); if(!geodetic)return;
    118 
     117       
    119118        /*Verbose: */
    120         if(VerboseSolution()) _printf0_("         computing geodetic sea level rise\n");
     119        if(VerboseSolution()) _printf0_("         computing GRD sea level patterns\n");
    121120
    122121        /*retrieve more parameters:*/
    … …  
    178177                RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea);
    179178
    180                 /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
     179                /*compute other elastic signatures, such as components of 3-D crustal motion: */
    181180                sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks);
    182181
    183                 /*compute viscosus (GIA) geodetic signatures:*/
     182                /*compute viscosus (GIA) signatures:*/
    184183                sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
    185184
    … …  
    269268        int  solution_type;
    270269        IssmDouble          dt;
    271         int  geodetic=0;
     270        int  grd=0;
    272271        int  step;
    273272        IssmDouble time;
    … …  
    281280        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    282281        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    283         femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum);
     282        femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
    284283        femmodel->parameters->FindParam(&cumgmtslr,CumGmtslrEnum);
    285284        femmodel->parameters->FindParam(&cumbslr,CumBslrEnum);
    … …  
    301300        GetStericRate(&steric_rate_g,femmodel);
    302301        GetDynamicRate(&dynamic_rate_g,femmodel);
    303         if(geodetic){
     302        if(grd){
    304303                GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
    305304                GetVectorFromInputsx(&U_gia_rate,femmodel,UGiaRateEnum,VertexSIdEnum);
    … …  
    318317
    319318        /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate dt*/
    320         if(geodetic){
     319        if(grd){
    321320                SL->AXPY(N_gia_rate,dt);
    322321                SL->AXPY(N_esa_rate,dt);
    … …  
    326325
    327326        /*compute new bedrock position: */
    328         if(geodetic){
     327        if(grd){
    329328                bedrock->AXPY(U_esa_rate,dt);
    330329                bedrock->AXPY(U_gia_rate,dt);
    … …  
    340339        delete steric_rate_g;
    341340        delete dynamic_rate_g;
    342         if(geodetic){
     341        if(grd){
    343342                delete U_esa_rate;
    344343                delete U_gia_rate;

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

    r25751 r25763  
    351351syn keyword cConstant SealevelriseGElasticEnum
    352352syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum
     353syn keyword cConstant SolidearthSettingsGRDEnum
    353354syn keyword cConstant SolidearthSettingsGlfractionEnum
    354355syn keyword cConstant SolidearthSettingsRunFrequencyEnum

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

    r25751 r25763  
    345345        SealevelriseGElasticEnum,
    346346        SolidearthSettingsComputesealevelchangeEnum,
     347        SolidearthSettingsGRDEnum,
    347348        SolidearthSettingsGlfractionEnum,
    348349        SolidearthSettingsRunFrequencyEnum,

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

    r25751 r25763  
    353353                case SealevelriseGElasticEnum : return "SealevelriseGElastic";
    354354                case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange";
     355                case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD";
    355356                case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction";
    356357                case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency";

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

    r25751 r25763  
    359359              else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
    360360              else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
     361              else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
    361362              else if (strcmp(name,"SolidearthSettingsGlfraction")==0) return SolidearthSettingsGlfractionEnum;
    362363              else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    … …  
    382383              else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
    383384              else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum;
    384               else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum;
     388              if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
     389              else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum;
    389390              else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum;
    390391              else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum;
    … …  
    505506              else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
    506507              else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum;
    507               else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
     511              if (strcmp(name,"Adjoint")==0) return AdjointEnum;
     512              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    512513              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
    513514              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
    … …  
    628629              else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    629630              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    630               else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
     634              if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
     635              else if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
    635636              else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum;
    636637              else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum;
    … …  
    751752              else if (strcmp(name,"SealevelriseGE")==0) return SealevelriseGEEnum;
    752753              else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum;
    753               else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
     757              if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
     758              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    758759              else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum;
    759760              else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum;
    … …  
    874875              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    875876              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
    876               else if (strcmp(name,"Area")==0) return AreaEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum;
     880              if (strcmp(name,"Area")==0) return AreaEnum;
     881              else if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum;
    881882              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
    882883              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    … …  
    997998              else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
    998999              else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
    999               else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
     1003              if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
     1004              else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
    10041005              else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
    10051006              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
    … …  
    11201121              else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
    11211122              else if (strcmp(name,"FemModel")==0) return FemModelEnum;
    1122               else if (strcmp(name,"FileParam")==0) return FileParamEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum;
     1126              if (strcmp(name,"FileParam")==0) return FileParamEnum;
     1127              else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum;
    11271128              else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
    11281129              else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum;
    … …  
    12431244              else if (strcmp(name,"MeshX")==0) return MeshXEnum;
    12441245              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    1245               else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"MinVx")==0) return MinVxEnum;
     1249              if (strcmp(name,"MinVel")==0) return MinVelEnum;
     1250              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    12501251              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    12511252              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
    … …  
    13661367              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    13671368              else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
    1368               else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
     1372              if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
     1373              else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
    13731374              else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
    13741375              else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum;

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

    r25758 r25763  
    66classdef solidearth
    77        properties (SetAccess=public)
    8                 sealevel               = NaN;
     8                initialsealevel               = NaN;
    99                settings               = solidearthsettings();
     10                external               = [];
    1011                surfaceload            = surfaceload();
    1112                lovenumbers            = lovenumbers();
    … …  
    3839                self.partitionhydro=[];
    3940
     41                %no external solutions by defalt:
     42                self.external=[];
     43
    4044                %earth radius
    4145                self.planetradius= planetradius('earth');
    … …  
    4852                        end
    4953
    50                         md = checkfield(md,'fieldname','solidearth.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
     54                        md = checkfield(md,'fieldname','solidearth.initialsealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    5155                        md = checkfield(md,'fieldname','solidearth.requested_outputs','stringrow',1);
    5256
    … …  
    5559                        self.lovenumbers.checkconsistency(md,solution,analyses);
    5660                        self.rotational.checkconsistency(md,solution,analyses);
     61                        if ~isempty(self.external),
     62                                if ~isa(self.external,'solidearthsolution'),
     63                                        error('solidearth consistency check: external field should be a solidearthsolution');
     64                                end
     65                                self.external.checkconsistency(md,solution,analyses);
     66                        end
    5767
    5868                       
    … …  
    6474                        disp(sprintf('   solidearth inputs, forcings and settings:'));
    6575
    66                         fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
     76                        fielddisplay(self,'initialsealevel','sea level at the start of computation) [m]');
    6777                        fielddisplay(self,'planetradius','planet radius [m]');
    6878                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    … …  
    7080                        fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
    7181                        fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
     82                        if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end
    7283                        self.settings.disp();
    7384                        self.surfaceload.disp();
    7485                        self.lovenumbers.disp();
    7586                        self.rotational.disp();
     87                        if ~isempty(self.external),
     88                                self.external.disp();
     89                        end
    7690
    7791                end % }}}
    7892                function marshall(self,prefix,md,fid) % {{{
    7993                       
    80                         WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     94                        WriteData(fid,prefix,'object',self,'fieldname','initialsealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    8195                        WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double');
    8296                        WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
    … …  
    103117                        self.lovenumbers.marshall(prefix,md,fid);
    104118                        self.rotational.marshall(prefix,md,fid);
     119                        if ~isempty(self.external),
     120                                self.external.marshall(prefix,md,fid);
     121                        end
    105122                       
    106123                        %process requested outputs
    … …  
    116133                function savemodeljs(self,fid,modelname) % {{{
    117134               
    118                         writejs1Darray(fid,[modelname '.solidearth.sealevel'],self.sealevel);
     135                        writejs1Darray(fid,[modelname '.solidearth.initialsealevel'],self.initialsealevel);
    119136                        self.settings.savemodeljs(fid,modelname);
    120137                        self.surfaceload.savemodeljs(fid,modelname);
    121138                        self.lovenumbers.savemodeljs(fid,modelname);
    122139                        self.rotational.savemodeljs(fid,modelname);
     140                        if ~isempty(self.external),
     141                                self.external.savemodeljs(fid,modelname);
     142                        end
    123143                        writejscellstring(fid,[modelname '.solidearth.requested_outputs'],self.requested_outputs);
    124144                        writejscellarray(fid,[modelname '.solidearth.transitions'],self.transitions);
    … …  
    126146                end % }}}
    127147                function self = extrude(self,md) % {{{
    128                         self.sealevel=project3d(md,'vector',self.sealevel,'type','node');
     148                        self.initialsealevel=project3d(md,'vector',self.initialsealevel,'type','node');
    129149                end % }}}
    130150        end

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

    r25758 r25763  
    1414                ocean_area_scaling     = 0;
    1515                runfrequency           = 1; %how many time steps we skip before we run grd_core
    16                 computesealevelchange  = 0; %will grd_core compute sea level?
     16                computesealevelchange  = 1; %will sea-level be coputed?
     17                isgrd                  = 1; %will GRD patterns be computed?
    1718                degacc                 = 0; %degree increment for resolution of Green tables
    1819                horiz                  = 0; %compute horizontal deformation
    … …  
    4243                self.rotation=1;
    4344                self.ocean_area_scaling=0;
    44                 self.computesealevelchange=0;
     45                self.isgrd=1;
     46                self.computesealevelchange=1;
    4547
    4648                %numerical discretization accuracy
    … …  
    7072
    7173                        %a coupler to a planet model is provided.
    72                         if self.computesealevelchange,
     74                        if self.isgrd,
    7375                                if md.transient.iscoupler,
    7476                                        %we are good;
    … …  
    9092                        fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
    9193                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
    92                         fielddisplay(self,'computesealevelchange','compute sealevel change from GRD in addition to steric?) default 0');
     94                        fielddisplay(self,'computesealevelchange','compute sealevel change (default 1)');
     95                        fielddisplay(self,'isgrd','compute GRD patterns (default 1)');
    9396                        fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
    9497                        fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
    … …  
    110113                        WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer');
    111114                        WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer');
     115                        WriteData(fid,prefix,'object',self,'fieldname','isgrd','name','md.solidearth.settings.isgrd','format','Integer');
    112116                        WriteData(fid,prefix,'object',self,'fieldname','glfraction','name','md.solidearth.settings.glfraction','format','Integer');
    113117                end % }}}

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

    r25753 r25763  
    3232                function md = checkconsistency(self,md,solution,analyses) % {{{
    3333
    34                         if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslr==0),
    35                                 return;
    36                         end
    37                         md = checkfield(md,'fieldname','solidearth.solution.displacementeast','Inf',1,'timeseries',1);
    38                         md = checkfield(md,'fieldname','solidearth.solution.displacementnorth','Inf',1,'timeseries',1);
    39                         md = checkfield(md,'fieldname','solidearth.solution.displacementup','Inf',1,'timeseries',1);
    40                         md = checkfield(md,'fieldname','solidearth.solution.geoid','Inf',1,'timeseries',1);
    41                         md = checkfield(md,'fieldname','solidearth.solution.barystaticsealevel','Inf',1,'timeseries',1);
     34                        md = checkfield(md,'fieldname','solidearth.external.displacementeast','Inf',1,'timeseries',1);
     35                        md = checkfield(md,'fieldname','solidearth.external.displacementnorth','Inf',1,'timeseries',1);
     36                        md = checkfield(md,'fieldname','solidearth.external.displacementup','Inf',1,'timeseries',1);
     37                        md = checkfield(md,'fieldname','solidearth.external.geoid','Inf',1,'timeseries',1);
     38                        md = checkfield(md,'fieldname','solidearth.external.barystaticsealevel','Inf',1,'timeseries',1);
    4239
    4340                end % }}}
    4441                function disp(self) % {{{
    45                         disp(sprintf('   solidearth solution:'));
    46 
    4742                        fielddisplay(self,'displacementeast','solid-Earth Eastwards bedrock motion (m/yr)');
    4843                        fielddisplay(self,'displacementnorth','solid-Earth Northwards bedrock motion (m/yr)');
    … …  
    5348                end % }}}
    5449                function marshall(self,prefix,md,fid) % {{{
    55                         WriteData(fid,prefix,'object',self,'data',1,'name','md.solidearth.external.nature','format','Integer'); %code 1 for father class
    5650                        yts=md.constants.yts;
    5751                        WriteData(fid,prefix,'object',self,'class','solidearthsolution','fieldname','displacementeast','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);

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

    r25181 r25763  
    88%solidearth loading:  {{{
    99md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    10 md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
     10md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    1111md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    1212md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);

  • issm/trunk-jpl/test/NightlyRun/test2002.py

    r25688 r25763  
    1919#solidearth loading:
    2020md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1))
    21 md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1))
     21md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1))
    2222md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
    2323md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))

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

    r25166 r25763  
    88%solidearth loading:  {{{
    99md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    10 md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
     10md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    1111md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    1212md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);

  • issm/trunk-jpl/test/NightlyRun/test2003.py

    r25182 r25763  
    1919#solidearth loading:  {{{
    2020md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, ))
    21 md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, ))
     21md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
    2222md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
    2323md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))

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

    r25627 r25763  
    167167                end
    168168
    169                 md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
     169                md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    170170
    171171                md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    … …  
    302302        end
    303303
    304         md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
     304        md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    305305
    306306        md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    … …  
    349349sl.transfer('mesh.long');
    350350sl.transfer('solidearth.surfaceload.icethicknesschange');
    351 sl.transfer('solidearth.sealevel');
     351sl.transfer('solidearth.initialsealevel');
    352352sl.transfer('dsl.sea_surface_height_change_above_geoid');
    353353sl.transfer('dsl.sea_water_pressure_change_at_sea_floor');

  • issm/trunk-jpl/test/NightlyRun/test2004.py

    r25517 r25763  
    206206            md.solidearth.surfaceload.icethicknesschange = np.mean(delHAIS[md.mesh.elements - 1], axis=1) / 100
    207207
    208         md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, ))
     208        md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
    209209
    210210        md.dsl.global_average_thermostatic_sea_level_change = np.zeros((2, 1))
    … …  
    335335        md.mask.ocean_levelset[pos] = 1
    336336
    337     md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, ))
     337    md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
    338338
    339339    md.dsl.global_average_thermostatic_sea_level_change = np.zeros((2, 1))
    … …  
    378378sl.transfer('mesh.long')
    379379sl.transfer('solidearth.surfaceload.icethicknesschange') #
    380 sl.transfer('solidearth.sealevel')
     380sl.transfer('solidearth.initialsealevel')
    381381sl.transfer('dsl.sea_surface_height_change_above_geoid')
    382382sl.transfer('dsl.sea_water_pressure_change_at_sea_floor')

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

    r25688 r25763  
    88%solidearth loading:  {{{
    99md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    10 md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
     10md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    1111md.dsl.global_average_thermosteric_sea_level_change=[1;0];
    1212md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);

  • issm/trunk-jpl/test/NightlyRun/test2005.py

    r25688 r25763  
    1919# Solidearth loading #{{{
    2020md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1))
    21 md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1))
     21md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1))
    2222md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
    2323md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))

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

    r25698 r25763  
    99%solidearth loading:  {{{
    1010md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    11 md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
     11md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    1212md.dsl.global_average_thermosteric_sea_level_change=[1;0];
    1313md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);

  • issm/trunk-jpl/test/NightlyRun/test2006.py

    r25698 r25763  
    2222# Solidearth loading #{{{
    2323md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1))
    24 md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1))
     24md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1))
    2525md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
    2626md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))

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

    r25627 r25763  
    1515md.solidearth.surfaceload.icethicknesschange(pos(6:7))=-1;
    1616
    17 md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
     17md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    1818md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    1919md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);

  • issm/trunk-jpl/test/NightlyRun/test2010.py

    r25651 r25763  
    2727md.solidearth.surfaceload.icethicknesschange[pos[5:7]] = -1
    2828
    29 md.solidearth.sealevel = np.zeros(md.mesh.numberofvertices)
     29md.solidearth.initialsealevel = np.zeros(md.mesh.numberofvertices)
    3030md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
    3131md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))

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

    r25124 r25763  
    2020
    2121md.timestepping.time_step=.1;
    22 md.solidearth.sealevel=newforcing(md.timestepping.start_time, md.timestepping.final_time, md.timestepping.time_step,-200,200,md.mesh.numberofvertices);
     22md.solidearth.initialsealevel=newforcing(md.timestepping.start_time, md.timestepping.final_time, md.timestepping.time_step,-200,200,md.mesh.numberofvertices);
    2323
    2424md.cluster=generic('name',oshostname(),'np',3);

  • issm/trunk-jpl/test/NightlyRun/test2424.py

    r25125 r25763  
    3030
    3131md.timestepping.time_step = .1
    32 md.solidearth.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time,
     32md.solidearth.initialsealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time,
    3333                             md.timestepping.time_step, -200., 200., md.mesh.numberofvertices)
    3434

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

    r25124 r25763  
    3131md.geometry.bed=md.geometry.bed+1000;
    3232md.geometry.surface=md.geometry.surface+1000;
    33 md.solidearth.sealevel=1000*ones(md.mesh.numberofvertices,1);
     33md.solidearth.initialsealevel=1000*ones(md.mesh.numberofvertices,1);
    3434
    3535md=solve(md,'Transient','checkconsistency','no');

  • issm/trunk-jpl/test/NightlyRun/test2425.py

    r25125 r25763  
    3939md.geometry.bed = md.geometry.bed + 1000
    4040md.geometry.surface = md.geometry.surface + 1000
    41 md.solidearth.sealevel = 1000 * np.ones((md.mesh.numberofvertices, ))
     41md.solidearth.initialsealevel = 1000 * np.ones((md.mesh.numberofvertices, ))
    4242
    4343md = solve(md, 'Transient', 'checkconsistency', 'no')
Note: See TracChangeset for help on using the changeset viewer.