Changeset 23518

Timestamp:

12/07/18 15:30:01 (6 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added communications to vertices

Location:

issm/trunk-jpl/src/c

Files:

• classes/Inputs/DatasetInput.cpp (modified) (1 diff)
• classes/Vertices.cpp (modified) (2 diffs)
• classes/Vertices.h (modified) (1 diff)
• modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp

@@ -73 +73 @@
         MARSHALLING(enum_type);
         MARSHALLING(numids);
-        MARSHALLING_DYNAMIC(ids,int,numids)
+        MARSHALLING_DYNAMIC(ids,int,numids);
         if (marshall_direction == MARSHALLING_BACKWARD) inputs = new Inputs();
         inputs->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);

issm/trunk-jpl/src/c/classes/Vertices.cpp

@@ -25 +25 @@
 /*Object constructors and destructor*/
 Vertices::Vertices(){/*{{{*/
-        enum_type=VerticesEnum;
+        this->enum_type       = VerticesEnum;
+        this->common_recv     = NULL;
+        this->common_recv_ids = NULL;
+        this->common_send     = NULL;
+        this->common_send_ids = NULL;
         return;
 }
 /*}}}*/
 Vertices::~Vertices(){/*{{{*/
+
+        int num_proc=IssmComm::GetSize();
+
+        if(common_recv); xDelete<int>(common_recv);
+        if(common_send); xDelete<int>(common_send);
+        if(common_recv_ids){
+                for(int i=0;i<num_proc;i++) if(common_recv_ids[i]) xDelete<int>(common_recv_ids[i]);
+                xDelete<int*>(common_recv_ids);
+        }
+        if(common_send_ids){
+                for(int i=0;i<num_proc;i++) if(common_send_ids[i]) xDelete<int>(common_send_ids[i]);
+                xDelete<int*>(common_send_ids);
+        }
         return;
 }
@@ -35 +52 @@
 
 /*Numerics management*/
