Changeset 26057


Ignore:
Timestamp:
03/09/21 12:43:42 (4 years ago)
Author:
Mathieu Morlighem
Message:

CHG: fixed some valgrind stuff

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

Legend:

Unmodified
Added
Removed

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

    r26047 r26057  
    972972        xDelete<IssmDouble>(newthickness);
    973973        xDelete<IssmDouble>(cumdeltathickness);
     974        xDelete<IssmDouble>(deltathickness);
    974975        xDelete<IssmDouble>(newbase);
    975976        xDelete<IssmDouble>(newsurface);

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

    r26047 r26057  
    23102310                }
    23112311                else{
    2312                         values[i]=meltratefactor[i]*tanh((base[i]-bed[i])/thresholdthickness)*(upperdepthmelt-base[i]);
     2312                        values[i]=meltratefactor[i]*tanh((base[i]-bed[i])/thresholdthickness);//*(upperdepthmelt-base[i]);
    23132313                }
    23142314        }

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

    r26047 r26057  
    29612961                                        else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
    29622962                                                //heaviside: 0 for floating, 1 for grounded
    2963                                                 heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
     2963                                                heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(M_PI*(groundedice-calvinghaf)/haf_eps)/(2.*M_PI);
    29642964                                                calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
    29652965                                                meltingrate=heaviside*meltingrate+0.;

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

    r26056 r26057  
    15121512        /*retrieve coordinates: lat,long,radius */
    15131513        ::GetVerticesCoordinates(&llr_list[0][0],vertices,NUMVERTICES,spherical);
    1514         x1=llr_list[0][0]/180*PI; y1=llr_list[0][1]/180*PI; z1=llr_list[0][2];
    1515         x2=llr_list[1][0]/180*PI; y2=llr_list[1][1]/180*PI; z2=llr_list[1][2];
    1516         x3=llr_list[2][0]/180*PI; y3=llr_list[2][1]/180*PI; z3=llr_list[2][2];
     1514        x1=llr_list[0][0]/180.*M_PI; y1=llr_list[0][1]/180.*M_PI; z1=llr_list[0][2];
     1515        x2=llr_list[1][0]/180.*M_PI; y2=llr_list[1][1]/180.*M_PI; z2=llr_list[1][2];
     1516        x3=llr_list[2][0]/180.*M_PI; y3=llr_list[2][1]/180.*M_PI; z3=llr_list[2][2];
    15171517
    15181518        /*compute great circle distance between vertices */
    1519         arc12=2.*asin(sqrt(pow(sin((x2-x1)/2),2.0)+cos(x1)*cos(x2)*pow(sin((y2-y1)/2),2)));
    1520         arc23=2.*asin(sqrt(pow(sin((x3-x2)/2),2.0)+cos(x2)*cos(x3)*pow(sin((y3-y2)/2),2)));
    1521         arc31=2.*asin(sqrt(pow(sin((x1-x3)/2),2.0)+cos(x3)*cos(x1)*pow(sin((y1-y3)/2),2)));
     1519        arc12=2.*asin(sqrt(pow(sin(0.5*(x2-x1)),2)+cos(x1)*cos(x2)*pow(sin(0.5*(y2-y1)),2)));
     1520        arc23=2.*asin(sqrt(pow(sin(0.5*(x3-x2)),2)+cos(x2)*cos(x3)*pow(sin(0.5*(y3-y2)),2)));
     1521        arc31=2.*asin(sqrt(pow(sin(0.5*(x1-x3)),2)+cos(x3)*cos(x1)*pow(sin(0.5*(y1-y3)),2)));
    15221522
    15231523        /*semi parameter */
    … …  
    21192119
    21202120        /*Assign output pointers*/
     2121        xDelete<IssmDouble>(timesteps);
    21212122        *pvalues=values;
    21222123        *ptimes=times;
    … …  
    36413642                                        else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
    36423643                                                //heaviside: 0 for floating, 1 for grounded
    3643                                                 heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
     3644                                                heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*M_PI);
    36443645                                                calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
    36453646                                                meltingrate=heaviside*meltingrate+0.;
    … …  
    51735174                dx = x_element - xx[i];         dy = y_element - yy[i];
    51745175                dist = sqrt(pow(dx,2)+pow(dy,2));                                               // distance between vertex and elemental centroid [m]
    5175                 alpha = dist*360.0/(2*PI*earth_radius) * PI/180.0;      // [in radians] 360 degree = 2*pi*earth_radius
     5176                alpha = dist*360.0/(2*M_PI*earth_radius) * M_PI/180.0;  // [in radians] 360 degree = 2*pi*earth_radius
    51765177
    51775178                /*Compute azimuths, both north and east components: */
    5178                 ang = PI/2 - atan2(dy,dx);              // this is bearing angle!
     5179                ang = M_PI/2 - atan2(dy,dx);            // this is bearing angle!
    51795180                Y_azim = cos(ang);
    51805181                X_azim = sin(ang);
    51815182
    51825183                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    5183                 int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
     5184                int index=reCast<int,IssmDouble>(alpha/M_PI*(M-1));
    51845185                U_elastic[i] += U_elastic_precomputed[index];
    51855186                Y_elastic[i] += H_elastic_precomputed[index]*Y_azim;
    … …  
    51875188
    51885189                /*Add all components to the pUp solution vectors:*/
    5189                 U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i];
    5190                 Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i];
    5191                 X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i];
     5190                U_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*U_elastic[i];
     5191                Y_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*Y_elastic[i];
     5192                X_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*X_elastic[i];
    51925193
    51935194                /*North-south, East-west components */
    51945195                if (hemi == -1) {
    5195                         ang2 = PI/2 - atan2(yy[i],xx[i]);
     5196                        ang2 = M_PI/2 - atan2(yy[i],xx[i]);
    51965197                }
    51975198                else if (hemi == 1) {
    5198                         ang2 = PI/2 - atan2(-yy[i],-xx[i]);
     5199                        ang2 = M_PI/2 - atan2(-yy[i],-xx[i]);
    51995200                }
    52005201                if (hemi != 0){
    … …  
    52035204                        N_elastic[i] += H_elastic_precomputed[index]*N_azim;
    52045205                        E_elastic[i] += H_elastic_precomputed[index]*E_azim;
    5205                         N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i];
    5206                         E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_elastic[i];
     5206                        N_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*N_elastic[i];
     5207                        E_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*E_elastic[i];
    52075208                }
    52085209        }
    … …  
    53085309        if(longe>180)longe=longe-360;
    53095310
    5310         late=late/180*PI;
    5311         longe=longe/180*PI;
     5311        late=late/180.*M_PI;
     5312        longe=longe/180.*M_PI;
    53125313        /*}}}*/
    53135314
    … …  
    53445345
    53455346                /*Compute alpha angle between centroid and current vertex: */
    5346                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     5347                lati=latitude[i]/180.*M_PI; longi=longitude[i]/180.*M_PI;
    53475348
    53485349                delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    … …  
    53625363
    53635364                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    5364                 int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
     5365                int index=reCast<int,IssmDouble>(alpha/M_PI*(M-1));
    53655366                U_elastic[i] += U_elastic_precomputed[index];
    53665367                N_elastic[i] += H_elastic_precomputed[index]*N_azim;
    … …  
    54555456        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
    54565457
    5457         late=90-late;
    5458         if(longe>180)longe=(longe-180)-180;
    5459 
    5460         late=late/180*PI;
    5461         longe=longe/180*PI;
     5458        late=90.-late;
     5459        if(longe>180.)longe=(longe-180.)-180.;
     5460
     5461        late=late/180.*M_PI;
     5462        longe=longe/180.*M_PI;
    54625463        /*}}}*/
    54635464        re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
    … …  
    54775478                 * ALL in geographic coordinates
    54785479                 * */
    5479                 dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
    5480                 dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
    5481                 dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
     5480                dI_list[0] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
     5481                dI_list[1] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
     5482                dI_list[2] = +4*M_PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
    54825483        }
    54835484        else if(masks->isiceonly[this->lid]){
    … …  
    54925493                deltathickness_input->GetInputAverage(&I);
    54935494
    5494                 dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
    5495                 dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
    5496                 dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
     5495                dI_list[0] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
     5496                dI_list[1] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
     5497                dI_list[2] = +4*M_PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
    54975498        }
    54985499
    … …  
    56215622
    56225623                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    5623                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    5624                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
     5624                lati=latitude[i]/180.*M_PI; longi=longitude[i]/180.*M_PI;
     5625                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    56255626                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    5626                 indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
     5627                indices[i]=alpha/M_PI*reCast<IssmDouble,int>(M-1);
    56275628                index=reCast<int,IssmDouble>(indices[i]);
    56285629
    … …  
    61456146
    61466147        /*element radius: */
    6147         IssmDouble re=sqrt(area/PI);
     6148        IssmDouble re=sqrt(area/M_PI);
    61486149
    61496150        /*figure out gravity center of our element: */

  • issm/trunk-jpl/src/c/main/esmfbinders.cpp

    r25827 r26057  
    153153        } /*}}}*/
    154154
     155        /*TODO: we need 2 initialize routines, the second one will be empty for now
     156         * In: ESMF config, ESMF Field bundle
     157         */
     158
    155159        /*FISOC*/
    156160        void InitializeISSM_FISOC(int argc, char** argv,ISSM_MPI_Comm comm_init){ /*{{{*/
    … …  
    165169                /*Initialize ESMC Mesh*/
    166170                int pdim;        /*parametric dimension is the same as the domain dimensions */
    167                 int sdim = 3;    /*coordinates of each vertex is always 2  (just x,y for now) */
     171                int sdim = 2;    /*coordinates of each vertex is always 2  (just x,y for now) */
    168172                ESMC_CoordSys_Flag coordsys = ESMC_COORDSYS_CART; /*Cartesian coordinate system by default */
    169173                femmodel->parameters->FindParam(&pdim,DomainDimensionEnum);
    … …  
    171175                if(rc!=ESMF_SUCCESS) _error_("could not create EMSC_Mesh");
    172176
    173                 /*How to install ESMG with lib/libesmf.so
    174                  * do indices need to be 1 based?
    175                  * ESMF installation, how can we get everything in /lib
     177                /*
    176178                 * What do we do with vertices at the boundary, declare twice?
    177                  * Do we need a restart file or save femmodel somewhere?*/
     179                 * */
    178180
    179181                /*Add nodes (which are ISSM Vertices)*/
    … …  
    182184                int        *nodeOwner = xNew<int>(numnodes);
    183185                IssmDouble *nodeCoord = xNew<IssmDouble>(sdim*numnodes);
    184                 for (int i=0;i<femmodel->elements->Size();i++){
     186                for (int i=0;i<femmodel->vertices->Size();i++){
    185187                        Vertex* vertex = xDynamicCast<Vertex*>(femmodel->vertices->GetObjectByOffset(i));
    186188                        nodeId[i]           = vertex->Sid()+1;
    … …  
    213215                xDelete<int>(elemConn);
    214216
    215                 /*DO we need to create fields here? https://earthsystemmodeling.org/docs/nightly/develop/ESMC_crefdoc/node5.html#SECTION05024400000000000000*/
    216                 //ESMC_InterArrayInt *gridToFieldMap,       // in
    217                 //ESMC_InterArrayInt *ungriddedLBound,      // in
    218                 //ESMC_InterArrayInt *ungriddedUBound,      // in
    219                 //ESMC_Field esmf_shelfbase = ESMC_FieldCreateMeshTypeKind(mesh,ESMC_TYPEKIND_R8,???,gridToFieldMap,ungriddedLBound,ungriddedUBound,"ShelfTopography",&rc);
    220                 if(rc!=ESMF_SUCCESS) _error_("could not create EMSC_Field");
    221 
    222217                /*Create restart file for later */
    223218                femmodel->Restart();
Note: See TracChangeset for help on using the changeset viewer.