Solving a Time-Dependent 3D Poisson Equation Solver Using Intel MKL

[AhmadOneAPI]

I want to solve a time-dependent 3D Poisson equation using the Intel MKL d_Helmholtz_3D function.

\nabla^2\phi(x,y,z,t)=f(x,y,z,t)\;\;\;(x,y,z)\in[0, L_x]\times[0, L_y]\times[0, L_z]. See the attached picture.


The domain is a fixed domain, boundary conditions are the same, and the only term that changes is the right-hand side of the Poisson equation f(x,y,z), hence solving for different phi(x,y,z) at different times, t0, t1, ..., tn.

I have successfully solved the Poisson equation for f(x,y,z) at a specific time step and I can loop through all the time steps and call all 4 functions below and solve for the solution at each timestep. However, I want to see if there is room for optimization and acceleration of this computation.

(1) d_init_Helmholtz_3D(&ax, &bx, &ay, &by, &az, &bz, &nx, &ny, &nz, BCtype, &q, ipar, dpar, &stat);

(2) d_commit_Helmholtz_3D(f, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, &xhandle, &yhandle, ipar, dpar, &stat);

(3) d_Helmholtz_3D(f, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, &xhandle, &yhandle, ipar, dpar, &stat);

(4) free_Helmholtz_3D(&xhandle, &yhandle, ipar, &stat);

Specifically, I wanted to know if it is possible to factor out the common tasks among solving the Poisson equation for all of f(x,y,z,t), calling the common tasks once, and then loop through f(x,y,z,t_i) i = 0, ...N for solving for phi?

For example, can I call d_init_Helmholtz_3D once and then loop through d_commit_Helmholtz_3D and d_Helmholtz_3D as in below? Or, is there a better way?

// One potential optimization:

// Setup the solver once, for all f(x,y,z,t_i)

d_init_Helmholtz_3D(&ax, &bx, &ay, &by, &az, &bz, &nx, &ny, &nz, BCtype, &q, ipar, dpar, &stat);

// Loop through the RHS functions

for(int i = 0; i < N; i++){

    // f_i = f(x,y,z,t_i)

    d_commit_Helmholtz_3D(f_i, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, &xhandle, &yhandle, ipar, dpar, &stat);

    d_Helmholtz_3D(f_i, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, &xhandle, &yhandle, ipar, dpar, &stat);

    // save the solution to a different variable

    // e.g., sol_i = f_i;

}

// Finally 

free_Helmholtz_3D(&xhandle, &yhandle, ipar, &stat);

As an additional acceleration and optimization, are these functions thread-safe and can I solve the equation at different times in parallel? 

Also, is there a code sample for a time-dependent Poisson equation solver?

I would appreciate any tips and help.

Best regards.

[ShanmukhS_Intel]

Hi Ahmad,

Thanks for posting in Intel communities.

Regarding the code sample, Could you please let us know if the example in the below path meets your requirement and get back to us if there is any difference in our understanding?

C:\Program Files(x86)\Intel\oneAPI\mkl\latest\examples\examples_core_c.zip\c\pdepoisson\source\d_Helmholtz_3D_c.c

Are these functions thread-safe and can I solve the equation at different times in parallel? 

>>Intel MKL functions are generally thread-safe, and you can often control the level of parallelism through configuration options like setting the number of threads. However, as the threading behavior can vary depending on the function, you could use the thread control functions like mkl_set_num_threads() or mkl_get_max_threads()  or an internal threading layer to manage parallelism.

Best Regards,

Shanmukh.SS

[ShanmukhS_Intel]

Hi Ahmad,

A gentle reminder:

Could you please let us know if the example provided in the earlier reply helped?

Best Regards,

Shanmukh.SS

[AhmadOneAPI]

Hi Shanmukh.SS,

I did find an example but it is not helpful. I was looking for an example for time dependent Poisson solver.

My question in my post is mainly for a case where there are multiple right hand functions for the same domain and same boundary conditions for all those right hand side function.

Thanks.

[ShanmukhS_Intel]

Hi Ahmad,

Thanks for your reply.

For the use case provide, You could try running d_init_Helmholtz_3D for initialization once and free_Helmholtz_3D once to free the resource while looping through the d_commit_Helmholtz_3D and d_Helmholtz_3D for different values as you mentioned by saving the result in a different variable. Please get back to us if you face any errors after implementing this.

Please refer to the below URL which helps you with more information regarding the sequence of Invoking Poisson Solver Routines

https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2023-1/sequence-of-invoking-poisson-solver-routines.html

Based on this, if you still face any challenges, we would like to assist you further on that by checking with the appropriate team internally.

Best Regards,

Shanmukh.SS

[ShanmukhS_Intel]

Hi Ahmad,

A gentle reminder:

Could you please let us know if there are any updates on your issue? It helps us in guiding you further regarding your issue.

Best Regards,

Shanmukh.SS