+void Vertices::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
+
+        int num_procs=IssmComm::GetSize();
+        int test = num_procs;
+        MARSHALLING_ENUM(VerticesEnum);
+        MARSHALLING(test);
+        if(test!=num_procs) _error_("number of cores is not the same as before");
+        MARSHALLING_DYNAMIC(this->common_recv,int,num_procs);
+        MARSHALLING_DYNAMIC(this->common_send,int,num_procs);
+        if(marshall_direction == MARSHALLING_BACKWARD){
+                this->common_recv_ids = xNew<int*>(num_procs);
+                this->common_send_ids = xNew<int*>(num_procs);
+        }
+        for(int i=0;i<num_procs;i++){
+                MARSHALLING_DYNAMIC(this->common_recv_ids[i],int,this->common_recv[i]);
+                MARSHALLING_DYNAMIC(this->common_send_ids[i],int,this->common_send[i]);
+        }
+        DataSet::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+}
+/*}}}*/
 void  Vertices::DistributePids(int numberofobjects){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Vertices.h

@@ -15 +15 @@
 
         public:
+                int*  common_recv;
+                int** common_recv_ids;
+                int*  common_send;
+                int** common_send_ids;
 
                 /*constructors, destructors:*/ 
                 Vertices();
                 ~Vertices();
+
+                /*Objects virtual functions*/
+                void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
 
                 /*numerics:*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

@@ -244 +244 @@
         iomodel->FindConstant(&isoceancoupling,"md.transient.isoceancoupling");
 
-        /*Fetch data:*/
+        /*Fetch data that will be used by the Vertex Constructor*/
         iomodel->FetchData(6,"md.mesh.x","md.mesh.y","md.mesh.z","md.geometry.base","md.geometry.thickness","md.mask.ice_levelset");
         if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
         else iomodel->FetchDataToInput(elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.);
         if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long");
-
         if (solution_type!=LoveSolutionEnum) CreateNumberNodeToElementConnectivity(iomodel);
 
-        int lid = 0;
+        const int MAXCONNECTIVITY = 5;
+        int*      epart = iomodel->epart;
+
+        /*Determine element width*/
+        int  elements_width;
+        switch(iomodel->meshelementtype){
+                case TriaEnum:  elements_width=3; break;
+                case TetraEnum: elements_width=4; break;
+                case PentaEnum: elements_width=6; break;
+                default: _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
+        }
+
+        /*Get my_rank:*/
+        int my_rank   = IssmComm::GetRank();
+        int num_procs = IssmComm::GetSize();
+
+        /*create matrix that keeps track of all ranks that have vertex i, and initialize as -1 (Common to all CPUs)*/
+        int* vertices_ranks = xNew<int>(MAXCONNECTIVITY*iomodel->numberofvertices);
+        for(int i=0;i<MAXCONNECTIVITY*iomodel->numberofvertices;i++) vertices_ranks[i] = -1;
+
+        /*For all vertices, cound how many cpus hold vertex i (initialize with 0)*/
+        int* vertices_proc_count = xNewZeroInit<int>(iomodel->numberofvertices);
+
+        /*Create vector of size total nbv, initialized with -1, that will keep track of local ids*/
+        int* vertices_lids  = xNew<int>(iomodel->numberofvertices);
+        for(int i=0;i<iomodel->numberofvertices;i++) vertices_lids[i] = -1;
+
+        /*Go through all elements and mark all vertices for all partitions*/
+        int  lid = 0;
+        for(int i=0;i<iomodel->numberofelements;i++){
+                for(int j=0;j<elements_width;j++){
+
+                        /*Get current vertex sid*/
+                        int vid = iomodel->elements[elements_width*i+j]-1;
+
+                        /*See if it has already been marked*/
+                        bool found = false;
+                        for(int k=0;k<vertices_proc_count[vid];k++){
+                                if(vertices_ranks[MAXCONNECTIVITY*vid+k] == epart[i]){
+                                        found = true;
+                                        break; 
+                                }
+                        }
+
+                        /*On go below if this vertex has not been seen yet in this partition*/
+                        if(!found){
+                                /*This rank has not been marked for this vertex just yet so go ahead and mark it*/
+                                vertices_ranks[MAXCONNECTIVITY*vid+vertices_proc_count[vid]] = epart[i];
+                                vertices_proc_count[vid]++;
+
+                                /*Keep track of local ids!*/
+                                if(epart[i]==my_rank){
+                                        vertices_lids[vid] = lid;
+                                        lid++;
+                                }
+
+                                /*Make sure we don't go too far in the table*/
+                                if(vertices_proc_count[vid]>MAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)");
+                        }
+                }
+        }
+
+        /*Now, Count how many clones we have with other partitions*/
+        int*  common_send = xNew<int>(num_procs);
+        int*  common_recv = xNew<int>(num_procs);
+        int** common_send_ids = xNew<int*>(num_procs);
+        int** common_recv_ids = xNew<int*>(num_procs);
+
+        /*First step: allocate, Step 2: populate*/
+        for(int step=0;step<2;step++){
+
+                if(step==1){
+                        /*Allocate send and receive arrays of ids now*/
+                        for(int i=0;i<num_procs;i++){
+                                _assert_(common_send[i]>=0 && common_recv[i]>=0);
+                                common_send_ids[i] = xNew<int>(common_send[i]);
+                                common_recv_ids[i] = xNew<int>(common_recv[i]);
+                        }
+                }
+
+                /*Re/Initialize counters to 0*/
+                for(int i=0;i<num_procs;i++){
+                        common_recv[i]=0;
+                        common_send[i]=0;
+                }
+
+                /*Go through table and find clones/masters etc*/
+                for(int i=0;i<iomodel->numberofvertices;i++){
+
+                        /*If we did not find this vertex in our current partition, go to next vertex*/
+                        if(vertices_lids[i] == -1) continue;
+
+                        /*Find in what column this rank belongs*/
+                        int col = -1;
+                        for(int j=0;j<MAXCONNECTIVITY;j++){
+                                if(vertices_ranks[MAXCONNECTIVITY*i+j] == my_rank){
+                                        col = j;
+                                        break;
+                                }
+                        }
+                        _assert_(col!=-1);
+
+                        /*If col==0, it is either not on boundary, or a master*/
+                        if(col==0){
+                                /*1. is this vertex on the boundary? Skip if not*/
+                                if(vertices_ranks[MAXCONNECTIVITY*i+col+1]==-1){
+                                        continue;
+                                }
+                                else{
+                                        for(int j=1;j<vertices_proc_count[i];j++){
+                                                _assert_(vertices_ranks[MAXCONNECTIVITY*i+j]>=0);
+                                                int rank = vertices_ranks[MAXCONNECTIVITY*i+j];
+                                                if(step==1){
+                                                        common_send_ids[rank][common_send[rank]] = vertices_lids[i];
+                                                }
+                                                common_send[rank]++;
+                                        }
+                                }
+                        }
+                        else{
+                                /*3. It is a slave, record that we need to receive for this cpu*/
+                                int rank = vertices_ranks[MAXCONNECTIVITY*i+0];
+                                if(step==1){
+                                        common_recv_ids[rank][common_recv[rank]] = vertices_lids[i];
+                                }
+                                common_recv[rank]++;
+                        }
+                }
+        }
+        xDelete<int>(vertices_proc_count);
+        xDelete<int>(vertices_ranks);
+
         for(int i=0;i<iomodel->numberofvertices;i++){
-                if(iomodel->my_vertices[i]) vertices->AddObject(new Vertex(i+1,i,lid++,i,iomodel));
-        }
+                if(vertices_lids[i]!=-1) vertices->AddObject(new Vertex(i+1,i,vertices_lids[i],i,iomodel));
+        }
+        xDelete<int>(vertices_lids);
 
         /*Free data: */
@@ -261 +392 @@
         if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->DeleteData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
         if (isoceancoupling) iomodel->DeleteData(2,"md.mesh.lat","md.mesh.long");
+
+        vertices->common_send=common_send;
+        vertices->common_recv=common_recv;
+        vertices->common_send_ids=common_send_ids;
+        vertices->common_recv_ids=common_recv_ids;
 }/*}}}*/