<div><span class="gmail_quote">On 8/20/09, <b class="gmail_sendername">nicolas aunai</b> <<a href="mailto:nicolas.aunai@gmail.com">nicolas.aunai@gmail.com</a>> wrote:</span>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">Hi,</blockquote>
<div> </div>
<div> </div>
<div>I will answer the questions that I can below:</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">I've written a little code based on the exercice 29, to solve my<br>equation (u - mu*laplacian(u) = f) with one DOF and with dirichlet<br>
boundary conditions.<br><br>my code is here :<br><br><a href="http://nau.cetp.ipsl.fr/inertia1.c">http://nau.cetp.ipsl.fr/inertia1.c</a><br><br>It seems that the code is *very* slow (never waited the end...) when<br>the dmmg_nlevels is > 4 (with a 3x3 coarsest grid).</blockquote>
<div> </div>
<div> </div>
<div>Send the output of -log_summary.</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">The refinement factor has the default value 2, and I want my fine grid<br>to be for example nx=513, ny=257 grid points.<br>
<br><br>I have some more questions :<br><br>1/ How can I change the linear solver and preconditionner *in* the<br>code, not with the command -ksp_type and -pc_type ?</blockquote>
<div> </div>
<div>You can use DMMGGetKSP() and then KSPGetPC(), and then extract</div>
<div>the various parts of the MG PC. However, I would never ever do this. First,</div>
<div>it precludes nice testing you can do by just setting the PC to LU and the</div>
<div>KSP to preonly. Furthermore, it seems to preclude nice optimizations to</div>
<div>be made later. I really really really recommend that you do not do this.</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">2/ What and Where is actually used the linear solver and the<br>preconditionners ? Is it used at each grid level as a smoother ? Can<br>
we change the pre/post smoothing parameters ?</blockquote>
<div> </div>
<div> </div>
<div>Yes, and yes. You can always use -help to see all the options that can</div>
<div>customize the solve.</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">3/ My code is written to work in parallel, can I use a direct solver<br>when I'm at the coarset grid (small matrix) ?</blockquote>
<div> </div>
<div> </div>
<div>Yes. The default is to solve the coarse problem redundantly using LU.</div>
<div>You can also use a parallel 3rd party LU.</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">4/ When I add -ksp_monitor, do I see the residual of the smoothing<br>iterations at each grid level ? If no, what is it ?</blockquote>
<div> </div>
<div>I believe you can use -dmmg_ksp_monitor, however I can't check right</div>
<div>now since I am at Domodedovo airport. You can check using -help.</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">5/ My code is just a test code with some analytical RHS. In the end,<br>the solver is supposed to be added in a simulation code where the RHS<br>
will be a/ the result of some other calculation, b/ changing with time<br>(the system being solved at each time step). Is it relevant to look at<br>the timing with this test code since the RHS is not the one that is<br>supposed to be (i.e. is the RHS very important for convergence speed)<br>
?</blockquote>
<div> </div>
<div>This depends completely on the particular problem. Usually not.</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">6/ When the solver will be added in the simulation, will it be much<br>more efficient to setup the solution vector with the previous solution<br>
as an initial guess ?</blockquote>
<div> </div>
<div>Yes, usually if you are solving for the full solution and not an update.</div>
<div> </div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">7/ My goal is not to do full dirichlet conditions like I have done in<br>the test code above, but rather DA_XPERIODIC and dirichlet in the Y<br>
direction. I have observed that with a refinement factor of 2, and a<br>coarsest grid (nx=3,ny=3), when dmmg_nlevels is set to 2 I don't have<br>a 5x5 grid but a 6x5 grid, why ? for dmmg_nlevels = 8 I have 384x257<br>
grid instead of the 257^2 grid I was expected... what is the reason<br>for this ?<br><br>8/ With DA_XPERIODIC, am I supposed to use the stencil as if I were<br>within the domain when I'm on periodic "borders", i.e. is col[1].i =<br>
i-1 is ok when i=0 (left border) ?</blockquote>
<div> </div>
<div>Yes.</div>
<div> </div>
<div> Matt</div><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">Ok that's it for the moment,<br>Thx<br><br><br>Nico<br><br><br><br>2009/8/19 Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>>:<br>
><br>> On Aug 19, 2009, at 2:04 AM, nicolas aunai wrote:<br>><br>>> Hi,<br>>><br>>> OK then, lets go for multigrid.<br>>><br>>><br>>> When I'm looking for example at exercice 32 :<br>
>><br>>> <a href="http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/ksp/ksp/examples/tutorials/ex32.c.html">http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/ksp/ksp/examples/tutorials/ex32.c.html</a><br>
>><br>>> I realize that I don't understand some things, in particular about the<br>>> DMMGSetKSP function.<br>>><br>>> I need to give the RHS function and the ComputeMatrix function.<br>>><br>
>> 1/ My operator is always the same during the simulation, I was<br>>> thinking to build the matrix once and for all at the begining of the<br>>> program. Why can't I do that ? Is it because at each multigrid level<br>
>> the required matrix is different (size) and it can only be built with<br>>> the information the DA is giving ?<br>><br>> Yes. The idea is rather than you giving the matrix. You give the function<br>
> that fills up the matrix.<br>>><br>>> 2/ I don't really understand the prototype of the callback function<br>>> ComputeMatrix. The first argument is the DMMG object, ok. But the<br>>> second and third arguments are matrices called 'J' and 'jac'. It<br>
>> seems (at least in ex32) that 'J' is not used, what is the goal of<br>>> this matrix in general ?<br>><br>> Sometimes one uses say a second order discretization to define the system<br>> but wants a cheap preconditioner built from the first order discretization.<br>
> Hence the two matrices. If you are not doing this then you pass the same<br>> matrix in both slots.<br>><br>>><br>>> Is 'jac' stands for jacobian ? Reading the ComputeMatrix() code, it<br>>> seems that the 'jac' matrix is built like a standard finite difference<br>
>> matrix, with the stencil corresponding to the operator I have to<br>>> solve. if so, why is it called 'jacobian' ? (sorry it looks more a<br>>> math question than a petsc question...)<br>><br>
> This is just brought over from the nonlinear solver code. It could just be<br>> called matrix here.<br>><br>> Barry<br>><br>>><br>>><br>>> Another quick question :<br>>><br>>> Why is there no documentation for the function 'DMMGGetr(dmmg)' which<br>
>> seems to be part of the petsc library since I can't see the definition<br>>> of the function in the exercice file. I suppose it computes the<br>>> residual ?<br>>><br>>><br>>> I guess this is it.... for the moment.<br>
>> Thx<br>>><br>>> Nico<br>>><br>>><br>>><br>>> 2009/8/19 Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>>:<br>>>><br>>>> On Aug 18, 2009, at 11:42 AM, nicolas aunai wrote:<br>
>>><br>>>>> Hello,<br>>>>><br>>>>> Yes it is a bit overkill indeed :-) I was looking at ex29 actually...<br>>>>> which is still linear and with DA's but scalar.<br>
>>>><br>>>>> But I'd rather like not use multigrids, however I have not seen an<br>>>>> example of a code that uses DA's and KSP solvers.<br>>>>><br>>>>> Every example use multigrids.<br>
>>><br>>>> Multigrid is the way to go for this problem, it will be faster than some<br>>>> generic preconditioned CG.<br>>>><br>>>>><br>>>>> When you have a multicomponent linear equation to solve on a<br>
>>>> rectangular uniform grid such as mine, is it good to :<br>>>>><br>>>>> 1/ use DA's to represent your grid<br>>>>><br>>>>> 2/ create the matrix corresponding to your operator once and for all<br>
>>>> (the operator neve changes).( In fact there are 2 matrices because two<br>>>>> cases possibles (Neumann/periodic or Dirichlet/Periodic, depending on<br>>>>> the component). but never minds)<br>
>>>><br>>>>><br>>>>> 3/ Call KSPSolve() 3 times (for each vector component associated with<br>>>>> the<br>>>>> DA) ?<br>>>>><br>>>> Yes. You can do all this using the DMMG object.<br>
>>><br>>>> Note that the DMMG object in your linear case just manages a KSP and the<br>>>> meshes for you. In the end it is a KSPSolve() that actually does the<br>>>> solve<br>>>> and you have full control over the solve including using -pc_type icc or<br>
>>> something else that does not use multigrid. Though multigrid will be<br>>>> faster.<br>>>><br>>>><br>>>> Barry<br>>>><br>>>><br>>>>><br>>>>> Thx<br>
>>>> Nico<br>>>>><br>>>>><br>>>>> 2009/8/18 Aron Ahmadia <<a href="mailto:aron.ahmadia@kaust.edu.sa">aron.ahmadia@kaust.edu.sa</a>>:<br>>>>>><br>>>>>> ex19 comes to mind, though it's a bit overkill for what you're doing...<br>
>>>>><br>>>>>><br>>>>>><br>>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/snes/examples/tutorials/ex19.c.html">http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/snes/examples/tutorials/ex19.c.html</a><br>
>>>>><br>>>>>> ex19 uses DAs and DMMG, which is kind of like a meta-DA for using<br>>>>>> multigrid-style solvers. Both work well with structured grids.<br>>>>>><br>
>>>>> You may also want to review the PETSc manual section on DA.<br>>>>>><br>>>>>> A<br>>>>>><br>>>>>> On Tue, Aug 18, 2009 at 3:30 AM, nicolas aunai<br>
>>>>> <<a href="mailto:nicolas.aunai@gmail.com">nicolas.aunai@gmail.com</a>><br>>>>>> wrote:<br>>>>>>><br>>>>>>> Dear all,<br>>>>>>><br>
>>>>>><br>>>>>>> I have a simulation code (particle in cell) in which I have to solve a<br>>>>>>> linear vector equation 4 times per time step. A typical run consists<br>
>>>>>> of 50 000 time steps. The grid is rectangular and uniform of size Lx,<br>>>>>>> Ly, with nx+1 and ny+1 points in the x and y direction respectively.,<br>>>>>>> (nx, ny) could be at max (1024,1024).<br>
>>>>>><br>>>>>>> The vector equation is the following :<br>>>>>>><br>>>>>>> B(x,y) - a*Laplacian(B(x,y)) = S(x,y)<br>>>>>>><br>>>>>>><br>
>>>>>> Where B and S are 2D vector fields with 3 components (Bx, By, Bz and<br>>>>>>> Sx, Sy, Sz), each depending on the x and y coordinates.<br>>>>>>><br>>>>>>><br>
>>>>>> 'a' is a positive constant, smaller than one, typically a= 0.02<br>>>>>>><br>>>>>>> Boundary conditions may depend on for which B component we are solving<br>
>>>>>> the equation. On the x=cst borders, the boundary is always periodic,<br>>>>>>> no matter what component is solved, but on y=cst borders, the boundary<br>>>>>>> condition can be either Neumann or Dirichlet.<br>
>>>>>><br>>>>>>> So far, I have written a small code that creates a vector solution, a<br>>>>>>> vector RHS and the matrix operator, and solve a scalar equation of<br>
>>>>>> this type. Solving my vector equation would then just call 3 times<br>>>>>>> this kind of code... but I believe there is another way to deal with<br>>>>>>> vector fields and linear systems with Petsc, using DAs no ? Could<br>
>>>>>> someone explain me how solve this the proper way and/or show me some<br>>>>>>> code solving linear system with vector fields ? (is there an example<br>>>>>>> in the exercices that I would not have seen ?)<br>
>>>>>><br>>>>>>><br>>>>>>><br>>>>>>> Thanks a lot<br>>>>>>> Nico<br>>>>>><br>>>>>><br>>>>>><br>
>>>>> --<br>>>>>> Aron Jamil Ahmadia<br>>>>>> Assistant Research Scientist<br>>>>>> King Abdullah University of Science and Technology<br>>>>>><br>
>>><br>>>><br>><br>><br></blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener