[petsc-users] petsc-3.5.2: ex12 with Neumann BC

Thu Nov 27 07:44:29 CST 2014

On Thu, Nov 27, 2014 at 7:06 AM, Olivier Bonnefon <
olivier.bonnefon at avignon.inra.fr> wrote:

>  Hello,
>
> Thanks you, it works. If I remove the NullSpaceCondition I get the correct
> solution also with Neumann BC.
> From my point of view, the linear system I get is not singular because
> there is a unique solution. Do you agree ?
>

Yes, the non-uniqueness is for the Poisson prblem only.

  Thanks,

     Matt

> Regards,
> Olivier Bonnefon
>
>
>
> On 11/25/2014 03:01 PM, Matthew Knepley wrote:
>
>  On Tue, Nov 25, 2014 at 4:08 AM, Olivier Bonnefon <
> olivier.bonnefon at avignon.inra.fr> wrote:
>
>>  Hello,
>>
>> Thank you for your answer, I agree with you.
>>
>> I have adapted the ex12 example for the system:
>>
>> -\nabla \nabla u + u + f = 0 in \Omega  (1)
>>
>> with f = 4-x^2-y^2
>> The solution is still x^2-y^2.
>>
>> It consists in adding the following gradient function:
>>
>> void g0_uu(const PetscScalar u[], const PetscScalar u_t[], const
>> PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const
>> PetscScalar a_x[], const PetscReal x[], PetscScalar g0[])
>> {
>>   g0[0] = 1.0;
>> }
>> ...
>> //and set gradient
>> ierr = PetscDSSetJacobian(prob, 0, 0, g0_uu, NULL, NULL,
>> g3_uu);CHKERRQ(ierr);
>> ...
>> Of course, I have modified the f0_u function.
>>
>>
>> With Dirichlet BC, I get the solution.
>>
>> With Neumann BC, the solution is still shifted down (-2/3). The problem
>> (1) with Neumann BC has a unique solution.
>> How remove the condition \int_\Omega u = 0 ?
>>
>
>  That condition is one way of imposing boundary conditions for pressure.
> When I call MatSetNullSpace(), I am implying
> this condition. You can take out that call, but the system will be
> singular. You would have to add a condition on the
> pressure somehow, like fixing it at a point. This tends to make the system
> more ill-conditioned, but if that is appropriate
> for you, then you make a similar call to
>
>        ierr = DMPlexAddBoundary(cdm, user->bcType == DIRICHLET, "wall",
> user->bcType == NEUMANN ? "boundary" : "marker", 0, user->exactFuncs[0], 1,
> &id, user);CHKERRQ(ierr);
>
>  but for field 1. To test , you can jsut change 0 to 1 in this call.
> However, this sets the entire boundary, and you really only
> need to set the pressure at one point, so you could make another label
> that only has a single point and use that instead.
>
>    Thanks,
>
>       Matt
>
>
>>  Thanks
>> Olivier Bonnefon
>>
>>
>>
>> On 11/06/2014 07:14 AM, Matthew Knepley wrote:
>>
>>  On Mon, Nov 3, 2014 at 9:46 AM, Olivier Bonnefon <
>> olivier.bonnefon at avignon.inra.fr> wrote:
>>
>>>  Hello,
>>>
>>> Thank for your answer, I'll explain my trouble:
>>>
>>> My problem is that the BC Neumann leads to wrong result.
>>>
>>> With Dirichlet BC, I get:
>>> -------------------------------
>>>
>>> > ./ex12 -run_type full -refinement_limit 0.0    -bc_type dirichlet
>>> -interpolate 0 -petscspace_order 1 -show_initial -show_solution
>>> -dm_plex_print_fem 1
>>> ...
>>> ...
>>> Solution
>>> Vec Object:potential 1 MPI processes
>>>   type: seq
>>> 0.5
>>>
>>> This result is correct.
>>>
>>> With Neuman BC, I get:
>>> --------------------------------
>>>
>>> > ./ex12 -run_type full -refinement_limit 0.0    -bc_type neumann
>>> -interpolate 1 -petscspace_order 1 -show_initial -dm_plex_print_fem 1
>>> -show_solution -bd_petscspace_order 1 -snes_linesearch_monitor
>>> -snes_monitor -ksp_monitor_true_residual -snes_converged_reason
>>> -ksp_converged_reason
>>> ....
>>> ....
>>>
>>>
>>> Solution
>>> Vec Object:potential 1 MPI processes
>>>   type: seq
>>> -0.75
>>> -0.583333
>>> 0.0833333
>>> -0.583333
>>> -0.333333
>>> 0.416667
>>> 0.0833333
>>> 0.416667
>>> 1.25
>>>
>>>
>>> That is not the values of the solution x*x+y*y.
>>>
>>
>>  Sorry, this is poor documentation on my part. I will fix it in the
>> example. The Naumann solution
>> explicitly discards the constant part, meaning
>>
>>    \int_\Omega u = 0
>>
>>  If we look at my solution
>>
>>    \int^1_0 dx \int^1_0 dy x^2 + y^2 = 2/3
>>
>>  However, this is for the continuum result, whereas we are enforcing
>> this in the discrete case.
>> Thus, what we really get is 0.75, since we are integrating linear patches
>> instead of quadratic
>> patches. After shifting by that, you still do not get exactly x^2 + y^2
>> since there is some
>> discretization error. If you run with P2 you should get the exact answer,
>> shifted down to eliminate
>> the DC component.
>>
>>    Thanks,
>>
>>       Matt
>>
>>
>>>
>>> I tried many ksp options.
>>> Moreover, the neumann BC with "-run_type full" is not cover in the list
>>> https://bitbucket.org/petsc/petsc/src/fced3c3f9e703542693913793d15321603e40fe6/config/builder.py?at=master#cl-257
>>>
>>> Do you know what is wrong ?
>>>
>>> Thanks,
>>>
>>> Olivier Bonnefon
>>>
>>> On 10/31/2014 06:50 PM, Matthew Knepley wrote:
>>>
>>>  On Fri, Oct 31, 2014 at 10:43 AM, Olivier Bonnefon <
>>> olivier.bonnefon at avignon.inra.fr> wrote:
>>>
>>>> Hello,
>>>>
>>>> I'm working on the snes/examples/tutorial/ex12 version 3.5.2.
>>>>
>>>> I didn't succed to run the simplest case with Neumann BC:
>>>>
>>>> ./ex12 -run_type full -refinement_limit 0.0    -bc_type neumann
>>>> -interpolate 1 -petscspace_order 1 -show_initial -dm_plex_print_fem 1
>>>> -show_solution -bd_petscspace_order 1 -snes_linesearch_monitor
>>>> -snes_monitor -ksp_monitor_true_residual -snes_converged_reason
>>>> -ksp_converged_reason
>>>>
>>>> This leads to dofs negatives values.
>>>>
>>>
>>>  I do not understand what you mean. Please always mail the full error
>>> message.
>>>
>>>
>>>> Do you know the options to get a correct result with Neumann BC ?
>>>>
>>>
>>>  There are some tests here:
>>>
>>>
>>> https://bitbucket.org/petsc/petsc/src/fced3c3f9e703542693913793d15321603e40fe6/config/builder.py?at=master#cl-257
>>>
>>>     Matt
>>>
>>>
>>>> Regards,
>>>> Olivier Bonnefon
>>>>
>>>>
>>>
>>>
>>>  --
>>> 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
>>>
>>>
>>>
>>> --
>>> Olivier Bonnefon
>>> INRA PACA-Avignon, Unité BioSP
>>> Tel: +33 (0)4 32 72 21 58
>>>
>>>
>>
>>
>>  --
>> 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
>>
>>
>>
>> --
>> Olivier Bonnefon
>> INRA PACA-Avignon, Unité BioSP
>> Tel: +33 (0)4 32 72 21 58
>>
>>
>
>
>  --
> 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
>
>
>
> --
> Olivier Bonnefon
> INRA PACA-Avignon, Unité BioSP
> Tel: +33 (0)4 32 72 21 58
>
>

-- 
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/20141127/c0fee150/attachment.html>