Changeset 24915


Ignore:
Timestamp:
05/28/20 21:38:38 (5 years ago)
Author:
Eric.Larour
Message:

CHG: 100 fold increase in speed of the sea level solver.

Location:
issm/trunk-jpl/src/c/classes
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r24889 r24915  
    377377                virtual IssmDouble    OceanArea(void)=0;
    378378                virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0;
    379                 virtual void          SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
    380                 virtual void          SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
    381                 virtual void          SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0;
     379                virtual void          SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
     380                virtual void          SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
     381                virtual void          SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0;
    382382                #endif
    383383

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

    r24790 r24915  
    215215                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
    216216                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    217                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    218                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    219                 void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     217                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     218                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
     219                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
    220220                #endif
    221221

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

    r24790 r24915  
    173173#ifdef _HAVE_SEALEVELRISE_
    174174                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    175                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    176                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    177                 void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     175                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     176                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
     177                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
    178178                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    179179                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};

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

    r24790 r24915  
    179179#ifdef _HAVE_SEALEVELRISE_
    180180                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    181                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    182                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    183                 void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     181                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     182                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
     183                void    SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
    184184                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    185185                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};

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

    r24889 r24915  
    55645564        return;
    55655565}/*}}}*/
    5566 void    Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5566void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    55675567
    55685568        /*Computational flags:*/
    … …  
    55755575                /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
    55765576                 *bottom pressure fingerprints:*/
    5577                 if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
     5577                if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
    55785578        }
    55795579        //if(!IsIceInElement()){
    55805580                /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
    5581                 this->SealevelriseEustaticIce(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
     5581                this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
    55825582        //}
    55835583
    55845584}
    55855585/*}}}*/
    5586 void    Tria::SealevelriseEustaticIce(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5586void    Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    55875587
    55885588        /*diverse:*/
    … …  
    57445744
    57455745        if(computeelastic | computerigid){
    5746                 int* indices=xNew<int>(gsize);
    5747                 IssmDouble* values=xNew<IssmDouble>(gsize);
    57485746                IssmDouble alpha;
    57495747                IssmDouble delPhi,delLambda;
    57505748                for(int i=0;i<gsize;i++){
    5751                         indices[i]=i;
    57525749
    57535750                        IssmDouble G_rigid=0;  //do not remove =0!
    … …  
    57575754                        lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    57585755
    5759                    delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     5756                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    57605757                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    57615758
    … …  
    57695766                        }
    57705767
    5771                         /*Add all components to the pSgi or pSgo solution vectors:*/
    5772                         values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
    5773                 }
    5774                 pSgi->SetValues(gsize,indices,values,ADD_VAL);
    5775 
    5776                 /*free ressources:*/
    5777                 xDelete<IssmDouble>(values);
    5778                 xDelete<int>(indices);
     5768                        /*Add all components to the Sgi or Sgo solution vectors:*/
     5769                        Sgi[i]+=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
     5770
     5771                }
    57795772        }
    57805773
    … …  
    57865779}
    57875780/*}}}*/
    5788 void    Tria::SealevelriseEustaticBottomPressure(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5781void    Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    57895782
    57905783        /*diverse:*/
    … …  
    58985891
    58995892        if(computeelastic | computerigid){
    5900                 int* indices=xNew<int>(gsize);
    5901                 IssmDouble* values=xNew<IssmDouble>(gsize);
    59025893                IssmDouble alpha;
    59035894                IssmDouble delPhi,delLambda;
    59045895                for(int i=0;i<gsize;i++){
    5905                         indices[i]=i;
    59065896
    59075897                        IssmDouble G_rigid=0;  //do not remove =0!
    … …  
    59245914
    59255915                        /*Add all components to the pSgi or pSgo solution vectors:*/
    5926                         values[i]=3*rho_water/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
    5927                 }
    5928                 pSgi->SetValues(gsize,indices,values,ADD_VAL);
    5929 
    5930                 /*free ressources:*/
    5931                 xDelete<IssmDouble>(values);
    5932                 xDelete<int>(indices);
     5916                        Sgi[i]+=3*rho_water/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
     5917                }
    59335918        }
    59345919
    … …  
    59405925}
    59415926/*}}}*/
    5942 void    Tria::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/
     5927void    Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/
    59435928
    59445929        /*diverse:*/
    … …  
    60466031        if(computerigid) G_rigid=xNewZeroInit<IssmDouble>(gsize);
    60476032
    6048         int* indices=xNew<int>(gsize);
    6049         IssmDouble* values=xNewZeroInit<IssmDouble>(gsize);
    60506033        IssmDouble alpha;
    60516034        IssmDouble delPhi,delLambda;
    60526035
    60536036        for(int i=0;i<gsize;i++){
    6054 
    6055                 indices[i]=i;
    60566037
    60576038                /*Compute alpha angle between centroid and current vertex : */
    … …  
    60646045                if(computerigid){
    60656046                        G_rigid[i]=1.0/2.0/sin(alpha/2.0);
    6066                         values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i];
     6047                        //values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i];
     6048                        Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i];
    60676049                }
    60686050
    … …  
    60716053                        int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    60726054                        G_elastic[i] += G_elastic_precomputed[index];
    6073                         values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
    6074                 }
    6075         }
    6076 
    6077         pSgo->SetValues(gsize,indices,values,ADD_VAL);
    6078 
    6079         /*free ressources:*/
    6080         xDelete<IssmDouble>(values);
    6081         xDelete<int>(indices);
     6055                        //values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
     6056                        Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
     6057                }
     6058        }
     6059
    60826060
    60836061        /*Free ressources:*/
    … …  
    60886066}
    60896067/*}}}*/
    6090 void    Tria::SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/
     6068void    Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/
    60916069
    60926070        /*diverse:*/
    … …  
    62176195        }
    62186196
    6219         int* indices=xNew<int>(gsize);
    6220         U_values=xNewZeroInit<IssmDouble>(gsize);
     6197        //int* indices=xNew<int>(gsize);
     6198        //U_values=xNewZeroInit<IssmDouble>(gsize);
    62216199        if(horiz){
    6222                 N_values=xNewZeroInit<IssmDouble>(gsize);
    6223                 E_values=xNewZeroInit<IssmDouble>(gsize);
     6200                //N_values=xNewZeroInit<IssmDouble>(gsize);
     6201                //E_values=xNewZeroInit<IssmDouble>(gsize);
    62246202        }
    62256203        IssmDouble alpha;
    … …  
    62306208        for(int i=0;i<gsize;i++){
    62316209
    6232                 indices[i]=i;
     6210                //indices[i]=i;
    62336211
    62346212                /*Compute alpha angle between centroid and current vertex: */
    … …  
    62626240                /*Add all components to the pUp solution vectors:*/
    62636241                if(IsIceOnlyInElement()){
    6264                         U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
     6242                        Up[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
    62656243                        if(horiz){
    6266                                 N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
    6267                                 E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
     6244                                North[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
     6245                                East[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
    62686246                        }
    62696247                }
    62706248                else if(IsOceanInElement()) {
    6271                         U_values[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i];
     6249                        Up[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i];
    62726250                        if(horiz){
    6273                                 N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i];
    6274                                 E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i];
     6251                                North[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i];
     6252                                East[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i];
    62756253                        }
    62766254                }
    62776255        }
    6278         pUp->SetValues(gsize,indices,U_values,ADD_VAL);
    6279         if(horiz){
    6280                 pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
    6281                 pEast->SetValues(gsize,indices,E_values,ADD_VAL);
    6282         }
    62836256
    62846257        /*free ressources:*/
    6285         xDelete<int>(indices);
    62866258        xDelete<IssmDouble>(U_values);
    62876259        xDelete<IssmDouble>(U_elastic);

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

    r24790 r24915  
    166166                IssmDouble OceanAverage(IssmDouble* Sg);
    167167                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea);
    168                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    169                 void    SealevelriseEustaticIce(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    170                 void    SealevelriseEustaticBottomPressure(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    171                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
    172                 void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz);
     168                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     169                void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     170                void    SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     171                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
     172                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz);
    173173                #endif
    174174                /*}}}*/

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

    r24861 r24915  
    46944694        IssmDouble  eartharea      = 0.;
    46954695        IssmDouble  eartharea_cpu  = 0.;
    4696         int         ns,nsmax;
    4697 
    4698         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    4699         ns = elements->Size();
     4696        IssmDouble* RSLgi  = NULL;
     4697        int* indices = NULL;
     4698        int         gsize;
     4699
     4700        /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior
     4701         * to assembly:*/
     4702        gsize = this->nodes->NumberOfDofs(GsetEnum);
     4703        RSLgi=xNewZeroInit<IssmDouble>(gsize);
     4704        indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
    47004705
    47014706        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
    4702         for(int i=0;i<ns;i++){
     4707        for(int i=0;i<elements->Size();i++){
    47034708                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    47044709                oceanarea_cpu += element->OceanArea();
    … …  
    47124717        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    47134718
    4714         /*Figure out max of ns: */
    4715         ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    4716         ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    4717 
    47184719        /*Call the sea level rise core: */
    4719         for(int i=0;i<nsmax;i++){
    4720                 if(i<ns){
    4721 
    4722                         if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%  ");
    4723 
    4724                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4725                         element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
    4726                         eustatic_cpu+=eustatic_cpu_e;
    4727                 }
    4728                 if(i%loop==0)pRSLgi->Assemble();
    4729         }
    4730         if(VerboseConvergence())_printf0_("\n");
    4731 
    4732         /*One last time: */
     4720        for(int i=0;i<elements->Size();i++){
     4721                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     4722                element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
     4723                eustatic_cpu+=eustatic_cpu_e;
     4724        }
     4725        /*Plug values once and assemble: */
     4726        pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL);
    47334727        pRSLgi->Assemble();
    4734 
     4728       
    47354729        /*Sum all eustatic components from all cpus:*/
    47364730        ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    … …  
    47384732        _assert_(!xIsNan<IssmDouble>(eustatic));
    47394733
     4734        /*Free ressources:*/
     4735        xDelete<int>(indices);
     4736        xDelete<IssmDouble>(RSLgi);
     4737
    47404738        /*Assign output pointers:*/
    47414739        *peustatic=eustatic;
    … …  
    47514749        IssmDouble  eartharea_cpu=0;
    47524750
    4753         int         ns,nsmax;
     4751        IssmDouble* RSLgo  = NULL;
     4752        int* indices = NULL;
     4753        int         gsize;
     4754
     4755        /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior
     4756         * to assembly:*/
     4757        gsize = this->nodes->NumberOfDofs(GsetEnum);
     4758        RSLgo=xNewZeroInit<IssmDouble>(gsize);
     4759        indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
    47544760
    47554761        /*Serialize vectors from previous iteration:*/
    47564762        RSLg_old=pRSLg_old->ToMPISerial();
    47574763
    4758         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    4759         ns = elements->Size();
    4760 
    47614764        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
    4762         for(int i=0;i<ns;i++){
     4765        for(int i=0;i<elements->Size();i++){
    47634766                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    47644767                eartharea_cpu += element->GetAreaSpherical();
    … …  
    47684771        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    47694772
    4770         /*Figure out max of ns: */
    4771         ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    4772         ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    4773 
    47744773        /*Call the sea level rise core: */
    4775         for(int i=0;i<nsmax;i++){
    4776                 if(i<ns){
    4777                         if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%   ");
    4778                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4779                         element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea);
    4780                 }
    4781                 if(i%loop==0)pRSLgo->Assemble();
    4782         }
    4783         if(verboseconvolution)if(VerboseConvergence())_printf0_("\n");
     4774        for(int i=0;i<elements->Size();i++){
     4775                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     4776                element->SealevelriseNonEustatic(RSLgo,RSLg_old,latitude,longitude,radius,eartharea);
     4777        }
     4778        pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
     4779        pRSLgo->Assemble();
    47844780
    47854781        /*Free ressources:*/
     4782        xDelete<int>(indices);
     4783        xDelete<IssmDouble>(RSLgo);
    47864784        xDelete<IssmDouble>(RSLg_old);
     4785
    47874786}
    47884787/*}}}*/
    … …  
    48814880        IssmDouble  eartharea_cpu=0;
    48824881
    4883         int         ns,nsmax;
     4882        IssmDouble* Up  = NULL;
     4883        IssmDouble* North  = NULL;
     4884        IssmDouble* East  = NULL;
     4885        int* indices = NULL;
     4886        int         gsize;
     4887
    48844888
    48854889        /*Serialize vectors from previous iteration:*/
    48864890        RSLg=pRSLg->ToMPISerial();
    48874891
    4888         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    4889         ns = elements->Size();
     4892        /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior
     4893         * to assembly:*/
     4894        gsize = this->nodes->NumberOfDofs(GsetEnum);
     4895        Up=xNewZeroInit<IssmDouble>(gsize);
     4896        North=xNewZeroInit<IssmDouble>(gsize);
     4897        East=xNewZeroInit<IssmDouble>(gsize);
     4898        indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
    48904899
    48914900        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
    4892         for(int i=0;i<ns;i++){
     4901        for(int i=0;i<elements->Size();i++){
    48934902                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    48944903                eartharea_cpu += element->GetAreaSpherical();
    … …  
    48974906        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    48984907
    4899         /*Figure out max of ns: */
    4900         ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    4901         ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    4902 
    49034908        /*Call the sea level rise core: */
    4904         for(int i=0;i<nsmax;i++){
    4905                 if(i<ns){
    4906                         if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%  ");
    4907                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4908                         element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
    4909                 }
    4910                 if(i%loop==0){
    4911                         pUp->Assemble();
    4912                         if (horiz){
    4913                                 pNorth->Assemble();
    4914                                 pEast->Assemble();
    4915                         }
    4916                 }
    4917         }
    4918 
    4919         /*One last time: */
     4909        for(int i=0;i<elements->Size();i++){
     4910                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     4911                element->SealevelriseGeodetic(Up,North,East,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
     4912        }
     4913
     4914        pUp->SetValues(gsize,indices,Up,ADD_VAL);
    49204915        pUp->Assemble();
    4921         if(horiz){
     4916        if (horiz){
     4917                pNorth->SetValues(gsize,indices,North,ADD_VAL);
    49224918                pNorth->Assemble();
     4919                pEast->SetValues(gsize,indices,East,ADD_VAL);
    49234920                pEast->Assemble();
    49244921        }
    4925         if(VerboseConvergence())_printf0_("\n");
    49264922
    49274923        /*Free ressources:*/
     4924        xDelete<IssmDouble>(Up);
     4925        xDelete<IssmDouble>(North);
     4926        xDelete<IssmDouble>(East);
     4927        xDelete<int>(indices);
    49284928        xDelete<IssmDouble>(RSLg);
    49294929}
Note: See TracChangeset for help on using the changeset viewer.