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

Tue Nov 25 08:01:41 CST 2014

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141125/499ae4a6/attachment-0001.html>