Community
cancel
Showing results for 
Search instead for 
Did you mean: 
88 Views

multiple calls to poisson solver produces odd results

Jump to solution
Greetings,
I'm have difficulties with the 3D Poisson solver. I'm trying to implement it in a particle-in-cell code that requires Poisson's equation to be solved repeatedly on a constant grid with changes to the boundary conditions (possibly) and charge densities on the internal mesh points.
According to the MKL guide, in the Sequence of Invoking PL Routines section, I should be able to initialize in the beginning, then loop the commit and solve routines, and clean up after the I'm done.
If I only perform one iteration, the correct solution is produced. However, if I try to loop, without changing any parameters, and look at the solution afterwards it is different if nx=ny=nz is not satisfied. If nx=ny=nz IS satisfied, then a given number of iterations produces the same answer (as expected). The code is given below. If anyone notices any problems, please let me know.

ThreeDCartPotSolver solver(64,64,64,1,1,1); //initialize solver class

solver.setY0BoundPotsConst(10); //set boundary points constant (all dirichlet)

solver.setY1BoundPotsConst(10);

solver.setX0BoundPotsConst(-10);

solver.setX1BoundPotsConst(-10);

for(int s=0; s<100; s++) //loop the solve routine (without altering boundary conditions).

{

solver.zeroSolutionVec(); //writes zeroes solution vector (f from MKL manual)

solver.commitSolver(); //commit

solver.calcPotential(); //solve

}

solver.writePotData("cartpot.dat"); //writes out f to data file

0 Kudos
1 Solution
88 Views
Well, I've reproduced your issue on MKL example by calling commit and solve step in loop. This issue could be resolved by 2 ways:
First variant:
loop
{
set f = f_old;
dpar[0]=(bx-ax);
dpar[1]=(by-ay);
dpar[2]=(bz-az);
commit();
solve();
}
Second variant:
loop
{
set f = f_old;
commit();
solve();
ipar[0]=9; // more detailed inIntel Math Kernel Library ReferenceManual
}
If these ways doesn't resolve your issue please wrote it here and I will try to find solution of problem on your example
With best regards,
Alexander Kalinkin

View solution in original post

8 Replies
88 Views
Hi,
As for me you have some problems with changed array ipar but I can't check it in your interfaces. Could you provide testcase that I could compiler and run on my side?
With best regards,
Alexander Kalinkin
88 Views
Hi Alexander,
I suspected that it might be something like that. I'll attach the files if you want to take a look at it. Thanks for your help.
Ryan
88 Views
Also check dpar[0]-dpar[2]. These values have different meaning before and after commit step, before commit step they store length of correspondent domain edge but after they store mesh size. So, if you call in loop {commit,solve} check that these parameters set as length of domain edge before each commit step.
With best regards,
Alexander Kalinkin
89 Views
Well, I've reproduced your issue on MKL example by calling commit and solve step in loop. This issue could be resolved by 2 ways:
First variant:
loop
{
set f = f_old;
dpar[0]=(bx-ax);
dpar[1]=(by-ay);
dpar[2]=(bz-az);
commit();
solve();
}
Second variant:
loop
{
set f = f_old;
commit();
solve();
ipar[0]=9; // more detailed inIntel Math Kernel Library ReferenceManual
}
If these ways doesn't resolve your issue please wrote it here and I will try to find solution of problem on your example
With best regards,
Alexander Kalinkin

View solution in original post

88 Views
Hi Alexander,
Thanks for your help. I chose the first variant, and it worked perfectly. That probably would have taken me quite some time to find on my own. Thanks, again!
Ryan
88 Views
Hi Ryan,
You are welcome, feel free to ask any question about Poisson library or other parts of MKL.
With best regards,
Alexander Kalinkin
88 Views
Hi Alexander,
I just wanted to update this thread with another issue I've found by looping over commit and solve. I noticed when I ran this loop thousands of times a memory leak developed. It disappeared when I changed the cycle to initialize->commit->solve->free. I wanted to avoid initializing at every iteration as I didn't know how much overhead was involved in this process.
The "Typical Order of Invoking PL Routines" figure in the MKL manual seems to indicate that initializing at the beginning of every iteration isn't required, so maybe this memory leak is a bug?
Thanks,
Ryan
88 Views

Hi Rian,

You are right, if you choose first variant (dpar[0]=bx-ax and so on) the memory leak appeared because it is not recommended way of use Poisson library in loop, it is just a trick. Please choose second variant (set ipar[0]=9) then there is not memory leaks and result of PL is correct.

With best regards,

Alexander Kalinkin

Reply