I have been writing some code in C for particle-in-cell simulation. One of the steps of the PIC algorithm requires to solve (numerically) Poisson's equation

$$ \Delta \varphi = - 4 \pi \rho. $$

Initially I want to limit the program to 2d case ( instead of 3d ) and use the finite difference method ( instead of finite elements ).

Much to my surprise, I was not able to find any free open source C library for this task ( i.e. to solve 2d Poisson's equation using the finite difference method ). The only thing I've found is a FISHPACK - collection of FORTRAN77 routines, part of a famous SLATEC library.

I'm fine with fortran, but a call to a fortran routine from C requires some unnecessary array transformations and "native" C solution seems more convenient. Another reason I hesitate to use FISHPACK is that it uses single precision floats with implicit type declaration. Using the gcc '-fdefault-real-8' flag I was able to convert it to double precision, but I'm not sure whether or not everything still works correctly after this.

So, my question is: are there any free ( as in freedom ) C/C++ libraries to solve 2d Poisson equation using the finite difference method? And since I do intend to switch to 3d and finite elements, advices on C library for this purpose are also welcome.

  • $\begingroup$ It looks like you want to solve for the gravity potential. This typically involves solving in all of ${\mathbb R}^3$, which presents a significant complication of your task. Is this so? $\endgroup$ Mar 22 '14 at 19:16
  • $\begingroup$ No, I want to do a plasma simulation where electromagnetic interaction is a primary force and gravity is typically neglected. The simulation domain is finite, not the whole ${\mathbb R}^3$. Besides, I want to start from 2D domain and later switch to 3D. $\endgroup$ Mar 23 '14 at 8:02
  • $\begingroup$ Is uniform orthogonal grid sufficient for you? In that case I am preparing to release such a library. Which is a single purpose thing that could be used instead of FISHPACK for people that do not want to use large frameworks. Fortran 2008 with C bindings... $\endgroup$
    – Vladimir F
    Mar 24 '14 at 11:44
  • $\begingroup$ Yes, uniform orthogonal grid is currently sufficient. Please, do post a link here when you make a release. $\endgroup$ Mar 29 '14 at 16:36
  • 2
    $\begingroup$ PoisFFT by @VladimirF. $\endgroup$ Aug 18 '15 at 21:08

2-D Poisson by finite differences looks to be a few examples in PETSc (KSP: ex29, ex32, ex50). I'd probably use a nonlinear solver (SNES) or time-stepper (TS) in PETSc anyway for the added flexibility. 2-D Poisson via finite differences is probably an example in other packages (for instance, probably Trilinos). I know it's an example given in Using MPI.

There are many free C/C++ libraries to solve Poisson via finite elements. You could try PETSc again, deal.II, DUNE, FEniCS, any other library in this list on Wikipedia, or (probably) Trilinos. Using a library like deal.II or FEniCS, it looks to be relatively easy to change the dimensionality of a code (i.e., go from 2-D to 3-D), so you might consider writing a 2-D finite element code to start, and then converting it to 3-D.

  • $\begingroup$ Yes, deal.II (dealii.org) has several examples solving the Poisson equation in 2d, 3d, with uniform or adaptive meshes, and also on a few 10,000 processors if you need to. $\endgroup$ Mar 22 '14 at 19:15
  • $\begingroup$ Thanks for the references to the examples. I'll take a closer look at PETSc. $\endgroup$ Mar 23 '14 at 7:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.