[petsc-users] Problem with MatZeroRowsColumnsIS()
Püsök, Adina-Erika
puesoek at uni-mainz.de
Fri May 9 09:50:36 CDT 2014
Yes, I tested the implementation with both MatZeroRowsIS() and MatZeroRowsColumnsIS(). But first, I will be more explicit about the problem I was set to solve:
We have a Dirichlet block of size (L,W,H) and centered (xc,yc,zc), which is much smaller than the model domain, and we set Vx = Vpush, Vy=0 within the block (Vz is let free for easier convergence).
As I said before, since the code does not have a monolithic matrix, but 4 submatrices (VV VP; PV PP), and the rhs has 2 sub vectors rhs=(f; g), my approach is to modify only (VV, VP, f) for the Dirichlet BC.
The way I tested the implementation:
1) Output (VV, VP, f, Dirichlet dofs) - unmodified (no Dirichlet BC)
2) Output (VV, VP, f, Dirichlet dofs) - a) modified with MatZeroRowsIS(),
- b) modified with MatZeroRowsColumnsIS() -> S_PETSc
Again, the only difference between a) and b) is:
// ierr = MatZeroRowsColumnsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
// ierr = MatZeroRowsColumnsIS(VV_MAT,isy,v_vv,x_push,rhs); CHKERRQ(ierr);
ierr = MatZeroRowsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
ierr = MatZeroRowsIS(VV_MAT,isy,v_vv,x_push,rhs); CHKERRQ(ierr);
3) Read them in Matlab and perform the exact same operations on the unmodified matrices and f vector. -> S_Matlab
4) Compare S_PETSc with S_Matlab. If the implementation is correct, they should be equal (VV, VP, f).
5) Check for 1 cpu and 4 cpus.
Now to answer your questions:
a,b,d) Yes, matrix modification is done correctly (check the spy diagrams below) in all cases: MatZeroRowsIS() and MatZeroRowsColumnsIS() on 1 and 4 cpus.
I should have said that in the piece of code above:
v_vv = 1.0;
v_vp = 0.0;
The vector x_push is a duplicate of rhs, with zero elements except the values for the Dirichlet dofs.
c) The rhs is a different matter. With MatZeroRows() there is no problem. The rhs is equivalent with the one in Matlab, sequential and parallel.
However, with MatZeroRowsColumns(), the residual contains nonzero elements, and in parallel the nonzero pattern is even bigger (1 cpu - 63, 4 cpu - 554). But if you look carefully, the values of the nonzero residuals are very small < +/- 1e-10.
So, I did a tolerance filter:
tol = 1e-10;
res = f_petsc - f_mod_matlab;
for i=1:length(res)
if abs(res(i))>0 & abs(res(i))<tol
res(i)=0;
end
end
and then the f_petsc and f_mod_matlab are equivalent on 1 and 4 cpus (figure 5). So it seems that MatZeroRowsColumnsIS() might give some nonzero residuals.
Calculating the norm2 of the residuals defined above in each case gives:
MatZeroRowsIS() 1cpu: norm(res,2) = 0
MatZeroRowsIS() 4cpu: norm(res,2) = 0
MatZeroRowsColumnsIS() 1cpu: norm(res,2) = 1.6880e-10
MatZeroRowsColumnsIS() 4cpu: norm(res,2) = 7.3786e+06
Since this is purely a problem of matrix and vector assembly/manipulation, I think the nonzero residuals of the rhs with MatZeroRowsColumnsIS() give the parallel artefacts that I showed last time.
If you need the raw data and the matlab scripts that I used for testing for your consideration, please let me know.
Thanks,
Adina
When performing the manual operations on the unmodified matrices and rhs vector in Matlab, I took into account:
- matlab indexing = petsc indexing +1;
- the vectors written to file for matlab (PETSC_VIEWER_BINARY_MATLAB) have the natural ordering, rather than the petsc ordering. On 1 cpu, they are equivalent, but on 4 cpus, the Dirichlet BC indices had to be converted to natural indices in order to perform the correct operations on the rhs.
[cid:8547B616-9585-4E29-BC0D-D86348480A10 at Geo.Uni-Mainz.DE]
[cid:A982EFF3-D817-404A-BCD3-E74CE9CB9EE4 at Geo.Uni-Mainz.DE]
[cid:ED674059-4A40-445D-8D43-FE23792E4B9F at Geo.Uni-Mainz.DE]
[cid:0C431827-2AA7-47EF-A158-4004D3C0D9CE at Geo.Uni-Mainz.DE]
[cid:1C25FB42-709A-4D7B-BE52-94879AB5FAC7 at Geo.Uni-Mainz.DE]
On May 6, 2014, at 4:22 PM, Matthew Knepley wrote:
On Tue, May 6, 2014 at 7:23 AM, Püsök, Adina-Erika <puesoek at uni-mainz.de<mailto:puesoek at uni-mainz.de>> wrote:
Hello!
I was trying to implement some internal Dirichlet boundary conditions into an aij matrix of the form: A=( VV VP; PV PP ). The idea was to create an internal block (let's say Dirichlet block) that moves with constant velocity within the domain (i.e. check all the dofs within the block and set the values accordingly to the desired motion).
Ideally, this means to zero the rows and columns in VV, VP, PV corresponding to the dirichlet dofs and modify the corresponding rhs values. However, since we have submatrices and not a monolithic matrix A, we can choose to modify only VV and PV matrices.
The global indices of the velocity points within the Dirichlet block are contained in the arrays rowid_array.
What I want to point out is that the function MatZeroRowsColumnsIS() seems to create parallel artefacts, compared to MatZeroRowsIS() when run on more than 1 processor. Moreover, the results on 1 cpu are identical.
See below the results of the test (the Dirichlet block is outlined in white) and the piece of the code involved where the 1) - 2) parts are the only difference.
I am assuming that you are showing the result of solving the equations. It would be more useful, and presumably just as easy
to say:
a) Are the correct rows zeroed out?
b) Is the diagonal element correct?
c) Is the rhs value correct?
d) Are the columns zeroed correctly?
If we know where the problem is, its easier to fix. For example, if the rhs values are
correct and the rows are zeroed, then something is wrong with the solution procedure.
Since ZeroRows() works and ZeroRowsColumns() does not, this is a distinct possibility.
Thanks,
Matt
Thanks,
Adina Pusok
// Create an IS required by MatZeroRows()
ierr = ISCreateGeneral(PETSC_COMM_WORLD,numRowsx,rowidx_array,PETSC_COPY_VALUES,&isx); CHKERRQ(ierr);
ierr = ISCreateGeneral(PETSC_COMM_WORLD,numRowsy,rowidy_array,PETSC_COPY_VALUES,&isy); CHKERRQ(ierr);
ierr = ISCreateGeneral(PETSC_COMM_WORLD,numRowsz,rowidz_array,PETSC_COPY_VALUES,&isz); CHKERRQ(ierr);
1) /* ierr = MatZeroRowsColumnsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
ierr = MatZeroRowsColumnsIS(VV_MAT,isy,v_vv,x_push,rhs); CHKERRQ(ierr);
ierr = MatZeroRowsColumnsIS(VV_MAT,isz,v_vv,x_push,rhs); CHKERRQ(ierr);*/
2) ierr = MatZeroRowsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
ierr = MatZeroRowsIS(VV_MAT,isy,v_vv,x_push,rhs); CHKERRQ(ierr);
ierr = MatZeroRowsIS(VV_MAT,isz,v_vv,x_push,rhs); CHKERRQ(ierr);
ierr = MatZeroRowsIS(VP_MAT,isx,v_vp,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
ierr = MatZeroRowsIS(VP_MAT,isy,v_vp,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
ierr = MatZeroRowsIS(VP_MAT,isz,v_vp,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
ierr = ISDestroy(&isx); CHKERRQ(ierr);
ierr = ISDestroy(&isy); CHKERRQ(ierr);
ierr = ISDestroy(&isz); CHKERRQ(ierr);
Results (velocity) with MatZeroRowsColumnsIS().
1cpu<r01_1cpu_rows_columns.png> 4cpu<r01_rows_columns.png>
Results (velocity) with MatZeroRowsIS():
1cpu<r01_1cpu_rows.png> 4cpu<r01_rows.png>
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140509/1f1c932e/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spy_zerorows_1cpu.png
Type: image/png
Size: 15916 bytes
Desc: spy_zerorows_1cpu.png
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140509/1f1c932e/attachment-0005.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spy_zerorows_4cpu.png
Type: image/png
Size: 17690 bytes
Desc: spy_zerorows_4cpu.png
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140509/1f1c932e/attachment-0006.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spy_zerorowscol_1cpu.png
Type: image/png
Size: 16300 bytes
Desc: spy_zerorowscol_1cpu.png
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140509/1f1c932e/attachment-0007.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spy_zerorowscol_4cpu.png
Type: image/png
Size: 18174 bytes
Desc: spy_zerorowscol_4cpu.png
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140509/1f1c932e/attachment-0008.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: residual_tol.png
Type: image/png
Size: 3577 bytes
Desc: residual_tol.png
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140509/1f1c932e/attachment-0009.png>
More information about the petsc-users
mailing list