Changeset 21344


Ignore:
Timestamp:
11/04/16 17:42:57 (8 years ago)
Author:
adhikari
Message:

CHG: rotational feedback works fine now

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

Legend:

Unmodified
Added
Removed

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

    r21331 r21344  
    7373        parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
    7474        parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_k",SealevelriseTidalLoveKEnum));
     75        parameters->AddObject(iomodel->CopyConstantObject("md.slr.fluid_love",SealevelriseFluidLoveEnum));
     76        parameters->AddObject(iomodel->CopyConstantObject("md.slr.equatorial_moi",SealevelriseEquatorialMoiEnum));
     77        parameters->AddObject(iomodel->CopyConstantObject("md.slr.polar_moi",SealevelrisePolarMoiEnum));
    7578        parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
    7679        parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));

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

    r21295 r21344  
    305305                virtual IssmDouble    OceanAverage(IssmDouble* Sg)=0;
    306306                virtual IssmDouble    OceanArea(void)=0;
     307                virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0;
    307308                virtual void          SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
    308309                virtual void          SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;

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

    r21295 r21344  
    190190                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    191191                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
     192                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    192193                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    193194                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};

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

    r21295 r21344  
    173173#endif
    174174#ifdef _HAVE_SEALEVELRISE_
     175                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    175176                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    176177                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};

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

    r21295 r21344  
    180180#endif
    181181#ifdef _HAVE_SEALEVELRISE_
     182                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    182183                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    183184                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};

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

    r21331 r21344  
    38703870}
    38713871/*}}}*/
     3872void    Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/
     3873
     3874        /*early return if we are not on an ice cap OR ocean:*/
     3875        if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return;
     3876
     3877        /*Compute area of element:*/
     3878        IssmDouble area;
     3879        area=GetAreaSpherical();
     3880
     3881        /*Compute lat,long,radius of elemental centroid: */
     3882        bool spherical=true;
     3883        IssmDouble llr_list[NUMVERTICES][3];
     3884        IssmDouble late,longe,re;
     3885        /* Where is the centroid of this element?:{{{*/
     3886        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
     3887
     3888        IssmDouble minlong=400;
     3889        IssmDouble maxlong=-20;
     3890        for (int i=0;i<NUMVERTICES;i++){
     3891                llr_list[i][0]=(90-llr_list[i][0]);
     3892                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
     3893                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
     3894                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
     3895        }
     3896        if(minlong==0 && maxlong>180){
     3897                if (llr_list[0][1]==0)llr_list[0][1]=360;
     3898                if (llr_list[1][1]==0)llr_list[1][1]=360;
     3899                if (llr_list[2][1]==0)llr_list[2][1]=360;
     3900        }
     3901
     3902        // correction at the north pole
     3903        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3904        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3905        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3906
     3907        //correction at the south pole
     3908        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3909        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3910        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3911
     3912        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
     3913        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
     3914
     3915        late=90-late;
     3916        if(longe>180)longe=(longe-180)-180;
     3917
     3918        late=late/180*PI;
     3919        longe=longe/180*PI;
     3920        /*}}}*/
     3921        re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
     3922
     3923        if(IsWaterInElement()){
     3924                IssmDouble rho_water, S;
     3925               
     3926                /*recover material parameters: */
     3927                rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     3928       
     3929                /*From Sg_old, recover water sea level rise:*/
     3930                S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
     3931               
     3932                dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
     3933                dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
     3934                dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(cos(late),2))/eartharea;
     3935        }
     3936        else if(this->inputs->Max(MaskIceLevelsetEnum)<0){
     3937                IssmDouble rho_ice, I;
     3938               
     3939                /*recover material parameters: */
     3940                rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     3941       
     3942                /*Compute ice thickness change: */
     3943                Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
     3944                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
     3945                deltathickness_input->GetInputAverage(&I);
     3946               
     3947                dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
     3948                dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
     3949                dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(cos(late),2))/eartharea;
     3950        }
     3951       
     3952        return;
     3953}/*}}}*/
    38723954void    Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    38733955
    … …  
    40264108        IssmDouble area;
    40274109        IssmDouble S;  //change in water water level(Farrel and Clarke, Equ. 4)
    4028         IssmDouble late,longe,re;
     4110        IssmDouble late,longe;
    40294111        IssmDouble lati,longi,ri;
    40304112        IssmDouble minlong=400;
    … …  
    40484130        bool computerigid = true;
    40494131        bool computeelastic= true;
    4050         bool computerotation= true;
    40514132
    40524133        /*early return if we are not on the ocean:*/
    … …  
    40564137        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
    40574138        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    4058         this->parameters->FindParam(&computerotation,SealevelriseRotationEnum);
    4059 
    4060         /*recover some parameters for rotational feedback: */
    4061         IssmDouble tide_love_h, tide_love_k, omega;
    4062         if(computerotation){
    4063                 this->parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum);
    4064                 this->parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum);
    4065                 this->parameters->FindParam(&omega,SealevelriseAngularVelocityEnum);
    4066         }
    40674139
    40684140        /*early return if rigid or elastic not requested:*/
    … …  
    41344206        IssmDouble delPhi,delLambda;
    41354207
    4136          /* if rotation, compute m1 m2 m3 from loading function */
    4137          IssmDouble m1, m2, m3;
    4138          if(computerotation){
    4139                  m1=0;
    4140                  m2=0;
    4141                  m3=0;
    4142          }
    4143 
    41444208        for(int i=0;i<gsize;i++){
    41454209
    … …  
    41524216                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    41534217
    4154                 //Rigid earth gravitational perturbation:
     4218                /*Rigid earth gravitational perturbation: */
    41554219                if(computerigid){
    41564220                        G_rigid[i]=1.0/2.0/sin(alpha/2.0);
    … …  
    41584222                }
    41594223
    4160                 //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
     4224                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    41614225                if(computeelastic){
    41624226                        int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    41634227                        G_elastic[i] += G_elastic_precomputed[index];
    41644228                        values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
    4165                 }
    4166 
    4167                 /*Rotational feedback: */
    4168                 if(computerotation){
    4169                         values[i]+=(1+tide_love_k-tide_love_h/9.81)*0.5*pow((omega*radius[i]),2)*(-m3*(2+m3)*(1+3*cos(2*lati))
    4170                                                 + 2*pow(m2,2)*(pow(cos(lati),2) + 0.5*(-1 + 3*cos(2*longi)) * pow(sin(lati),2))
    4171                                                 + 2*pow(m1,2)*(pow(cos(lati),2) - 0.5*(1+3*cos(2*longi)) * pow(sin(lati),2))
    4172                                                 - m1*(1+m3)*cos(longi)*sin(2*lati)
    4173                        - 2*m1*m2*cos(longi)*pow(sin(lati),2)*sin(longi)
    4174                        - m2*(1+m3)*sin(2*lati)*sin(longi));
    41754229                }
    41764230        }

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

    r21295 r21344  
    151151                IssmDouble OceanArea(void);
    152152                IssmDouble OceanAverage(IssmDouble* Sg);
     153                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea);
    153154                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    154155                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r21295 r21344  
    24592459        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    24602460        ns = elements->Size();
    2461        
     2461
    24622462        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
    24632463        for(int i=0;i<ns;i++){
    … …  
    24652465                eartharea_cpu += element->GetAreaSpherical();
    24662466        }
     2467       
    24672468        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    24682469        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    … …  
    24752476        for(int i=0;i<nsmax;i++){
    24762477                if(i<ns){
    2477 
    24782478                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
    2479 
    24802479                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2481                         element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,eartharea);
     2480                        element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
    24822481                }
    24832482                if(i%100==0)pSgo->Assemble();
    24842483        }
    24852484        if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
    2486 
     2485       
    24872486        /*Free ressources:*/
    24882487        xDelete<IssmDouble>(Sg_old);
     2488}
     2489/*}}}*/
     2490void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
     2491
     2492        /*serialized vectors:*/
     2493        IssmDouble* Sg_old=NULL;
     2494        IssmDouble  eartharea=0;
     2495        IssmDouble  eartharea_cpu=0;
     2496        IssmDouble      tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g;
     2497        IssmDouble      load_love_k2 = -0.30922675; //degree 2 load Love number
     2498        IssmDouble      m1, m2, m3;
     2499        IssmDouble      lati, longi, value;
     2500
     2501        /*Serialize vectors from previous iteration:*/
     2502        Sg_old=pSg_old->ToMPISerial();
     2503
     2504        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
     2505        for(int i=0;i<elements->Size();i++){
     2506                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2507                eartharea_cpu += element->GetAreaSpherical();
     2508        }
     2509        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2510        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2511
     2512        IssmDouble moi_list[3]={0,0,0};
     2513        IssmDouble moi_list_cpu[3]={0,0,0};
     2514        for(int i=0;i<elements->Size();i++){
     2515                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2516                element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
     2517                moi_list_cpu[0] += moi_list[0];
     2518                moi_list_cpu[1] += moi_list[1];
     2519                moi_list_cpu[2] += moi_list[2];
     2520        }
     2521        ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2522        ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2523        //     
     2524        ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2525        ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2526        //     
     2527        ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2528        ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2529       
     2530        /*pull out some useful parameters: */
     2531        parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum);
     2532        parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum);
     2533        parameters->FindParam(&fluid_love,SealevelriseFluidLoveEnum);
     2534        parameters->FindParam(&moi_e,SealevelriseEquatorialMoiEnum);
     2535        parameters->FindParam(&moi_p,SealevelrisePolarMoiEnum);
     2536        parameters->FindParam(&omega,SealevelriseAngularVelocityEnum);
     2537
     2538        /*compute perturbation terms for angular velocity vector: */
     2539        m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0];
     2540        m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1];
     2541        m3 = -(1+load_love_k2)/moi_p * moi_list[2];
     2542
     2543        for(int i=0;i<vertices->Size();i++){
     2544                //Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
     2545                Vertex* vertex=xDynamicCast<Vertex*>(vertices->GetObjectByOffset(i));
     2546               
     2547                lati=latitude[i]/180*PI;        longi=longitude[i]/180*PI;
     2548
     2549                /*only first order terms are considered now: */
     2550                value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radius[i],2.0)*
     2551                                                (-m3/6.0 - m3*cos(2.0*lati)/2.0 - sin(lati)*cos(lati)*(m1*cos(longi)+m2*sin(longi)));
     2552       
     2553                pSgo_rot->SetValue(vertex->Sid(),value,INS_VAL); //INS_VAL ensures that you don't add several times
     2554        }
     2555
     2556        /*Assemble mesh velocity*/
     2557        pSgo_rot->Assemble();
     2558       
     2559        /*Free ressources:*/
     2560        xDelete<IssmDouble>(Sg_old);
     2561       
    24892562}
    24902563/*}}}*/

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r21272 r21344  
    120120                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    121121                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
     122                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    122123                void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    123124                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
     125                IssmDouble SealevelriseRotationalFeedback(IssmDouble* moi_list);
    124126                #endif
    125127                void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);

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

    r21079 r21344  
    1818
    1919        Vector<IssmDouble> *Sgo    = NULL; //ocean convolution of the perturbation to gravity potential.
     20        Vector<IssmDouble> *Sgo_rot= NULL; // rotational feedback
    2021        IssmDouble          Sgo_oceanaverage = 0;  //average of Sgo over the ocean.
    2122
    … …  
    2728        bool spherical=true;
    2829        bool converged=true;
     30        bool rotation=true;
    2931        bool verboseconvolution=true;
    3032        int max_nonlinear_iterations;
    … …  
    4143        femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum);
    4244        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     45       
     46        /*computational flag: */
     47        femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum);
    4348
    4449        /*first, recover lat,long and radius vectors from vertices: */
    … …  
    7580                Sgo->Assemble();
    7681
     82                if(rotation){
     83                        /*call rotational feedback  module: */
     84                        Sgo_rot = new Vector<IssmDouble>(gsize); Sgo_rot->Assemble();
     85                        femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,latitude,longitude,radius);
     86                        Sgo_rot->Assemble();
     87
     88                        Sgo->AXPY(Sgo_rot,1);
     89                }
     90               
    7791                /*we need to average Sgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    7892                Sgo_oceanaverage=femmodel->SealevelriseOceanAverage(Sgo);

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

    r21331 r21344  
    775775        SealevelriseRotationEnum,
    776776        SealevelriseTidalLoveHEnum,
    777    SealevelriseTidalLoveKEnum,
     777        SealevelriseTidalLoveKEnum,
     778        SealevelriseFluidLoveEnum,
     779        SealevelriseEquatorialMoiEnum,
     780        SealevelrisePolarMoiEnum,
    778781        SealevelriseAngularVelocityEnum,
    779782        SealevelriseOceanAreaScalingEnum,

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

    r21331 r21344  
    757757                case SealevelriseTidalLoveHEnum : return "SealevelriseTidalLoveH";
    758758                case SealevelriseTidalLoveKEnum : return "SealevelriseTidalLoveK";
     759                case SealevelriseFluidLoveEnum : return "SealevelriseFluidLove";
     760                case SealevelriseEquatorialMoiEnum : return "SealevelriseEquatorialMoi";
     761                case SealevelrisePolarMoiEnum : return "SealevelrisePolarMoi";
    759762                case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity";
    760763                case SealevelriseOceanAreaScalingEnum : return "SealevelriseOceanAreaScaling";

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

    r21331 r21344  
    775775              else if (strcmp(name,"SealevelriseTidalLoveH")==0) return SealevelriseTidalLoveHEnum;
    776776              else if (strcmp(name,"SealevelriseTidalLoveK")==0) return SealevelriseTidalLoveKEnum;
     777              else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum;
     778              else if (strcmp(name,"SealevelriseEquatorialMoi")==0) return SealevelriseEquatorialMoiEnum;
     779              else if (strcmp(name,"SealevelrisePolarMoi")==0) return SealevelrisePolarMoiEnum;
    777780              else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum;
    778781              else if (strcmp(name,"SealevelriseOceanAreaScaling")==0) return SealevelriseOceanAreaScalingEnum;
    … …  
    872875              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    873876              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    874               else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    875               else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
    876               else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
     880              if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
     881              else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
     882              else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
     883              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    881884              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
    882885              else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
Note: See TracChangeset for help on using the changeset viewer.