linear elliptic vector equation

[nicolas aunai]

Hi,

I've written a little code based on the exercice 29, to solve my
equation (u - mu*laplacian(u) = f) with one DOF and with dirichlet
boundary conditions.

my code is here :

http://nau.cetp.ipsl.fr/inertia1.c

It seems that the code is *very* slow (never waited the end...) when
the dmmg_nlevels is > 4 (with a 3x3 coarsest grid).
The refinement factor has the default value 2, and I want my fine grid
to be for example nx=513, ny=257 grid points.


I have some more questions :

1/ How can I change the linear solver and preconditionner *in* the
code, not with the command -ksp_type and -pc_type ?

2/ What and Where is actually used the linear solver and the
preconditionners ? Is it used at each grid level as a smoother ? Can
we change the pre/post smoothing parameters ?

3/ My code is written to work in parallel, can I use a direct solver
when I'm at the coarset grid (small matrix) ?

4/ When I add -ksp_monitor, do I see the residual of the smoothing
iterations at each grid level ? If no, what is it ?

5/ My code is just a test code with some analytical RHS. In the end,
the solver is supposed to be added in a simulation code where the RHS
will be a/ the result of some other calculation, b/ changing with time
(the system being solved at each time step). Is it relevant to look at
the timing with this test code since the RHS is not the one that is
supposed to be (i.e. is the RHS very important for convergence speed)
?

6/ When the solver will be added in the simulation, will it be much
more efficient to setup the solution vector with the previous solution
as an initial guess ?

7/ My goal is not to do full dirichlet conditions like I have done in
the test code above, but rather DA_XPERIODIC and dirichlet in the Y
direction. I have observed that with a refinement factor of 2, and a
coarsest grid (nx=3,ny=3), when dmmg_nlevels is set to 2 I don't have
a 5x5 grid but a 6x5 grid, why ? for dmmg_nlevels = 8 I have 384x257
grid instead of the 257^2 grid I was expected... what is the reason
for this ?

8/ With DA_XPERIODIC, am I supposed to use the stencil as if I were
within the domain when I'm on periodic "borders", i.e. is col[1].i =
i-1 is ok when i=0 (left border) ?


Ok that's it for the moment,
Thx


Nico



