[petsc-users] Memory Distribution and Mesh Partition in Finite Element Method

[Tabrez Ali]

I do something similar, e.g., see lines 41-138 in 
https://bitbucket.org/stali/defmod/src/default/src/main.F90 which you 
can just copy/paste.

Using DMPlex might be a better option in the long run.

Tabrez

On 04/26/2016 06:50 PM, Jie Cheng wrote:
> Hi Barry
>
> Thanks for your answer. But given the mesh file I have, how do I partition the elements and vertices? MatCreateMPIAdj requires the graph information on adjacent elements (to set up ia and ja). Do I have to use more sophisticated meshing software to create mesh files that have graph information in it?
>
> Jie
>
>> On Apr 26, 2016, at 2:18 PM, Barry Smith<bsmith at mcs.anl.gov>  wrote:
>>
>>
>>    The "trick" is that first partition the element across processes, then partition the vertices (nodal values) subservient to the partitioning of the elements and then you "renumber" the elements and vertices so that the elements on the first process are numbered first, followed by the elements on the second process etc and similarly the vertices on the first process are numbered before the vertices on the second processes etc.
>>
>>   Now each process just needs to loop over its elements compute the element stiffness/load and call MatSetValues/VecSetValues() the the "new" numbering of the vertices.  The "old" numbering that was on the disk is simply not used in communicating with PETSc, you only use the new PETSc numbering.
>>
>>    Barry
>>
>>> On Apr 26, 2016, at 1:03 PM, Jie Cheng<chengj5 at rpi.edu>  wrote:
>>>
>>> Hello everyone
>>>
>>> I have a finite element code to solve nonlinear solid mechanics (using Newton-Raphson iteration) implemented with PETSc. The code is serial, organized as following:
>>>
>>> 1) Read the connectivity of the unstructured mesh, coordinates of nodes from individual txt files.
>>> 1.1) Connectivity file contains [# of elements] rows, and each row lists the global number of nodes on that element. Take 3d hexahedral elements for instance:
>>>       223 224 298 297 1 2 76 75
>>>       224 225 299 298 2 3 77 76
>>>       … …
>>> 1.2) Coordinates file contains [# of nodes] rows, and each row lists the coordinates of a node, for example:
>>>       0	            0.0011			3.9e-5
>>>       2.3677e-5     0.001.9975		3.9e-5
>>>       … …
>>>
>>> 2) Create the global stiffness A matrix with MatCreateSeqAIJ since the dimensions and nonzero pattern are known from the connectivity.
>>> 3) Loop over the element to compute the element stiffness matrix and right hand side. Then assemble the global stiffness matrix and right hand side.
>>> 4) Solve the linear equation with KSPsolve for the displacement increment, then go back to Step 3.
>>>
>>> The code works fine in serial, now I’m trying to parallelize it. To partition the mesh, I can use partdmesh from METIS, or let PETSc calls it. Either way I will find a way to assign different elements and nodes to different ranks. My question is: since PETSc does not allow us to control the memory distribution of the parallel matrix/vector, how do I make sure the rank happens to have all/most memory it needs for the specific elements? For example, rank 0 is in charged of element n, and needs to insert values to A[ i ][ j ], how do I make sure the i-th row is assigned to rank 0?
>>>
>>> This question is fundamental for people work with finite element methods. I checked the tutorial codes but did not find an example to follow. Section 3.5 of the manual talks about partition, but it does not say anything about linking the partition with the memory distribution. Could anyone give me some instructions on this issue? Thank you in advance!
>>>
>>> Best
>>> Jie Cheng