Changeset 21031


Ignore:
Timestamp:
07/28/16 16:52:19 (9 years ago)
Author:
adhikari
Message:

CHG: minor changes

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r21002 r21031  
    36993699                int* indices=xNew<int>(gsize);
    37003700                IssmDouble* values=xNew<IssmDouble>(gsize);
     3701                IssmDouble alpha;
     3702                IssmDouble delPhi,delLambda;
    37013703                for(int i=0;i<gsize;i++){
    37023704                        indices[i]=i;
    37033705
    3704                         IssmDouble alpha;
    37053706                        IssmDouble G_rigid=0;  //do not remove =0!
    37063707                        IssmDouble G_elastic=0;  //do not remove =0!
    3707                         IssmDouble delPhi,delLambda;
    37083708
    37093709                        /*Compute alpha angle between centroid and current vertex : */
    37103710                        lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    37113711
    3712                     delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     3712                   delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    37133713                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    37143714
    … …  
    37443744        IssmDouble llr_list[NUMVERTICES][3];
    37453745        IssmDouble area;
    3746         IssmDouble I;  //change in water water level(Farrel and Clarke, Equ. 4)
     3746        IssmDouble S;  //change in water water level(Farrel and Clarke, Equ. 4)
    37473747        IssmDouble late,longe,re;
    37483748        IssmDouble lati,longi,ri;
    … …  
    37863786
    37873787        /*From Sg_old, recover water sea level rise:*/
    3788         I=0; for(int i=0;i<NUMVERTICES;i++) I+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
     3788        S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    37893789
    37903790        /*Compute area of element:*/
    … …  
    38403840        int* indices=xNew<int>(gsize);
    38413841        IssmDouble* values=xNewZeroInit<IssmDouble>(gsize);
     3842        IssmDouble alpha;
     3843        IssmDouble delPhi,delLambda;
    38423844
    38433845        for(int i=0;i<gsize;i++){
    38443846
    38453847                indices[i]=i;
    3846                 if(computeelastic | computerigid){
    3847        
    3848                         IssmDouble alpha;
    3849                         IssmDouble delPhi,delLambda;
    3850 
    3851                         /*Compute alpha angle between centroid and current vertex : */
    3852                         lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    3853 
    3854                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    3855                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    3856 
    3857                         //Rigid earth gravitational perturbation:
    3858                         if(computerigid)G_rigid[i]=1.0/2.0/sin(alpha/2.0);
    3859 
    3860                         //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
    3861                         if(computeelastic){
    3862                                 int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    3863                                 G_elastic[i] += G_elastic_precomputed[index];
    3864                         }
    3865                 }
    3866 
    3867                 /*Add all components to the pSgo solution vectors:*/
    3868                 if(computerigid)values[i]+=3*rho_water/rho_earth*area/eartharea*I*G_rigid[i];
    3869                 if(computeelastic)values[i]+=3*rho_water/rho_earth*area/eartharea*I*G_elastic[i];
     3848
     3849                /*Compute alpha angle between centroid and current vertex : */
     3850                lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     3851
     3852                delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     3853                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     3854
     3855                //Rigid earth gravitational perturbation:
     3856                if(computerigid){
     3857                        G_rigid[i]=1.0/2.0/sin(alpha/2.0);
     3858                        values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i];
     3859                }
     3860
     3861                //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
     3862                if(computeelastic){
     3863                        int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
     3864                        G_elastic[i] += G_elastic_precomputed[index];
     3865                        values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
     3866                }
    38703867        }
    38713868       
    … …  
    38913888        IssmDouble xyz_list[NUMVERTICES][3];
    38923889        IssmDouble area;
    3893         IssmDouble I;           //change in relative sea level or ice thickness
     3890        IssmDouble I, S;                //change in relative ice thickness and sea level
    38943891        IssmDouble late,longe,re;
    38953892        IssmDouble lati,longi,ri;
    … …  
    39223919        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    39233920       
    3924         /*early return if rigid or elastic not requested:*/
    3925         if(!computerigid && !computeelastic) return;
     3921        /*early return if elastic not requested:*/
     3922        if(!computeelastic) return;
    39263923
    39273924        /*recover material parameters: */
    … …  
    39803977        z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
    39813978
    3982         if(computeelastic){
     3979        /*recover elastic Green's functions for displacement:*/
     3980        DoubleVecParam* U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter);
     3981        DoubleVecParam* H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter);
     3982        U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
     3983        H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
     3984
     3985        /*From Sg, recover water sea level rise:*/
     3986        S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
    39833987       
    3984                 /*recover elastic Green's functions for displacement:*/
    3985                 DoubleVecParam* U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter);
    3986                 DoubleVecParam* H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter);
    3987                 U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
    3988                 H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
    3989 
    3990                 /*initialize: */
    3991                 U_elastic=xNewZeroInit<IssmDouble>(gsize);
    3992                 N_elastic=xNewZeroInit<IssmDouble>(gsize);
    3993                 E_elastic=xNewZeroInit<IssmDouble>(gsize);
    3994         }
     3988        /*Compute ice thickness change: */
     3989        Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
     3990        if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
     3991        deltathickness_input->GetInputAverage(&I);
     3992               
     3993        /*initialize: */
     3994        U_elastic=xNewZeroInit<IssmDouble>(gsize);
     3995        N_elastic=xNewZeroInit<IssmDouble>(gsize);
     3996        E_elastic=xNewZeroInit<IssmDouble>(gsize);
    39953997
    39963998        int* indices=xNew<int>(gsize);
    … …  
    39984000        IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
    39994001        IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
     4002        IssmDouble alpha;
     4003        IssmDouble delPhi,delLambda;
     4004        IssmDouble dx, dy, dz, x, y, z;
     4005        IssmDouble N_azim, E_azim;
    40004006
    40014007        for(int i=0;i<gsize;i++){
    40024008
    40034009                indices[i]=i;
    4004                 if(computeelastic){
    4005        
    4006                         IssmDouble alpha;
    4007                         IssmDouble delPhi,delLambda;
    4008 
    4009                         /*Compute alpha angle between centroid and current vertex: */
    4010                         lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    4011 
    4012                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    4013                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    4014 
    4015                         /*Compute azimuths, both north and east components: */
    4016                         IssmDouble dx, dy, dz, x, y, z;
    4017                         IssmDouble N_azim, E_azim;
    4018                         x = xx[i]; y = yy[i]; z = zz[i];
    4019                         if(latitude[i]==90){
    4020                                 x=1e-12; y=1e-12;
    4021                         }
    4022                         if(latitude[i]==-90){
    4023                                 x=1e-12; y=1e-12;
    4024                         }
    4025                         dx = x_element-x; dy = y_element-y; dz = z_element-z;
    4026                         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);
    4027          E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    4028                        
    4029                         /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    4030                         int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    4031                         U_elastic[i] += U_elastic_precomputed[index];
    4032                         N_elastic[i] += H_elastic_precomputed[index]*N_azim;
    4033                         E_elastic[i] += H_elastic_precomputed[index]*E_azim;
    4034                 }
     4010
     4011                /*Compute alpha angle between centroid and current vertex: */
     4012                lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     4013
     4014                delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     4015                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     4016
     4017                /*Compute azimuths, both north and east components: */
     4018                x = xx[i]; y = yy[i]; z = zz[i];
     4019                if(latitude[i]==90){
     4020                        x=1e-12; y=1e-12;
     4021                }
     4022                if(latitude[i]==-90){
     4023                        x=1e-12; y=1e-12;
     4024                }
     4025                dx = x_element-x; dy = y_element-y; dz = z_element-z;
     4026                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);
     4027                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     4028               
     4029                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
     4030                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
     4031                U_elastic[i] += U_elastic_precomputed[index];
     4032                N_elastic[i] += H_elastic_precomputed[index]*N_azim;
     4033                E_elastic[i] += H_elastic_precomputed[index]*E_azim;
    40354034
    40364035                /*Add all components to the pUp solution vectors:*/
    4037                 if(computerigid){
    4038                         U_values[i]+=0; N_values[i]+=0; E_values[i]+=0;
    4039                 }
    4040                 if(computeelastic){
    4041                        
    4042                         if(this->inputs->Max(MaskIceLevelsetEnum)<0){
    4043                                
    4044                                 /*Compute ice thickness change: */
    4045                                 Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    4046                                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
    4047                                 deltathickness_input->GetInputAverage(&I);
    4048                        
    4049                                 U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
    4050                                 N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
    4051                                 E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
    4052                         }
    4053                         else if(IsWaterInElement()) {
    4054                        
    4055                                 /*From Sg, recover water sea level rise:*/
    4056                                 I=0; for(int i=0;i<NUMVERTICES;i++) I+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
    4057 
    4058                                 U_values[i]+=3*rho_water/rho_earth*area/eartharea*I*U_elastic[i];
    4059                                 N_values[i]+=3*rho_water/rho_earth*area/eartharea*I*N_elastic[i];
    4060                                 E_values[i]+=3*rho_water/rho_earth*area/eartharea*I*E_elastic[i];
    4061                         }
    4062                 }
     4036                if(this->inputs->Max(MaskIceLevelsetEnum)<0){
     4037                        U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
     4038                        N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
     4039                        E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
     4040                }
     4041                else if(IsWaterInElement()) {
     4042                        U_values[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i];
     4043                        N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i];
     4044                        E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i];
     4045                }
    40634046        }
    40644047        pUp->SetValues(gsize,indices,U_values,ADD_VAL);
    … …  
    40694052        xDelete<int>(indices);
    40704053        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
    4071 
    4072         /*Free ressources:*/
    4073         if(computeelastic) {
    4074                 xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
    4075         }
     4054        xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
    40764055
    40774056        return;
Note: See TracChangeset for help on using the changeset viewer.