Changeset 23602


Ignore:
Timestamp:
01/05/19 22:36:12 (6 years ago)
Author:
Mathieu Morlighem
Message:

NEW: Mergesolution from f to g is now working without MPI calls, for speed and efficiency, and cleaned up many things in the code'

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r23586 r23602  
    239239                                        ./modules/ResetFSBasalBoundaryConditionx/ResetFSBasalBoundaryConditionx.cpp\
    240240                                        ./modules/Solverx/Solverx.cpp\
    241                                         ./modules/VecMergex/VecMergex.cpp\
    242241                                        ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
    243242                                        ./cores/ProcessArguments.cpp\
    … …  
    352351                                        ./toolkits/petsc/patches/VecToMPISerial.cpp\
    353352                                        ./toolkits/petsc/patches/MatToSerial.cpp\
    354                                         ./toolkits/petsc/patches/VecMerge.cpp\
    355353                                        ./toolkits/petsc/patches/NewVec.cpp\
    356354                                        ./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\

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

    r23600 r23602  
    590590}
    591591/*}}}*/
    592 void Node::VecMerge(Vector<IssmDouble>* ug, IssmDouble* vector_serial,int setenum){/*{{{*/
    593 
    594         IssmDouble *values  = NULL;
    595         int        *indices = NULL;
    596         int         count   = 0;
    597         int         i;
    598 
    599         if(setenum==FsetEnum){
    600                 if(this->indexing.fsize){
    601                         indices=xNew<int>(this->indexing.fsize);
    602                         values=xNew<IssmDouble>(this->indexing.fsize);
    603 
    604                         for(i=0;i<this->indexing.gsize;i++){
    605                                 if(this->indexing.f_set[i]){
    606                                         _assert_(vector_serial);
    607                                         values[count]=vector_serial[this->indexing.fdoflist[count]];
    608                                         indices[count]=this->indexing.gdoflist[i];
    609                                         count++;
    610                                 }
    611                         }
    612 
    613                         /*Add values into ug: */
    614                         ug->SetValues(this->indexing.fsize,indices,values,INS_VAL);
    615                 }
    616         }
    617         else if(setenum==SsetEnum){
    618                 if(this->indexing.ssize){
    619                         indices=xNew<int>(this->indexing.ssize);
    620                         values=xNew<IssmDouble>(this->indexing.ssize);
    621 
    622                         for(i=0;i<this->indexing.gsize;i++){
    623                                 if(this->indexing.s_set[i]){
    624                                         _assert_(vector_serial);
    625                                         values[count]=vector_serial[this->indexing.sdoflist[count]];
    626                                         indices[count]=this->indexing.gdoflist[i];
    627                                         count++;
    628                                 }
    629                         }
    630 
    631                         /*Add values into ug: */
    632                         ug->SetValues(this->indexing.ssize,indices,values,INS_VAL);
    633                 }
    634         }
    635         else _error_("VecMerge can only merge from the s or f-set onto the g-set!");
    636 
    637         /*Free ressources:*/
    638         xDelete<IssmDouble>(values);
    639         xDelete<int>(indices);
     592void Node::VecMerge(Vector<IssmDouble>* ug,IssmDouble* local_uf,int* indices_uf,int* pindex_uf,IssmDouble* local_ys,int* indices_ys,int* pindex_ys){/*{{{*/
     593
     594        /*Only perform operation if not clone*/
     595        if(this->IsClone()) return;
     596
     597        /*Recover indices*/
     598        int ind_uf = *pindex_uf;
     599        int ind_ys = *pindex_ys;
     600
     601        if(this->indexing.fsize){
     602                int*        indices = xNew<int>(this->indexing.fsize);
     603                IssmDouble* values  = xNew<IssmDouble>(this->indexing.fsize);
     604
     605                int count = 0;
     606                for(int i=0;i<this->indexing.gsize;i++){
     607                        if(this->indexing.f_set[i]){
     608                                _assert_(local_uf);
     609                                _assert_(this->indexing.fdoflist[count]==indices_uf[ind_uf]);
     610
     611                                values[count]=local_uf[ind_uf];
     612                                indices[count]=this->indexing.gdoflist[i];
     613                                count++;
     614                                ind_uf++;
     615                        }
     616                }
     617                ug->SetValues(this->indexing.fsize,indices,values,INS_VAL);
     618                /*Free ressources:*/
     619                xDelete<IssmDouble>(values);
     620                xDelete<int>(indices);
     621        }
     622        if(this->indexing.ssize){
     623                int*        indices = xNew<int>(this->indexing.ssize);
     624                IssmDouble* values  = xNew<IssmDouble>(this->indexing.ssize);
     625
     626                int count = 0;
     627                for(int i=0;i<this->indexing.gsize;i++){
     628                        if(this->indexing.s_set[i]){
     629                                _assert_(local_ys);
     630                                _assert_(this->indexing.sdoflist[count]==indices_ys[ind_ys]);
     631
     632                                values[count]=local_ys[ind_ys];
     633                                indices[count]=this->indexing.gdoflist[i];
     634                                count++;
     635                                ind_ys++;
     636                        }
     637                }
     638                ug->SetValues(this->indexing.ssize,indices,values,INS_VAL);
     639                /*Free ressources:*/
     640                xDelete<IssmDouble>(values);
     641                xDelete<int>(indices);
     642        }
     643
     644        /*Update index values*/
     645        *pindex_uf = ind_uf;
     646        *pindex_ys = ind_ys;
    640647}
    641648/*}}}*/

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

    r23600 r23602  
    8181                int   Sid(void);
    8282                void  UpdateCloneDofs(int* alltruerows,int setenum);
    83                 void  VecMerge(Vector<IssmDouble>* ug, IssmDouble* vector_serial,int setenum);
     83                void  VecMerge(Vector<IssmDouble>* ug,IssmDouble* local_uf,int* indices_uf,int* pindex_uf,IssmDouble* local_ys,int* indices_ys,int* pindex_ys);
    8484                void  VecReduce(Vector<IssmDouble>* vector, IssmDouble* ug_serial,int setnum);
    8585                void  SetApproximation(int in_approximation);

  • issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp

    r23587 r23602  
    33 */
    44
    5 #include "../VecMergex/VecMergex.h"
    65#include "../../shared/io/io.h"
    76#include "./Mergesolutionfromftogx.h"
    … …  
    98void    Mergesolutionfromftogx( Vector<IssmDouble>** pug, Vector<IssmDouble>* uf, Vector<IssmDouble>* ys, Nodes* nodes, Parameters* parameters, bool flag_ys0){
    109
    11         /*output: */
    12         Vector<IssmDouble>* ug=NULL;
    13 
    14         /*intermediary: */
    15         int gsize,fsize,ssize;
    16 
    1710        /*Display message*/
    1811        if(VerboseModule()) _printf0_("   Merging solution vector from fset to gset\n");
    1912
    2013        /*first, get gsize, fsize and ssize: */
    21         gsize=nodes->NumberOfDofs(GsetEnum);
    22         fsize=nodes->NumberOfDofs(FsetEnum);
    23         ssize=nodes->NumberOfDofs(SsetEnum);
     14        int gsize=nodes->NumberOfDofs(GsetEnum);
     15        int fsize=nodes->NumberOfDofs(FsetEnum);
     16        int ssize=nodes->NumberOfDofs(SsetEnum);
    2417
    2518        /*serialize uf and ys: those two vectors will be indexed by the nodes, who are the only ones
    … …  
    3124        }
    3225
     26        /*Get local vectors ys and uf*/
     27        int        *indices_ys = NULL;
     28        IssmDouble *local_ys   = NULL;
     29        ys->GetLocalVector(&local_ys,&indices_ys);
     30        int        *indices_uf = NULL;
     31        IssmDouble *local_uf   = NULL;
     32        uf->GetLocalVector(&local_uf,&indices_uf);
     33
    3334        /*initialize ug: */
    34         ug=new Vector<IssmDouble>(gsize);
     35        Vector<IssmDouble>* ug=new Vector<IssmDouble>(gsize);
    3536
    36         /*Merge f set back into g set: */
    37         if(fsize){
    38                 VecMergex(ug,uf,nodes,parameters,FsetEnum);
     37        /*Let nodes figure it out*/
     38        int index_uf = 0;
     39        int index_ys = 0;
     40        for(int i=0;i<nodes->Size();i++){
     41                Node* node=(Node*)nodes->GetObjectByOffset(i);
     42                node->VecMerge(ug,local_uf,indices_uf,&index_uf,local_ys,indices_ys,&index_ys);
    3943        }
    4044
    41         /*Merge s set back into g set: */
    42         if(ssize){
    43                 VecMergex(ug,ys,nodes,parameters,SsetEnum);
    44         }
     45        /*Assemble vector: */
     46        ug->Assemble();
    4547
    46         /*Assign correct pointer*/
     48        /*Cleanup and assounf output pointer*/
     49        xDelete<int>(indices_uf);
     50        xDelete<int>(indices_ys);
     51        xDelete<IssmDouble>(local_uf);
     52        xDelete<IssmDouble>(local_ys);
    4753        *pug=ug;
    4854}

  • issm/trunk-jpl/src/c/modules/modules.h

    r23586 r23602  
    9999#include "./UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h"
    100100#include "./VertexCoordinatesx/VertexCoordinatesx.h"
    101 #include "./VecMergex/VecMergex.h"
    102101#endif

  • issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h

    r19977 r23602  
    3737                virtual void GetSize(int* pM)=0;
    3838                virtual void GetLocalSize(int* pM)=0;
     39                virtual void GetLocalVector(double** pvector,int** pindices)=0;
    3940                virtual IssmAbsVec<doubletype>* Duplicate(void)=0;
    4041                virtual void Set(doubletype value)=0;

  • issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h

    r22738 r23602  
    337337                }
    338338                /*}}}*/
     339                void GetLocalVector(double** pvector,int** pindices){/*{{{*/
     340
     341                        _assert_(this->vector);
     342
     343                        /*First, check that vector size is not 0*/
     344                        int vector_size;
     345                        this->GetSize(&vector_size);
     346                        if(vector_size==0){
     347                                *pvector=NULL;
     348                                *pindices=NULL;
     349                                return;
     350                        }
     351
     352                        /*Get Ownership range*/
     353                        int lower_row,upper_row;
     354                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
     355                        int range=upper_row-lower_row;   
     356
     357                        /*return NULL if no range*/
     358                        if(range==0){
     359                                *pvector=NULL;
     360                                *pindices=NULL;
     361                                return;
     362                        }
     363
     364                        /*Build indices*/
     365                        int* indices=xNew<int>(range);
     366                        for(int i=0;i<range;i++) indices[i]=lower_row+i;
     367
     368                        /*Get vector*/
     369                        _assert_(range==this->m);
     370                        IssmDouble* values =xNew<double>(range);
     371                        xMemCpy<doubletype>(values,this->vector,this->m);
     372
     373                        *pvector  = values;
     374                        *pindices = indices;
     375                } /*}}}*/
    339376                IssmMpiVec<doubletype>* Duplicate(void){/*{{{*/
    340377

  • issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h

    r22737 r23602  
    141141                /*}}}*/
    142142                void GetLocalSize(int* pM){/*{{{*/
    143 
    144143                        *pM=this->M;
    145 
    146                 }
    147                 /*}}}*/
     144                }/*}}}*/
     145                void GetLocalVector(double** pvector,int** pindices){/*{{{*/
     146                        _error_("not implemented");
     147                } /*}}}*/
    148148                IssmSeqVec<doubletype>* Duplicate(void){/*{{{*/
    149149

  • issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h

    r19977 r23602  
    153153                }
    154154                /*}}}*/
     155                void GetLocalVector(double** pvector,int** pindices){/*{{{*/
     156                        vector->GetLocalVector(pvector,pindices);
     157                } /*}}}*/
    155158                IssmVec<doubletype>* Duplicate(void){/*{{{*/
    156159

  • issm/trunk-jpl/src/c/toolkits/objects/Vector.h

    r22558 r23602  
    217217                }
    218218                /*}}}*/
     219                void GetLocalVector(double** pvector,int** pindices){_assert_(this);/*{{{*/
     220
     221                        if(type==PetscVecType){
     222                                #ifdef _HAVE_PETSC_
     223                                this->pvector->GetLocalVector(pvector,pindices);
     224                                #endif
     225                        }
     226                        else this->ivector->GetLocalVector(pvector,pindices);
     227
     228                }
     229                /*}}}*/
    219230                Vector<doubletype>* Duplicate(void){_assert_(this);/*{{{*/
    220231

  • issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp

    r23255 r23602  
    113113        _assert_(this->vector);
    114114        VecGetSize(this->vector,pM);
    115 
    116 }
    117 /*}}}*/
    118 void PetscVec::GetLocalSize(int* pM){/*{{{*/
    119 
    120         _assert_(this->vector);
    121         VecGetLocalSize(this->vector,pM);
    122 
    123 }
    124 /*}}}*/
     115}
     116/*}}}*/
     117void PetscVec::GetLocalVector(double** pvector,int** pindices){/*{{{*/
     118
     119        _assert_(this->vector);
     120
     121        /*First, check that vector size is not 0*/
     122        int vector_size;
     123        this->GetSize(&vector_size);
     124        if(vector_size==0){
     125                *pvector=NULL;
     126                *pindices=NULL;
     127                return;
     128        }
     129
     130        /*Get Ownership range*/
     131        PetscInt lower_row,upper_row;
     132        VecGetOwnershipRange(this->vector,&lower_row,&upper_row);
     133        int range=upper_row-lower_row;   
     134
     135        /*return NULL if no range*/
     136        if(range==0){
     137                *pvector=NULL;
     138                *pindices=NULL;
     139                return;
     140        }
     141
     142        /*Build indices*/
     143        int* indices=xNew<int>(range);
     144        for(int i=0;i<range;i++) indices[i]=lower_row+i;
     145        /*Get vector*/
     146        IssmDouble* values =xNew<double>(range);
     147        VecGetValues(this->vector,range,indices,values);
     148
     149        *pvector  = values;
     150        *pindices = indices;
     151} /*}}}*/
    125152PetscVec* PetscVec::Duplicate(void){/*{{{*/
    126153

  • issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h

    r23255 r23602  
    4545                void        GetSize(int* pM);
    4646                void        GetLocalSize(int* pM);
     47                void        GetLocalVector(double** pvector,int** pindices);
    4748                PetscVec*   Duplicate(void);
    4849                void        Set(IssmDouble value);

  • issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h

    r22580 r23602  
    3030                double* col_partition_vector,int col_partition_vector_size);
    3131void PetscOptionsDetermineSolverType(int* psolver_type);
    32 void VecMerge(Vec A, Vec B, double* row_partition_vector,int row_partition_size);
    3332void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm);
    3433void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
Note: See TracChangeset for help on using the changeset viewer.