Changeset 3378


Ignore:
Timestamp:
04/02/10 15:06:53 (15 years ago)
Author:
Mathieu Morlighem
Message:

fixed DG in parallel

Location:
issm/trunk/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp

    r3373 r3378  
    1818        extern int my_rank;
    1919       
    20         _printf_("      Configuring elements...\n");
     20        //_printf_("      Configuring elements...\n");
    2121        elements->Configure(elements,loads,nodes,materials,parameters);
    22         _printf_("      Configuring loads...\n");
    23         //if (my_rank==2){
    24         //      nodes->Echo();
    25         //}
     22        //_printf_("      Configuring loads...\n");
    2623        loads->Configure(elements,loads,nodes,materials,parameters);
    27         _printf_("      Configuring nodes...\n");
     24        //_printf_("      Configuring nodes...\n");
    2825        nodes->Configure(elements,loads,nodes,materials,parameters);
    29         _printf_("      Configuring parameters...\n");
     26        //_printf_("      Configuring parameters...\n");
    3027        parameters->Configure(elements,loads, nodes, materials,parameters);
    3128

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp

    r3373 r3378  
    22 * CreateElementsNodesAndMaterialsPrognostic2.c:
    33 */
    4 
    54
    65#include "../../DataSet/DataSet.h"
    … …  
    109#include "../../shared/shared.h"
    1110#include "../../include/macros.h"
     11#include "../../include/typedefs.h"
    1212#include "../../MeshPartitionx/MeshPartitionx.h"
    1313#include "../IoModel.h"
    14 
    1514
    1615void    CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    … …  
    2120
    2221        /*DataSets: */
    23         DataSet*    elements  = NULL;
    24         DataSet*    nodes = NULL;
    25         DataSet*    materials = NULL;
     22        DataSet* elements  = NULL;
     23        DataSet* nodes = NULL;
     24        DataSet* materials = NULL;
    2625       
    2726        /*Objects: */
    28         Node*       node   = NULL;
    29         Tria*       tria = NULL;
    30         Penta*      penta = NULL;
    31         Matice*     matice  = NULL;
    32         Matpar*     matpar  = NULL;
     27        Node*   node   = NULL;
     28        Tria*   tria = NULL;
     29        Penta*  penta = NULL;
     30        Matice* matice  = NULL;
     31        Matpar* matpar  = NULL;
    3332
    3433        /*output: */
    … …  
    124123        int  gridcount;
    125124        int  count;
     125        /*Nodes cloning*/
     126        double e1,e2;
     127        int i1,i2;
     128        int pos;
    126129        #endif
    127130        int  first_grid_index;
    … …  
    183186                #endif
    184187                       
    185                        
    186188                        /*ids: */
    187189                        tria_id=i+1; //matlab indexing.
    … …  
    233235
    234236                }//for (i=0;i<numberofelements;i++)
    235 
    236237       
    237238                /*Free data : */
    238                 xfree((void**)&iomodel->elements);
    239239                xfree((void**)&iomodel->thickness);
    240240                xfree((void**)&iomodel->surface);
    … …  
    249249
    250250        #ifdef _PARALLEL_
    251                 /*From the element partitioning, we can determine which grids are on the inside of this cpu's
    252                  *element partition, and which are on its border with other nodes:*/
    253                 gridborder=NewVec(3*iomodel->numberofelements);
    254 
    255                 for (i=0;i<3*iomodel->numberofelements;i++){
    256                         if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
     251        /*If we are in parallel, we must add the nodes that are not in this partition
     252         * but share edges of elements that are in this partition. This is needed for
     253         * the loads as numerical fluxes involve the dofs of all nodes on a given edge*/
     254
     255        /*Get edges and elements*/
     256        IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
     257
     258        /*!All elements have been partitioned above, only create elements for this CPU: */
     259        for (i=0;i<iomodel->numberofedges;i++){
     260
     261                /*Get left and right elements*/
     262                e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
     263                e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
     264
     265                /* 1) If the element e1 is in the current partition
     266                 * 2) and if the edge of the element is shared by another element
     267                 * 3) and if this element is not in the same partition:
     268                 * we must clone the nodes on this partition so that the loads (Numericalflux)
     269                 * will have access to their properties (dofs,...)*/
     270                if(my_rank==epart[(int)e1] && !isnan(e2) && my_rank!=epart[(int)e2]){
     271
     272                        /*1: Get vertices ids*/
     273                        i1=(int)iomodel->edges[4*i+0];
     274                        i2=(int)iomodel->edges[4*i+1];
     275
     276                        /*2: Get the column where these ids are located in the index*/
     277                        pos==UNDEF;
     278                        for(j=0;j<3;j++){
     279                                if (iomodel->elements[3*(int)e2+j]==i1) pos=j+1;
     280                        }
     281                        ISSMASSERT(pos!=UNDEF);
     282
     283                        /*3: We have the id of the elements and the position of the vertices in the index
     284                         * we can now create the corresponding nodes:*/
     285                        my_grids[(int)e2*3+pos-1]=1;
     286                        my_grids[(int)e2*3+((pos+1)%3)]=1;
    257287                }
    258                 VecAssemblyBegin(gridborder);
    259                 VecAssemblyEnd(gridborder);
    260 
    261                 #ifdef _ISSM_DEBUG_
    262                 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
    263                 #endif
    264                
    265                 VecToMPISerial(&my_bordergrids,gridborder);
    266 
    267                 #ifdef _ISSM_DEBUG_
    268                 if(my_rank==0){
    269                         for (i=0;i<iomodel->numberofnodes;i++){
    270                                 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
    271                         }
    272                 }
    273                 #endif
     288        }
     289
     290        /*Free data: */
     291        xfree((void**)&iomodel->edges);
     292
     293        /*From the element partitioning, we can determine which nodes are on the inside of this cpu's
     294         *element partition, and which are on its border with other nodes:*/
     295        gridborder=NewVec(3*iomodel->numberofelements);
     296
     297        for (i=0;i<3*iomodel->numberofelements;i++){
     298                if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
     299        }
     300        VecAssemblyBegin(gridborder);
     301        VecAssemblyEnd(gridborder);
     302
     303        VecToMPISerial(&my_bordergrids,gridborder);
     304
    274305        #endif
    275306
    … …  
    295326        materials->AddObject(matpar);
    296327
    297         /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
     328        /*Ok, let's summarise. Now, every CPU has the following array: my_grids
    298329         We can therefore determine  which grids are internal to this node's partition
    299330         and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
    … …  
    307338                IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
    308339        }
    309         IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    310340        IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
    311341        IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
    … …  
    318348        IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
    319349
    320 
    321350        /*Get number of dofs per node: */
    322351        DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
    … …  
    328357                        #ifdef _PARALLEL_
    329358                        /*!All elements have been partitioned above, only create elements for this CPU: */
    330                         if(my_rank==epart[i]){
     359                        if(my_grids[3*i+j]){
    331360                        #endif
    332361
    333                         //Get id of the vertex on which the current node is located
    334                         k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing)
    335                         ISSMASSERT(k>=0 && k<iomodel->numberofnodes);
    336 
    337                         //Get node properties
    338                         node_id=i*3+j+1;
    339                         node_partitionborder=0;
    340                         node_x[0]=iomodel->x[k];
    341                         node_x[1]=iomodel->y[k];
    342                         node_x[2]=iomodel->z[k];
    343                         node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]);
    344                         node_onbed=(int)iomodel->gridonbed[k];
    345                         node_onsurface=(int)iomodel->gridonsurface[k]; 
    346                         node_onshelf=(int)iomodel->gridoniceshelf[k];   
    347                         node_onsheet=(int)iomodel->gridonicesheet[k];   
    348                         if (strcmp(iomodel->meshtype,"3d")==0){
    349                                 if (isnan(iomodel->uppernodes[k])){
    350                                         node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
     362                                //Get id of the vertex on which the current node is located
     363                                k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing)
     364                                ISSMASSERT(k>=0 && k<iomodel->numberofnodes);
     365
     366                                //Get node properties
     367                                node_id=i*3+j+1;
     368                                #ifdef _PARALLEL_
     369                                if(my_bordergrids[node_id-1]>1.0) { //this grid belongs to a partition border
     370                                        node_partitionborder=1;
    351371                                }
    352372                                else{
    353                                         node_upper_node_id=(int)iomodel->uppernodes[k];
     373                                        node_partitionborder=0;
    354374                                }
    355                         }
    356                         else{
    357                                 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
    358                                 node_upper_node_id=node_id;
    359                         }
    360 
    361                         /*Create node using its constructor: */
    362                         node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
    363 
    364                         /*Add node to nodes dataset: */
    365                         nodes->AddObject(node);
    366 
     375                                #else
     376                                        node_partitionborder=0;
     377                                #endif
     378                                node_x[0]=iomodel->x[k];
     379                                node_x[1]=iomodel->y[k];
     380                                node_x[2]=iomodel->z[k];
     381                                node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]);
     382                                node_onbed=(int)iomodel->gridonbed[k];
     383                                node_onsurface=(int)iomodel->gridonsurface[k]; 
     384                                node_onshelf=(int)iomodel->gridoniceshelf[k];   
     385                                node_onsheet=(int)iomodel->gridonicesheet[k];   
     386                                if (strcmp(iomodel->meshtype,"3d")==0){
     387                                        if (isnan(iomodel->uppernodes[k])){
     388                                                node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
     389                                        }
     390                                        else{
     391                                                node_upper_node_id=(int)iomodel->uppernodes[k];
     392                                        }
     393                                }
     394                                else{
     395                                        /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
     396                                        node_upper_node_id=node_id;
     397                                }
     398
     399                                /*Create node using its constructor: */
     400                                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
     401
     402                                /*Add node to nodes dataset: */
     403                                nodes->AddObject(node);
    367404                        #ifdef _PARALLEL_
    368405                        }

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp

    r3371 r3378  
    2828        int  numericalflux_id;
    2929        int  numericalflux_node_ids[MAX_NUMERICALFLUX_NODES];
    30         int  numericalflux_elem_ids[MAX_NUMERICALFLUX_ELEMS];
     30        int  numericalflux_elem_id;
    3131
    3232        /*Create loads: */
    … …  
    5858                        strcpy(numericalflux_type,"internal");
    5959
    60                         numericalflux_elem_ids[0]=(int)e1+1;//id is in matlab index
    61                         numericalflux_elem_ids[1]=(int)e2+1;//id is in matlab index
     60                        numericalflux_elem_id=(int)e1+1;//id is in matlab index
    6261
    6362                        /*Now, we must get the nodes of the 4 nodes located on the edge*/
    … …  
    8584                        strcpy(numericalflux_type,"boundary");
    8685
    87                         numericalflux_elem_ids[0]=(int)e1+1;
     86                        numericalflux_elem_id=(int)e1+1;
    8887
    8988                        /*1: Get vertices ids*/
    … …  
    104103                }
    105104
    106                 numericalflux = new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_ids);
     105                numericalflux = new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id);
    107106
    108107                loads->AddObject(numericalflux);

  • issm/trunk/src/c/objects/Node.cpp

    r3372 r3378  
    220220        idxm=(id-1);
    221221        value=(double)doflist1[0];
     222        ISSMASSERT(value>=0);
    222223
    223224        VecSetValues(partition,1,&idxm,&value,INSERT_VALUES);

  • issm/trunk/src/c/objects/Numericalflux.cpp

    r3371 r3378  
    1717#include "./objects.h"
    1818
     19extern int my_rank;
     20
    1921/*Object constructors and destructor*/
    2022/*FUNCTION Numericalflux::Numericalflux(){{{1*/
    … …  
    2426/*}}}*/
    2527/*FUNCTION Numericalflux::Numericalflux(char numericalflux_type[NUMERICALFLUXSTRING],int numericalflux_fill...){{{1*/
    26 Numericalflux::Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_ids[MAX_NUMERICALFLUX_ELEMS]){
     28Numericalflux::Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_id){
    2729
    2830        int i;
    … …  
    3032        strcpy(type,numericalflux_type);
    3133        id=numericalflux_id;
    32        
    33         for(i=0;i<MAX_NUMERICALFLUX_ELEMS;i++){
    34                 element_ids[i]=numericalflux_element_ids[i];
    35                 element_offsets[i]=UNDEF;
    36                 elements[i]=NULL;
    37         }
     34
     35        element_id=numericalflux_element_id;
     36        element_offset=UNDEF;
     37        element=NULL;
    3838
    3939        for(i=0;i<MAX_NUMERICALFLUX_NODES;i++){
    … …  
    7373        memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
    7474
    75         memcpy(&element_ids,marshalled_dataset,sizeof(element_ids));marshalled_dataset+=sizeof(element_ids);
    76         memcpy(&element_offsets,marshalled_dataset,sizeof(element_offsets));marshalled_dataset+=sizeof(element_offsets);
    77         for(i=0;i<MAX_NUMERICALFLUX_ELEMS;i++)elements[i]=NULL;
     75        memcpy(&element_id,marshalled_dataset,sizeof(element_id));marshalled_dataset+=sizeof(element_id);
     76        memcpy(&element_offset,marshalled_dataset,sizeof(element_offset));marshalled_dataset+=sizeof(element_offset);
     77        element=NULL;
    7878       
    7979        memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
    … …  
    105105        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
    106106       
    107         memcpy(marshalled_dataset,&element_ids,sizeof(element_ids));marshalled_dataset+=sizeof(element_ids);
    108         memcpy(marshalled_dataset,&element_offsets,sizeof(element_offsets));marshalled_dataset+=sizeof(element_offsets);
     107        memcpy(marshalled_dataset,&element_id,sizeof(element_id));marshalled_dataset+=sizeof(element_id);
     108        memcpy(marshalled_dataset,&element_offset,sizeof(element_offset));marshalled_dataset+=sizeof(element_offset);
    109109
    110110        memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
    … …  
    120120        return sizeof(type)+
    121121                sizeof(id)+
    122                 sizeof(element_ids)+
    123                 sizeof(element_offsets)+
     122                sizeof(element_id)+
     123                sizeof(element_offset)+
    124124                sizeof(node_ids)+
    125125                sizeof(node_offsets)+
    … …  
    142142        if (strcmp(type,"internal")==0){
    143143                ResolvePointers((Object**)nodes,node_ids,node_offsets,4,nodesin);
    144                 ResolvePointers((Object**)elements,element_ids,element_offsets,2,elementsin);
    145144        }
    146145        else if (strcmp(type,"boundary")==0){
    147146                ResolvePointers((Object**)nodes,node_ids,node_offsets,2,nodesin);
    148                 ResolvePointers((Object**)elements,element_ids,element_offsets,1,elementsin);
    149147        }
    150148        else ISSMERROR("type not supported yet");
     149        ResolvePointers((Object**)&element,&element_id,&element_offset,1,elementsin);
    151150
    152151}
    … …  
    176175        printf("   id: %i\n",id);
    177176       
     177        printf("   element_id=%i\n",element_id);
     178        printf("   element_offset=%i\n",element_offset);
     179        if(element)element->Echo();
    178180        if (strcmp(type,"internal")==0){
    179                 printf("   element_ids=[%i,%i]\n",element_ids[0],element_ids[1]);
    180                 printf("   element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]);
    181                 for(i=0;i<2;i++){
    182                         if(elements[i])elements[i]->Echo();
    183                 }
    184181                printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
    185182                printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
    … …  
    189186        }
    190187        else{
    191                 printf("   element_ids=[%i,%i]\n",element_ids[0],element_ids[1]);
    192                 printf("   element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]);
    193                 for(i=0;i<1;i++){
    194                         if(elements[i])elements[i]->Echo();
    195                 }
    196188                printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
    197189                printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
    … …  
    215207        printf("   id: %i\n",id);
    216208
     209        printf("   element_id=%i\n",element_id);
     210        printf("   element_offset=%i]\n",element_offset);
     211
    217212        if (strcmp(type,"internal")==0){
    218                 printf("   element_ids=[%i,%i]\n",element_ids[0],element_ids[1]);
    219                 printf("   element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]);
    220213                printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
    221214                printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
    222215        }
    223216        else{
    224                 printf("   element_ids=%i\n",element_ids[0]);
    225                 printf("   element_offsets=%i\n",element_offsets[0]);
    226217                printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
    227218                printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);

  • issm/trunk/src/c/objects/Numericalflux.h

    r3371 r3378  
    2424       
    2525                /*elements: */
    26                 Element* elements[MAX_NUMERICALFLUX_ELEMS];
    27                 int      element_ids[MAX_NUMERICALFLUX_ELEMS];
    28                 int      element_offsets[MAX_NUMERICALFLUX_ELEMS];
     26                Element* element;
     27                int      element_id;
     28                int      element_offset;
    2929
    3030                /*nodes: */
    … …  
    3636
    3737                Numericalflux();
    38                 Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_ids[MAX_NUMERICALFLUX_ELEMS]);
     38                Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_id);
    3939                Object* copy();
    4040                ~Numericalflux();

  • issm/trunk/src/c/parallel/prognostic2.cpp

    r3373 r3378  
    6161        MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
    6262
    63         _printf_("recover , input file name and output file name:\n");
     63        _printf_("recover input file name and output file name:\n");
    6464        inputfilename=argv[2];
    6565        outputfilename=argv[3];

  • issm/trunk/src/c/parallel/prognostic2_core.cpp

    r3373 r3378  
    2828        int numberofdofspernode;
    2929        int numberofnodes;
    30         int dofs[1]={0};
     30        int dofs[1]={1};
    3131
    3232        /*fem prognostic model: */
    3333        FemModel* fem_p=NULL;
    34 
    3534
    3635        /*recover fem model: */

  • issm/trunk/src/m/classes/public/mesh/meshconvert.m

    r3354 r3378  
    3131md.y=bamgmesh_out.Vertices(:,2);
    3232md.elements=bamgmesh_out.Triangles(:,1:3);
     33md.edges=bamgmesh_out.ElementEdges;
    3334md.segments=bamgmesh_out.Segments(:,1:3);
    3435md.segmentmarkers=bamgmesh_out.Segments(:,4);

  • issm/trunk/src/m/solutions/jpl/prognostic2_core.m

    r3372 r3378  
    2323
    2424        displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']);
    25         results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');
     25        %results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');
    2626
    2727end %end function
Note: See TracChangeset for help on using the changeset viewer.