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
<[email protected]
<mailto:[email protected]>> 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
<[email protected]
<mailto:[email protected]>> 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
<[email protected]
<mailto:[email protected]>> 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
