linear elliptic vector equation

Mon Aug 24 01:24:32 CDT 2009

On 8/20/09, nicolas aunai <nicolas.aunai at gmail.com> wrote:
>
> Hi,

I will answer the questions that I can below:

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).

Send the output of -log_summary.

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 ?

You can use DMMGGetKSP() and then KSPGetPC(), and then extract
the various parts of the MG PC. However, I would never ever do this. First,
it precludes nice testing you can do by just setting the PC to LU and the
KSP to preonly. Furthermore, it seems to preclude nice optimizations to
be made later. I really really really recommend that you do not do this.

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 ?

Yes, and yes. You can always use -help to see all the options that can
customize the solve.

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

Yes. The default is to solve the coarse problem redundantly using LU.
You can also use a parallel 3rd party LU.

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

I believe you can use -dmmg_ksp_monitor, however I can't check right
now since I am at Domodedovo airport. You can check using -help.

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)
> ?

This depends completely on the particular problem. Usually not.

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 ?

Yes, usually if you are solving for the full solution and not an update.

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) ?

Yes.

   Matt

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
> >>>>>
> >>>
> >>>
> >
> >
>

-- 
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/20090824/114608ab/attachment.htm>