[petsc-users] FEM on 2D poisson equation

[Barry Smith]

  Olivier,

    I concur with Mark that the most likely problem is a wrong Jacobian. You can look at http://www.mcs.anl.gov/petsc/documentation/faq.html#newton  and it gives you a series of things to check. They also pretty much apply to your case, likely if you work through them you'll determine the problem.  

     If you continue to have difficulties do not hesitate to contact us again,

   Barry

On Jul 30, 2013, at 11:38 AM, Mark F. Adams <mfadams at lbl.gov> wrote:

> 
> On Jul 30, 2013, at 10:06 AM, Olivier Bonnefon <olivier.bonnefon at avignon.inra.fr> wrote:
> 
>> Hello,
>> 
>> I want to use PETSC for large problem about diffusive model in epidemiological field.
>> I'm using the slides 'Advanced PETSc Tutorial, Maison de la Simulation, Orsay, France, June 2013 (Matt)'.
>> 
>> My first step is to simulated a problem 0=-\nabla u + f(u), where f is a linear or non-linear function, with FEM. To do this I'm adapting the example ex12.c for the linear problem:
>> 
>> 0=-\nabla u + w*w*u
>> with the exact solution u(x,y)=exp(w*x)+exp(w*y). The Dirichlet boundary condition are defined from the exact solution, like in example ex12.
>> 
>> To do this, I change only the two following functions:
>> 
>> double WW=1.0;
>> void quadratic_u_2d(const PetscReal x[], PetscScalar *u)
>> {
>>   *u = exp(WW*(x[0])) +exp(WW*(x[1]))  ;
>> }
>> 
>> void f0_u(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
>> {
>>   const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
>>   PetscInt       comp;
>> 
>>   for (comp = 0; comp < Ncomp; ++comp) f0[comp] = WW*WW*u[comp] ;
>> }
>> 
>> 
>> The result is :
>> $ ./ex12 -refinement_limit 0.01 -snes_monitor_short -snes_converged_reason
>>   0 SNES Function norm 22.1518 
>>   1 SNES Function norm 0.364312 
>>   2 SNES Function norm 0.0165162 
>>   3 SNES Function norm 0.000792446 
>>   4 SNES Function norm 3.81143e-05 
>>   5 SNES Function norm 1.83353e-06 
>>   6 SNES Function norm 8.8206e-08 
>> Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 6
>> Number of SNES iterations = 6
>> L_2 Error: 0.00511
>> 
> 
> -ksp_monitor will show the linear solver residuals and these should match from one nonlinear iteration to the next if the problem is linear.  Your function and Jacobean might not be consistent (my guess).
> 
>> 
>> Something is wrong because it is a linear problem, why the snes didn't converge in one iteration ?
>> 
>> Thanks a lot.
>> 
>> Olivier B
>> 
>> On 07/18/2013 03:39 PM, Olivier Bonnefon wrote:
>>> On 07/18/2013 03:26 PM, Matthew Knepley wrote:
>>>> On Thu, Jul 18, 2013 at 8:17 AM, Olivier Bonnefon <olivier.bonnefon at avignon.inra.fr> wrote:
>>>> It is what I wanted, it works. 
>>>> If I well understand the code, ex12.h contains the P1 implementation. To simulate an other system, with time dependences for examples (du/dt), I have to adapt the plugin functions.
>>>> 
>>>> The way I would add time dependence is to convert this from a SNES example into a TS example. I can help you
>>>> do this since I want to start using TS by default. Does this sound reasonable?
>>> Yes, of course. My goal is to simulate diffusive equation with non linear sources, for example Lotka-Voltera competion. 
>>> 
>>> Olivier B
>>>> 
>>>>   Thanks,
>>>> 
>>>>      Matt
>>>>  
>>>> Thanks a lot.
>>>> 
>>>> Olivier B
>>>> 
>>>> On 07/18/2013 01:12 PM, Matthew Knepley wrote:
>>>>> On Thu, Jul 18, 2013 at 5:08 AM, Olivier Bonnefon <olivier.bonnefon at avignon.inra.fr> wrote:
>>>>> Hello,
>>>>> 
>>>>> I have a 2-d heat equation that I want to simulate with Finit Element Method, to do this, I'm looking for an example solving 2D poisson equation with FEM (DMDA or DMPlex). Is there an example like this ?
>>>>> 
>>>>> There is, but there it is still somewhat problematic. I use FIAT to generate the basis function tabulation,
>>>>> so you have to configure with
>>>>> 
>>>>>   --download-fiat --download-scientificpython --download-generator
>>>>> 
>>>>> and you need mesh generation and partitioning
>>>>> 
>>>>>   --download-triangle --download-chaco
>>>>> 
>>>>> and then you can run SNES ex12 using Builder (which will make the header file)
>>>>> 
>>>>>   python2.7 ./config/builder2.py check src/snes/examples/tutorials/ex12.c
>>>>> 
>>>>> Jed and I are working on an all C version of tabulation which would mean that you could bypass
>>>>> the Python code generation step. Once the header is generated for the element you want, then
>>>>> you can just run the example as normal.
>>>>> 
>>>>>     Matt
>>>>>  
>>>>> Thanks a lot.
>>>>> 
>>>>> Olivier Bonnefon
>>>>> 
>>>>> -- 
>>>>> 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
>>> 
>> 
>> 
>> -- 
>> Olivier Bonnefon
>> INRA PACA-Avignon, Unité BioSP
>> Tel: +33 (0)4 32 72 21 58
>> 
>