Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6956 Discussions

2D Poisson Intel MKL - Set Fixed Internal Nodes

JP12321
Beginner
745 Views

Hi All,

 

I've recently got the MKL 2D Poisson Solver up and running and I believe it's working properly for my simulations. 

 

What I would like to know, is there a way of setting fixed values for internal nodes in the domain? 

 

I am using the Poisson equation to calculate electrostatic potential, so am going from a charge density distribution to a potential array. 

 

I cannot disclose my source code as I am under an NDA, so I will try to explain my thoughts here:

 

1) All boundaries of my domain are = 0V (Dirichlet) 

2) There are two vertical lines where there is an applied voltage (Call it 1V)

 

So if I had a 9x9 array, the initial inputs for the potential array would look something like:

0      0      0      0      0      0      0      0      0  

0      0      0      0      0      0      0      0      0 

0      0      0      0      0      0      0      0      0 

0      0      0      1      0      1      0      0      0 

0      0      0      1      0      1      0      0      0 

0      0      0      1      0      1      0      0      0 

0      0      0      0      0      0      0      0      0 

0      0      0      0      0      0      0      0      0 

0      0      0      0      0      0      0      0      0 

So ignoring any charge distribution for now, I would expect the updated solution, just due to permittivity, to resemble (arbitrary values of course) something along the lines of this:

0      0      0      0      0      0      0      0      0  

0      0      0      0      0      0      0      0      0 

0      0      0    1/2    0    1/2     0     0      0 

0      0    1/2    1   1/2     1    1/2    0      0 

0      0    1/2    1   1/2     1    1/2    0      0 

0      0    1/2    1   1/2     1    1/2    0      0 

0      0      0    1/2    0    1/2     0     0      0 

0      0      0      0      0      0      0      0      0 

0      0      0      0      0      0      0      0      0 

 

As the potential fields will decrease in value the further from the fixed potential region they go. 

 

I would consider this the 'static' potential field - it comes from the fixed potential nodes and will be constant throughout the simulation. 

 

Additionally, there will be charged particles in the domain, giving a charge density, which through the Poisson solver will create their own potential fields. As the particles are moving around, I would consider this a 'dynamic' potential field. 

 

I believe I have modified the MKL 2D Poisson solver example properly, and it is now rapidly and accurately evaluating the 'dynamic' potential. However, I can't begin to think how to evaluate the 'static' field. Would anyone be able to advise on this? 

 

I have thought about things such as setting charge density to zero, etc, however I don't actually know how to prescribe fixed potential values into the solver. The only values I know how to set are the boundary conditions, which for my case are all zero. Therefore, if I have a zero charge density there will be no regions of potential at all! 

 

I have alternatively thought of simulating the potential fields in a separate FEM program and importing it in, however this seems quite cumbersome and I intend on performing many simulations with differing grid parameters. 

 

Does anyone know of a way of achieving this within the MKL? 

Once I have the 'static' field, I can just add it to the 'dynamic' field, giving me the total field array. 

 

Please ask any more questions if anything needs to be clarified - I do appreciate that not being able to include code is not very helpful! 

Many thanks in advance. 

 

