[petsc-users] petsc-3.5.2: ex12 with Neumann BC
Olivier Bonnefon
olivier.bonnefon at avignon.inra.fr
Thu Nov 27 07:06:58 CST 2014
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 ?
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
> <mailto: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
>> <mailto: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
>>> <mailto: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 <tel:%2B33%20%280%294%2032%2072%2021%2058>
>>
>>
>>
>>
>> --
>> 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 <tel:%2B33%20%280%294%2032%2072%2021%2058>
>
>
>
>
> --
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141127/9f1e4505/attachment-0001.html>
More information about the petsc-users
mailing list