Changeset 26230


Ignore:
Timestamp:
05/03/21 12:34:53 (4 years ago)
Author:
Eric.Larour
Message:

CHG: speed up of the convolutions by 2/3.

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

Legend:

Unmodified
Added
Removed

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

    r26227 r26230  
    397397                virtual void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0;
    398398                virtual void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom)=0;
    399                 virtual void       SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
     399                virtual void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
     400                virtual void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom)=0;
    400401                virtual void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
    401402                virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;

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

    r26227 r26230  
    230230                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
    231231                void       SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
    232                 void       SealevelchangeConvolution(GrdLoads* loads,  Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     232                void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     233                void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
    233234                void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    234235                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};

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

    r26227 r26230  
    185185                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
    186186                void       SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
    187                 void       SealevelchangeConvolution(GrdLoads* loads,  Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     187                void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     188                void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
    188189                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    189190                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};

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

    r26227 r26230  
    192192                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
    193193                void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
    194                 void       SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     194                void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     195                void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
    195196                void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    196197                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};

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

    r26227 r26230  
    63056305        #endif
    63066306        /*}}}*/
    6307 
    6308        
    63096307        return;
    63106308
    … …  
    63516349        slgeom->isoceanin[this->lid]=isocean; //keep track for later.
    63526350        area=areae[this->sid];
     6351
     6352        /*Compute element ids, used to speed up computations in convolution phase:{{{*/
     6353        for(int i=0;i<NUMVERTICES;i++){
     6354                slgeom->lids[this->vertices[i]->lid]=this->lid;
     6355        }
     6356        /*}}}*/
    63536357
    63546358        /*set barycentre for all elements, to be updated for fractional loads in the next routine: */
    … …  
    68846888}
    68856889/*}}}*/
    6886 void       Tria::SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
     6890void       Tria::SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
    68876891
    68886892        /*sal green function:*/
    … …  
    68966900        bool sal = false;
    68976901        bool rotation= false;
     6902        bool percpu= false;
    68986903        int  size;
    68996904        int  nel,nbar;
    … …  
    69076912        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    69086913
     6914        if(sal){
     6915                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
     6916                this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
     6917                this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
     6918                this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
     6919
     6920                this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true);
     6921        }
     6922
    69096923        if(rotation){
    69106924                Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
    … …  
    69126926                Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
    69136927               
    6914                 for(int i=0;i<NUMVERTICES;i++) SealevelGRD[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
    6915         }
    6916 
    6917         if(sal){
    6918                 this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    6919                 this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
    6920                 this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
    6921                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
    6922 
    6923                 this->SealevelchangeGxL(&SealevelGRD[0],G, Gsub, loads, slgeom, nel);
    6924         }
     6928                for(int i=0;i<NUMVERTICES;i++){
     6929                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     6930                                sealevelpercpu[this->vertices[i]->lid]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
     6931                        }
     6932                }
     6933        }
     6934        return;
     6935} /*}}}*/
     6936void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
     6937
     6938        /*sal green function:*/
     6939        IssmDouble* G=NULL;
     6940        IssmDouble* Gsub[SLGEOM_NUMLOADS];
     6941        IssmDouble oceanaverage=0;
     6942        IssmDouble oceanarea=0;
     6943        IssmDouble rho_water;
     6944       
     6945        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    69256946
    69266947        /*retrieve ocean average and area:*/
    69276948        for(int i=0;i<NUMVERTICES;i++){
    6928                 oceanaverage+=SealevelGRD[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
     6949                oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
    69296950        }
    69306951        #ifdef _ISSM_DEBUG_
    … …  
    69406961        else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    69416962
     6963        /*add ocean area into a global oceanareas vector:*/
    69426964        if(!loads->sealevelloads){
    69436965                oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
    … …  
    69476969                }
    69486970        }
    6949 
    69506971       
    69516972        return;
    … …  
    69766997        bool rotation= false;
    69776998        bool elastic=false;
     6999        bool percpu=false;
    69787000
    69797001        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    … …  
    69837005        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    69847006
    6985         if(rotation){
    6986                 Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
    6987                 Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
    6988                 Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
    6989                
    6990                 for(int i=0;i<NUMVERTICES;i++) SealevelRSL[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
    6991         }
     7007       
    69927008               
    69937009        if(sal){
    … …  
    70167032                }
    70177033
    7018                 this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel);
     7034                this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false);
    70197035
    70207036                if(elastic){
    7021                         this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel);
     7037                        this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false);
    70227038                        if(horiz ){
    7023                                 this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel);
    7024                                 this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel);
     7039                                this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false);
     7040                                this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false);
     7041                        }
     7042                }
     7043        }
     7044
     7045        if(rotation){
     7046                Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
     7047                Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
     7048                Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
     7049               
     7050                for(int i=0;i<NUMVERTICES;i++){
     7051                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     7052                                SealevelRSL[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
    70257053                        }
    70267054                }
    … …  
    70407068
    70417069} /*}}}*/
    7042 void       Tria::SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel) { /*{{{*/
     7070void       Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu) { /*{{{*/
     7071
     7072        IssmDouble sealevel[3]={0};
    70437073
    70447074        if(loads->sealevelloads){
    70457075                for(int i=0;i<NUMVERTICES;i++) {
     7076                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    70467077                        for (int e=0;e<nel;e++){
    70477078                                sealevel[i]+=G[i*nel+e]*(loads->sealevelloads[e]+loads->loads[e]);
    … …  
    70627093        else{  //this is the initial convolution where only loads are provided
    70637094                for(int i=0;i<NUMVERTICES;i++) {
     7095                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    70647096                        for (int e=0;e<nel;e++){
    70657097                                sealevel[i]+=G[i*nel+e]*(loads->loads[e]);
    … …  
    70727104                        }
    70737105                }
     7106        }
     7107
     7108        /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/
     7109        if(percpu){
     7110                for(int i=0;i<NUMVERTICES;i++){
     7111                        if(slgeom->lids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i];
     7112                }
     7113        }
     7114        else{
     7115                for(int i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i];
    70747116        }
    70757117

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

    r26227 r26230  
    174174                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
    175175                void       SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom);
    176                 void       SealevelchangeConvolution(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom);
     176                void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
     177                void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom);
    177178                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
    178179                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom);
    … …  
    242243                void           UpdateConstraintsExtrudeFromBase(void);
    243244                void           UpdateConstraintsExtrudeFromTop(void);
    244                 void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel);
     245                void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool optimize);
    245246                /*}}}*/
    246247

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

    r26223 r26230  
    1717
    1818/*Object constructors and destructor*/
    19 SealevelGeometry::SealevelGeometry(int localnelin){ /*{{{*/
     19SealevelGeometry::SealevelGeometry(int localnelin,int localnodsin){ /*{{{*/
    2020        localnel=localnelin;
    2121        for(int i=0;i<SLGEOM_NUMLOADS;i++){
    … …  
    3939        longe=xNew<IssmDouble>(localnel);
    4040        isoceanin=xNew<bool>(localnel);
     41        lids=xNew<int>(localnodsin);
    4142
    4243}; /*}}}*/
    … …  
    5960        xDelete<IssmDouble>(longe);
    6061        xDelete<bool>(isoceanin);
     62        xDelete<int>(lids);
    6163}; /*}}}*/
    6264

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

    r26223 r26230  
    3535                int         nsubel[SLGEOM_NUMLOADS];
    3636                int         nbar[SLGEOM_NUMLOADS];
     37                int*        lids;
    3738               
    38                 SealevelGeometry(int localnel);
     39                SealevelGeometry(int localnel,int localnods);
    3940                ~SealevelGeometry();
    4041                void InitializeMappingsAndBarycentres(void);

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

    r26227 r26230  
    245245        int  grdmodel;
    246246        int  computesealevel=0;
     247        IssmDouble*           sealevelpercpu=NULL;
    247248
    248249        /*}}}*/
    … …  
    292293        subelementoceanareas=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
    293294        oceanareas=new Vector<IssmDouble>(nel);
     295        sealevelpercpu=xNewZeroInit<IssmDouble>(femmodel->vertices->Size());
    294296
    295297        if(VerboseSolution()) _printf0_("         starting  GRD convolutions\n");
    … …  
    322324                for(Object* & object : femmodel->elements->objects){
    323325                        Element* element = xDynamicCast<Element*>(object);
    324                         element->SealevelchangeConvolution(loads, oceanareas, subelementoceanareas, rotationaxismotionvector,slgeom);
     326                        element->SealevelchangeConvolution(sealevelpercpu, loads , rotationaxismotionvector,slgeom);
     327                }
     328               
     329                /*retrieve sea level average  and ocean area:*/
     330                for(Object* & object : femmodel->elements->objects){
     331                        Element* element = xDynamicCast<Element*>(object);
     332                        element->SealevelchangeOceanAverage(loads, oceanareas, subelementoceanareas, sealevelpercpu, slgeom);
    325333                }
    326334
    … …  
    535543        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
    536544        nel=femmodel->elements->NumberOfElements();
    537        
     545               
    538546        /*early return?:*/
    539547        if(grdmodel==IvinsEnum) return;
    … …  
    594602       
    595603        /*initialize SealevelMasks structure: */
    596         slgeom=new SealevelGeometry(femmodel->elements->Size());
     604        slgeom=new SealevelGeometry(femmodel->elements->Size(),femmodel->vertices->Size());
    597605       
    598606        /*Verbose: */
Note: See TracChangeset for help on using the changeset viewer.