(Probably worth mentioning that I'm using Fortran 90, but this may not make a difference)

0 Kudos
1 Solution
Khang_N_Intel
Employee
548 Views

Hi Joseff,


Since the solutions have been provided, I would like to close this ticket.

There will be no more dialogue on this thread.


Best regards,

Khang


View solution in original post

0 Kudos
7 Replies
VidyalathaB_Intel
Moderator
697 Views

Hi Joseff,


Thanks for reaching out to us.


>>is there a way of setting fixed values for internal nodes in the domain? 

I have forwarded this to the concerned team and they are looking into this issue.


Additionally, If you have active priority support, we would recommend you submit a feature request at https://supporttickets.intel.com/?lang=en-US


Please refer to the below link for more information

https://software.intel.com/content/www/us/en/develop/articles/how-to-create-a-support-request-at-online-service-center.html


Please do let us know if you do not have active priority support.


Regards,

Vidya.


0 Kudos
Khang_N_Intel
Employee
680 Views

Hi Joseff,


We will take a look and will let you know how we would do nexrt.


Best,

Khang


0 Kudos
JP12321
Beginner
660 Views

Hi Both,

 

Thanks for your help. I'm still actively trying to work this out, however am still not having much luck.

I personally do not have active priority support, however my research institute may. Could you send me a PM?

Thanks

0 Kudos
Khang_N_Intel
Employee
646 Views

Hi Joseff,


That should be fine if you don't have paid support. We can work through the forum.


Best regards,

Khangf


0 Kudos
Khang_N_Intel
Employee
587 Views

Hi Joseff,


I apologize for the delay in response!

The mkl team will see what they can do. Currently, they are short-staffed.

I will let you know when I receive any news.


Best regards,

Khang


0 Kudos
Khang_N_Intel
Employee
583 Views

Hi Joseff,


This is the response from the mkl team.


Note that a mixture of \LaTeX symbols and English are being used:


Problem statement:


Solve for the electrical potential $\phi(\mathbf{x})$ for $\mathbf{x}\in\Omega$ (a bounded subset of $R^{2 or 3}$)

such that 


1. $-\Delta phi = f$ for $\mathbf{x}\in\Omega$

2. $\phi(\mathbf{x}) = 0$ for $\mathbf{x}\in\partial\Omega$  (the boundary of $\Omega$)

3. And subject to the internal constraint $\phi(\mathbf{x}) = V(\mathbf{x})$  for $\mathbf{x}\in\Gamma$ which is internal to $\Omega$.



This is a rather tricky problem to solve and if $\Gamma$ is too close to the boundaries, the boundary condition setting of $\phi=0$ may be incompatible… I will provide two suggestions for it, the first a more theoretical one, and may hopefully work well in this case. The second one is a "quick and dirty" solution which may lead to a wider field of solutions that I will just mention, as it represents a nice world of complicated but generally quite useful solutions to know and implement.



A theory based solution:


The theoretical solution is to use the superposition principle from electrical fields and essentially break your subdomain down into several sub problems. Taking $\Gamma$ consisting of two vertical strips of domain with the applied voltage = 1 as the author visually described in the problem statement, the contribution from each strip is independent and with the right breakdown of the domain, additive. (This is assuming that the effect of the voltage on one strip doesn’t pass through the fixed voltage on the other strip but can affect everything around it, and this assumption is definitely up for debate.)


So what is that breakdown?  


A. The first sub problem consists of the left domain over to the first strip. You want to solve the pde 1. On this sub problem with new boundary conditions of Dirichlet boundary conditions = 1 on the strip and normal derivative = 0 (neumann bc) for the degrees of freedom between the strip and the corners on the same side. Then Dirichlet boundary conditions = 0 for the remaining 3 sides.

B. The second sub problem is a little more complicated. It is a cut again down the left most strip of $\Gamma$ and then cuts from the top and bottom of the right strip over to the far most right edge. It looks like a rectangle with two prongs shooting off to the right. Now you will solve the pde on this domain with same combination of Dirichlet and neumann bc on the left side, and 0 Dirichlet bc on the rest of the sides.

C. The third and fourth are the same thing but using the second strip as base so by symmetry it is similar to A and B.


Get your 4 solutions and then add them together to get your final solution. Now, I don’t believe that the fast poisson solver in oneMKL can accept mixed boundary conditions like this on a rectangular domain so again you may need to discretize the pde yourself and solve using your favorite iterative or direct solver technique.


Side note: Now, if you want to remove the assumption we made of effect of voltage not passing through the other strip, then you can put the B. part as the rest of the domain left after A which is a much simpler convex rectangular system to solve… however when you go to add everything back up, you will find that the values of the electric field on your two strips is no longer = 1, but something larger… I am not sure how to go about trying to remove that issue


Quick and Dirty Solution:


The dirty solution is to set up a finite difference scheme for the poisson system in eq 1 and 2 without regard to the constraint 3. Then once you have the linear algebra system, since finite difference unknowns(degrees of freedom) are values at the nodes they represent in the domain, you can find the degrees of freedom that correspond to the nodes on Gamma and change the linear algebra system for those degrees of freedom (nodes) to be the identity * phi = 1 (or better yet, to keep scaling more consistent use a scale of h^2 or 1/n^2 on both rhs and in matrix, that is, in the linear algebra system remove the nonzeros of the row and just putting an h^2 or (1/n^2) on the main diagonal and target value times the same thing in the rhs spot. Now solve the system using your favorite iterative solver like (preconditioned conjugate gradient with SSOR preconditioner or anything more fancy) or with a direct solver like oneMKL Pardiso and you will get something that resembles your desired goal.


This is considered dirty since there is may not be a good stability of such linear systems and you have possibly broken some of the assumptions inherent in a finite difference approximation of a pde. 



Making the Dirty solution a little more clean but getting much more advanced:


The way to bring this type of constraint into a more stable approach is to delve into the world of constrained optimization, particularly introducing Lagrange multipliers into your pde system that help to enforce such an internal constraint. You end up with a pde with constraints that you will need to discretize and use some sort of saddle point iterative solver to alternately solve for solution and lagrange multipliers to enforce such constraints in a controlled and stable way. Note that the dirty solution is equivalent to pushing the Lagrange multipliers to infinity (or extremely large) to enforce the constraint on the solution. 


Please let me know if you have any questions.


Best regards,

Khang



0 Kudos
Khang_N_Intel
Employee
549 Views

Hi Joseff,


Since the solutions have been provided, I would like to close this ticket.

There will be no more dialogue on this thread.


Best regards,

Khang


0 Kudos
Reply