linear elliptic vector equation
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 ? 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 ? 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...) 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
cast PetscScalar to double
Hi I'm using a Petsc 3 installation compiled with support for complex numbers. However, in the code I have made frequent use of the stl library and the standard type double. This now gives a compilation error since PetscScalar can not be casted into double. Do you know of a quick workaround or do I have to replace every instance of double with PetscScalar? Regards, Fredrik Bengzon
cast PetscScalar to double
you can replace (double)(var) with one of the following [approprate to your code] thingy. PetscAbsScalar(var) PetscRealPart(var) Satish On Wed, 19 Aug 2009, Fredrik Bengzon wrote: Hi I'm using a Petsc 3 installation compiled with support for complex numbers. However, in the code I have made frequent use of the stl library and the standard type double. This now gives a compilation error since PetscScalar can not be casted into double. Do you know of a quick workaround or do I have to replace every instance of double with PetscScalar? Regards, Fredrik Bengzon
linear elliptic vector equation
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