2009/8/19 Barry Smith <bsmith at mcs.anl.gov>:
>
> On Aug 19, 2009, at 2:04 AM, nicolas aunai wrote:
>
>> Hi,
>>
>> OK then, lets go for multigrid.
>>
>>
>> When I'm looking for example at exercice 32 :
>>
>> http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/ksp/ksp/examples/tutorials/ex32.c.html
>>
>> I realize that I don't understand some things, in particular about the
>> DMMGSetKSP function.
>>
>> I need to give the RHS function and the ComputeMatrix function.
>>
>> 1/ My operator is always the same during the simulation, I was
>> thinking to build the matrix once and for all at the begining of the
>> program. Why can't I do that ? Is it because at each multigrid level
>> the required matrix is different (size) and it can only be built with
>> the information the DA is giving ?
>
>    Yes. The idea is rather than you giving the matrix. You give the function
> that fills up the matrix.
>>
>> 2/ I don't really understand the prototype of the callback function
>> ComputeMatrix. The first argument is the DMMG object, ok. But the
>> second and third arguments are matrices called 'J' and 'jac'.  It
>> seems (at least in ex32) that 'J' is not used, what is the goal of
>> this matrix in general ?
>
>   Sometimes one uses say a second order discretization to define the system
> but wants a cheap preconditioner built from the first order discretization.
> Hence the two matrices. If you are not doing this then you pass the same
> matrix in both slots.
>
>>
>> Is 'jac' stands for jacobian ? Reading the ComputeMatrix() code, it
>> seems that the 'jac' matrix is built like a standard finite difference
>> matrix, with the stencil corresponding to the operator I have to
>> solve. if so, why is it called 'jacobian' ? (sorry it looks more a
>> math question than a petsc question...)
>
>   This is just brought over from the nonlinear solver code. It could just be
> called matrix here.
>
>   Barry
>
>>
>>
>> Another quick question :
>>
>> Why is there no documentation for the function 'DMMGGetr(dmmg)' which
>> seems to be part of the petsc library since I can't see the definition
>> of the function in the exercice file. I suppose it computes the
>> residual ?
>>
>>
>> I guess this is it.... for the moment.
>> Thx
>>
>> Nico
>>
>>
>>
>> 2009/8/19 Barry Smith <bsmith at mcs.anl.gov>:
>>>
>>> On Aug 18, 2009, at 11:42 AM, nicolas aunai wrote:
>>>
>>>> Hello,
>>>>
>>>> Yes it is a bit overkill indeed :-) I was looking at ex29 actually...
>>>> which is still linear and with DA's but scalar.
>>>>
>>>> But I'd rather like not use multigrids, however I have not seen an
>>>> example of a code that uses DA's and KSP solvers.
>>>>
>>>> Every example use multigrids.
>>>
>>>  Multigrid is the way to go for this problem, it will be faster than some
>>> generic preconditioned CG.
>>>
>>>>
>>>> When you have a multicomponent linear equation to solve on a
>>>> rectangular uniform grid such as mine, is it good to :
>>>>
>>>> 1/  use DA's to represent your grid
>>>>
>>>> 2/ create the matrix corresponding to your operator once and for all
>>>> (the operator neve changes).( In fact there are 2 matrices because two
>>>> cases possibles (Neumann/periodic or Dirichlet/Periodic, depending on
>>>> the component). but never minds)
>>>>
>>>>
>>>> 3/ Call KSPSolve() 3 times (for each vector component associated with
>>>> the
>>>> DA) ?
>>>>
>>>  Yes. You can do all this using the DMMG object.
>>>
>>>  Note that the DMMG object in your linear case just manages a KSP and the
>>> meshes for you. In the end it is a KSPSolve() that actually does the
>>> solve
>>> and you have full control over the solve including using -pc_type icc or
>>> something else that does not use multigrid. Though multigrid will be
>>> faster.
>>>
>>>
>>>  Barry
>>>
>>>
>>>>
>>>> Thx
>>>> Nico
>>>>
>>>>
>>>> 2009/8/18 Aron Ahmadia <aron.ahmadia at kaust.edu.sa>:
>>>>>
>>>>> ex19 comes to mind, though it's a bit overkill for what you're doing...
>>>>>
>>>>>
>>>>>
>>>>> http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/snes/examples/tutorials/ex19.c.html
>>>>>
>>>>> ex19 uses DAs and DMMG, which is kind of like a meta-DA for using
>>>>> multigrid-style solvers.  Both work well with structured grids.
>>>>>
>>>>> You may also want to review the PETSc manual section on DA.
>>>>>
>>>>> A
>>>>>
>>>>> On Tue, Aug 18, 2009 at 3:30 AM, nicolas aunai
>>>>> <nicolas.aunai at gmail.com>
>>>>> wrote:
>>>>>>
>>>>>> Dear all,
>>>>>>
>>>>>>
>>>>>> I have a simulation code (particle in cell) in which I have to solve a
>>>>>> linear vector equation 4 times per time step. A typical run consists
>>>>>> of 50 000 time steps. The grid is rectangular and uniform of size Lx,
>>>>>> Ly, with nx+1 and ny+1 points in the x and y direction respectively.,
>>>>>> (nx, ny) could be at max (1024,1024).
>>>>>>
>>>>>> The vector equation is the following :
>>>>>>
>>>>>> B(x,y) - a*Laplacian(B(x,y)) = S(x,y)
>>>>>>
>>>>>>
>>>>>> Where B and S are 2D vector fields with 3 components (Bx, By, Bz and
>>>>>> Sx, Sy, Sz), each depending on the x and y coordinates.
>>>>>>
>>>>>>
>>>>>> 'a' is a positive constant, smaller than one, typically a= 0.02
>>>>>>
>>>>>> Boundary conditions may depend on for which B component we are solving
>>>>>> the equation. On the x=cst borders, the boundary is always periodic,
>>>>>> no matter what component is solved, but on y=cst borders, the boundary
>>>>>> condition can be either Neumann or Dirichlet.
>>>>>>
>>>>>> So far, I have written a small code that creates a vector solution, a
>>>>>> vector RHS and the matrix operator, and solve a scalar equation of
>>>>>> this type. Solving my vector equation would then just call 3 times
>>>>>> this kind of code... but I believe there is another way to deal with
>>>>>> vector fields and linear systems with Petsc, using DAs no ? Could
>>>>>> someone explain me how solve this the proper way and/or show me some
>>>>>> code solving linear system with vector fields ? (is there an example
>>>>>> in the exercices that I would not have seen ?)
>>>>>>
>>>>>>
>>>>>>
>>>>>> Thanks a lot
>>>>>> Nico
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Aron Jamil Ahmadia
>>>>> Assistant Research Scientist
>>>>> King Abdullah University of Science and Technology
>>>>>
>>>
>>>
>
>