[petsc-users] Problem with fortran version of ex29 in ksp

TAY wee-beng zonexo at gmail.com
Wed May 2 15:42:39 CDT 2012


On 2/5/2012 10:11 PM, Matthew Knepley wrote:
> On Wed, May 2, 2012 at 1:55 PM, TAY wee-beng <zonexo at gmail.com 
> <mailto:zonexo at gmail.com>> wrote:
>
>     Hi,
>
>     I did a MatView and VecView on both C and Fortran, right after Mat
>     and Vec assembly. I have attached the printout below. They are
>     exactly the same, but yet the result is different in Neumann
>     condition. However, the dirichlet condition gives the correct ans.
>     Is there anything else that could be wrong even if the Mat and Vec
>     are the same?
>
>
> Did you set the null space for the matrix when you have Neumann 
> conditions?
Yes, for the matrix, I set as:

call 
MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)

     call MatSetNullSpace(jac,nullspace,ierr)

     call MatNullSpaceDestroy(nullspace,ierr)

for the Vec,

call 
MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)

*!call MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)*

     call MatNullSpaceDestroy(nullspace,ierr)

MatNullSpaceRemove was comment out because there's error during linking


>
>    Matt
>
>     Thanks!
>
>     Fortran:
>
>     Matrix Object: 1 MPI processes
>       type: seqaij
>     row 0: (0, 2)  (1, -1)  (3, -1)
>     row 1: (0, -1)  (1, 3)  (2, -1)  (4, -1)
>     row 2: (1, -1)  (2, 2)  (5, -1)
>     row 3: (0, -1)  (3, 3)  (4, -1)  (6, -1)
>     row 4: (1, -1)  (3, -1)  (4, 4)  (5, -1)  (7, -1)
>     row 5: (2, -1)  (4, -1)  (5, 3)  (8, -1)
>     row 6: (3, -1)  (6, 2)  (7, -1)
>     row 7: (4, -1)  (6, -1)  (7, 3)  (8, -1)
>     row 8: (5, -1)  (7, -1)  (8, 2)
>     Vector Object:Vec_0000000084000000_0 1 MPI processes
>       type: mpi
>     Process [0]
>     0.25
>     0.0205213
>     1.135e-005
>     0.0205213
>     0.00168449
>     9.31663e-007
>     1.135e-005
>     9.31663e-007
>     5.15289e-010
>     Vector Object:Vec_0000000084000000_1 1 MPI processes
>       type: mpi
>     Process [0]
>     0.14924
>     0.0242397
>     -0.0260347
>     0.0242397
>     -0.0256192
>     -0.0400102
>     -0.0260347
>     -0.0400102
>     -0.0400102
>     Press any key to continue . . .
>
>     C:
>
>     Matrix Object: 1 MPI processes
>       type: seqaij
>     row 0: (0, 2)  (1, -1)  (3, -1)
>     row 1: (0, -1)  (1, 3)  (2, -1)  (4, -1)
>     row 2: (1, -1)  (2, 2)  (5, -1)
>     row 3: (0, -1)  (3, 3)  (4, -1)  (6, -1)
>     row 4: (1, -1)  (3, -1)  (4, 4)  (5, -1)  (7, -1)
>     row 5: (2, -1)  (4, -1)  (5, 3)  (8, -1)
>     row 6: (3, -1)  (6, 2)  (7, -1)
>     row 7: (4, -1)  (6, -1)  (7, 3)  (8, -1)
>     row 8: (5, -1)  (7, -1)  (8, 2)
>     Vector Object:Vec_0x1d3b000_0 1 MPI processes
>       type: mpi
>     Process [0]
>     0.25
>     0.0205212
>     1.135e-05
>     0.0205212
>     0.00168449
>     9.31663e-07
>     1.135e-05
>     9.31663e-07
>     5.15288e-10
>     Vector Object:Vec_0x1d3b000_1 1 MPI processes
>       type: mpi
>     Process [0]
>     0.139311
>     0.0305751
>     -0.0220633
>     0.0305751
>     -0.0135158
>     -0.042185
>     -0.0220633
>     -0.042185
>     -0.058449
>
>
>
>     Yours sincerely,
>
>     TAY wee-beng
>
>
>     On 1/5/2012 11:54 PM, Matthew Knepley wrote:
>>     On Tue, May 1, 2012 at 5:48 PM, TAY wee-beng <zonexo at gmail.com
>>     <mailto:zonexo at gmail.com>> wrote:
>>
>>         Hi,
>>
>>         Do you mean my method is wrong?
>>
>>         I am following the template of ex22f,
>>
>>         where the variables are declared as :
>>
>>         PetscScalar  v(5)
>>
>>         MatStencil   row(4),col(4,5)
>>
>>         Hence,
>>
>>         for the neumann BC
>>
>>         num = 1
>>
>>                        if (j/=0) then
>>
>>                            v(num) = -rho*HxdHy
>>
>>                            col(MatStencil_i,num) = i
>>
>>                            col(MatStencil_j,num) = j-1
>>
>>                            num = num + 1
>>
>>                        end if
>>
>>                        if (i/=0) then
>>
>>                            v(num) = -rho*HydHx
>>
>>                            col(MatStencil_i,num) = i-1
>>
>>                            col(MatStencil_j,num) = j
>>
>>                            num = num + 1
>>
>>                        end if
>>
>>                        if (i/=mx-1) then
>>
>>                            v(num) = -rho*HydHx
>>
>>                            col(MatStencil_i,num) = i+1
>>
>>                            col(MatStencil_j,num) = j
>>
>>                            num = num + 1
>>
>>                        end if
>>
>>                        if (j/=my-1) then
>>
>>                            v(num) = -rho*HxdHy
>>
>>                            col(MatStencil_i,num) = i
>>
>>                            col(MatStencil_j,num) = j+1
>>
>>                            num = num + 1
>>
>>                        end if
>>
>>                        v(num) = ((num-1)/2.0)*rho*(HxdHy + HydHx)
>>
>>                        print *, v
>>
>>                        col(MatStencil_i,num) = i
>>
>>                        col(MatStencil_j,num) = j
>>
>>                        !num = num + 1
>>
>>                        call
>>         MatSetValuesStencil(jac,i1,row,num,col,v,INSERT_VALUES,ierr)
>>
>>         I do not get any more out of range error. However,my ans is
>>         still different from that of ex29 in C.
>>
>>
>>     This is very simple. You have an error in your code. Checking it
>>     is very simple: run the code and
>>     break in MatSetValues(). Make sure ex29 makes calls with exactly
>>     the same indices as your ex29f.
>>
>>         Matt
>>
>>         Yours sincerely,
>>
>>         TAY wee-beng
>>
>>
>>     -- 
>>     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
>
>
>
>
> -- 
> 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/20120502/a845553f/attachment.htm>


More information about the petsc-users mailing list