Matrix assembly in fem using mpi
following my previous question regarding using MPI for finite element
http://www.cplusplus.com/forum/general/278132/ I have made a few changes in my code. First of all instead of dividing mesh into pieces, I copied entire mesh into all processors in order to avoid using hole data. I know this is not an optimized way for doing MPI but at the time being it works for me, because in the next step I am going to use CSR format for matrix assembly. After assembly I get strange number in global stiffness matrix which is -1.25549e+67 and I don't know why. I would be grateful if somebody could help me.
http://www.cplusplus.com/forum/general/278132/ I have made a few changes in my code. First of all instead of dividing mesh into pieces, I copied entire mesh into all processors in order to avoid using hole data. I know this is not an optimized way for doing MPI but at the time being it works for me, because in the next step I am going to use CSR format for matrix assembly. After assembly I get strange number in global stiffness matrix which is -1.25549e+67 and I don't know why. I would be grateful if somebody could help me.
|
|
Last edited on
@resabzr
The whole point of using distributed-memory computing is that you don't send everything to all processors.
Also, you have already been told (several times) on StackOverflow to use contiguous memory - not your (current) **double pointer arrangement. The latter will hopelessly slow down any MPI transfers. Either flatten your dynamic arrays to 1-d, or use the data buffer of std::vector.
If you absolutely must access elements in a 2-d array manner then, for contiguous memory you will have to allocate memory as @dutch does in this post:
http://www.cplusplus.com/forum/beginner/267829/#msg1152277
Note that this is essentially a flattened 1-d array "under the hood".
You shouldn't use MPI_Bcast node-by-node; that would be horribly slow. You also can't use whole-array MPI transfers unless your data is in contiguous memory, which it won't be with your (current) double-pointer arrangement. Similarly, I wouldn't guarantee that your MPI_Reduce will work (though I can't see what it would be used for in this code, anyway).
FWIW you have a data indexing error in two of the terms on line 249:
but there's no way of knowing whether this is the cause of your error, because there is no way of testing your code on an input file.
If you want to learn how to program the finite-element method then you can read
"Programming the Finite Element Method" by Smith, Griffiths and Margetts:
https://onlinelibrary.wiley.com/doi/epub/10.1002/9781119189237
and yeah, two of the authors are my (retired or current) work colleagues.
The whole point of using distributed-memory computing is that you don't send everything to all processors.
Also, you have already been told (several times) on StackOverflow to use contiguous memory - not your (current) **double pointer arrangement. The latter will hopelessly slow down any MPI transfers. Either flatten your dynamic arrays to 1-d, or use the data buffer of std::vector.
If you absolutely must access elements in a 2-d array manner then, for contiguous memory you will have to allocate memory as @dutch does in this post:
http://www.cplusplus.com/forum/beginner/267829/#msg1152277
Note that this is essentially a flattened 1-d array "under the hood".
You shouldn't use MPI_Bcast node-by-node; that would be horribly slow. You also can't use whole-array MPI transfers unless your data is in contiguous memory, which it won't be with your (current) double-pointer arrangement. Similarly, I wouldn't guarantee that your MPI_Reduce will work (though I can't see what it would be used for in this code, anyway).
FWIW you have a data indexing error in two of the terms on line 249:
local_NN[node_1][node_1] = NN[node_1][node_1] + NN_e[1][2];but there's no way of knowing whether this is the cause of your error, because there is no way of testing your code on an input file.
If you want to learn how to program the finite-element method then you can read
"Programming the Finite Element Method" by Smith, Griffiths and Margetts:
https://onlinelibrary.wiley.com/doi/epub/10.1002/9781119189237
and yeah, two of the authors are my (retired or current) work colleagues.
Last edited on
I do appreciate your answer. But I think you are mistaken. I have NOT been told in stuckoverflow several times about 1d contiguous array. anyway, I need the 2d array in advance. If you look at my previous code I defined 1d contiguous array, but I needed a 2d array, that's why I am using double pointer. About mpi reduce, imagine we have 2 matrix in 2 process,
A={{1,2},{0,0}};
B={{0,0},{1,2}};
I want mpi reduce do this for me,
A+B={{1,2},{1,2}}
A={{1,2},{0,0}};
B={{0,0},{1,2}};
I want mpi reduce do this for me,
A+B={{1,2},{1,2}}
Last edited on
@resabzr
No, you do not understand.
Instead of a 2-d array, say
A = [ [ 1, 2 ], [ 3, 4 ] ]
you can simply "flatten" this to
A = [ 1, 2, 3, 4 ]
and if you want the equivalent 2-d indices then, as you know the number of columns, then you can easily work them out. But the key advantage is that the latter form can be MPI-transferred in a SINGLE communication - you don't have to do it row by row. MPI communications are slow. Imagine what would happen if you were working in 3-d, not 2-d, with your multiple-pointer notation.
I can see that your MPI_Reduce will (in this code) do a sum. However, at the moment, all the processors seem to know everything, so all you will be doing is summing the same thing up lots of times. You haven't decomposed your domain.
I assume that the index problem that I pointed out didn't solve the problem? Does the problem still occur if you just use a single processor?
I had a look at the various answers you've had in Stack Overflow (to the same codes that you have posted here). Those answers include these gems:
BTW, the term is "halo data", not "hole data" (this thread).
Likewise, "contiguous", not "contagious" (one of your SO posts).
No, you do not understand.
Instead of a 2-d array, say
A = [ [ 1, 2 ], [ 3, 4 ] ]
you can simply "flatten" this to
A = [ 1, 2, 3, 4 ]
and if you want the equivalent 2-d indices then, as you know the number of columns, then you can easily work them out. But the key advantage is that the latter form can be MPI-transferred in a SINGLE communication - you don't have to do it row by row. MPI communications are slow. Imagine what would happen if you were working in 3-d, not 2-d, with your multiple-pointer notation.
I can see that your MPI_Reduce will (in this code) do a sum. However, at the moment, all the processors seem to know everything, so all you will be doing is summing the same thing up lots of times. You haven't decomposed your domain.
I assume that the index problem that I pointed out didn't solve the problem? Does the problem still occur if you just use a single processor?
| resabzr wrote: |
|---|
| But I think you are mistaken. I have NOT been told in stuckoverflow several times about 1d contiguous array |
I had a look at the various answers you've had in Stack Overflow (to the same codes that you have posted here). Those answers include these gems:
all these raw poiters in modern C++ make my head spin... The matrices are not contiguous At the very least, you need to allocate your 2D array(s) in contiguous memory. Assuming your matrix is in contiguous memory, all you need is MPI_Reduce(..., MPI_SUM, ...) |
BTW, the term is "halo data", not "hole data" (this thread).
Likewise, "contiguous", not "contagious" (one of your SO posts).
Last edited on
I got your point about single block of memory, I understand this
And I used this method http://www.cplusplus.com/forum/beginner/267829/#msg1152277.
I fixed that index problem but it did not make any change in the results. I wrote the same code in c++ and I'm getting the right answer, but here even with 1 process everything is different. Even when I print out the area of each element it gives me zero, or when I print out the gradient of shape functions it gives me non(ind). How it can be possible even on one process results change like this in comparison with c++ code without mpi implementation.
And I used this method http://www.cplusplus.com/forum/beginner/267829/#msg1152277.
I fixed that index problem but it did not make any change in the results. I wrote the same code in c++ and I'm getting the right answer, but here even with 1 process everything is different. Even when I print out the area of each element it gives me zero, or when I print out the gradient of shape functions it gives me non(ind). How it can be possible even on one process results change like this in comparison with c++ code without mpi implementation.
Supply:
(1) your latest code (below this post)
(2) the input file (preferably small)
Make sure that any hard-coded cell counts correspond to the numbers in the input file.
(1) your latest code (below this post)
(2) the input file (preferably small)
Make sure that any hard-coded cell counts correspond to the numbers in the input file.
this is the final code
|
|
Last edited on
and this the mesh and connectivity
MESH dimension 3 ElemType Triangle Nnode 3
Coordinates
MESH dimension 3 ElemType Triangle Nnode 3
Coordinates
|
|
Last edited on
What exactly do your files connec.txt and mesh.txt contain. These are the files being read by your program.
Note, you should remove the lines
and the corresponding end brace, or the code will attempt (and fail) to read your coordinates a second time.
Note, you should remove the lines
|
|
and the corresponding end brace, or the code will attempt (and fail) to read your coordinates a second time.
Last edited on
the data I put here as coordinate and element actually are two separate file named connec.txt and mesh.txt. I just copy them here so you can see the data, I didn't know how put files here. I just did what you said but no difference.
Last edited on
Just cut and paste it from notepad or whatever.
Add outpad tags around it (second item on the format menu).
Add outpad tags around it (second item on the format menu).
OK, change the lines
to
Make sure that the number of elements corresponds to that in the input file.
|
|
to
|
|
Make sure that the number of elements corresponds to that in the input file.
Last edited on
not working, I did these changes but not working. I think the problem start when area of each element yield zero, when I print the element area it's all zero
Last edited on
The areas will come out 0 because you don't set arrays x[] or y[]. Did you mean something to do with coordinate[][]?
With two processors your code runs and produces (from your 14-element files):
What other output do you expect?
Note that some of your printf statements have the wrong number of parameters.
With two processors your code runs and produces (from your 14-element files):
C:\c++>"C:\Program Files\Microsoft MPI\bin"\mpiexec -n 2 temp 0 11 9 12 1 0 2 1 2 8 10 7 3 7 10 12 4 8 7 3 5 3 7 4 6 8 3 6 7 6 3 2 8 4 1 3 9 5 6 2 10 13 12 10 11 9 4 7 12 3 1 2 13 9 7 12 0 2 0 0 1 1.5 0 0 2 2 0.5 0 3 1.28472 0.501861 0 4 1 0 0 5 2 1 0 6 1.5 1 0 7 0.708333 0.511165 0 8 1 1 0 9 0.5 0 0 10 0.5 1 0 11 0 0 0 12 0 0.5 0 13 0 1 0 process 0 prainting connectivity: 0 11 9 12 1 0 2 1 2 8 10 7 3 7 10 12 4 8 7 3 5 3 7 4 6 8 3 6 7 6 3 2 8 4 1 3 9 5 6 2 10 13 12 10 11 9 4 7 12 3 1 2 13 9 7 12 process 0 prainting coordinate: 0 2 0 0 1 1.5 0 0 2 2 0.5 0 3 1.28472 0.501861 0 4 1 0 0 5 2 1 0 6 1.5 1 0 7 0.708333 0.511165 0 8 1 1 0 9 0.5 0 0 10 0.5 1 0 11 0 0 0 12 0 0.5 0 13 0 1 0 process 1 prainting connectivity: 0 11 9 12 1 0 2 1 2 8 10 7 3 7 10 12 4 8 7 3 5 3 7 4 6 8 3 6 7 6 3 2 8 4 1 3 9 5 6 2 10 13 12 10 11 9 4 7 12 3 1 2 13 9 7 12 process 1 prainting coordinate: 0 2 0 0 1 1.5 0 0 2 2 0.5 0 3 1.28472 0.501861 0 4 1 0 0 5 2 1 0 6 1.5 1 0 7 0.708333 0.511165 0 8 1 1 0 9 0.5 0 0 10 0.5 1 0 11 0 0 0 12 0 0.5 0 13 0 1 0 |
What other output do you expect?
Note that some of your printf statements have the wrong number of parameters.
Last edited on
The different number is because I changed the mesh to post here. It was not 14 points, it is 246 as you can see from the link I've shared. The problem is not the Broadcasting part. It works properly. The problem is in matrix assembly I get strange number. like this -6.27744e+66. All the elements area are zero, all the shape function gradient are non(ind). and I don't know why. The same code without MPI is working properly
Last edited on
Your number of elements (hard-coded) must match that in the input file.
Your code now runs.
Areas will work out as 0 (or, possibly, garbage) because you have not put anything in arrays x[] and y[].
Your code now runs.
Areas will work out as 0 (or, possibly, garbage) because you have not put anything in arrays x[] and y[].
Topic archived. No new replies allowed.