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

Matthew Knepley knepley at gmail.com
Fri May 4 08:05:00 CDT 2012


On Fri, May 4, 2012 at 8:59 AM, TAY wee-beng <zonexo at gmail.com> wrote:

>
> Hi,
>
> Is there anything else I can try to get it working right?
>

The MatGetNullSpaceRemove() is missing.

   Matt


> Thanks!
>
>
> 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> 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> 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
>
>


-- 
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/20120504/da7c7fc7/attachment-0001.htm>


More information about the petsc-users mailing list