Fix issue #100: get_val with HAVE_MPI

[thchr]

Hi Steven,

I finally found sufficient incentive to attempt solving #100, i.e. making real-space interpolation work with MPI.
I pretty much went about it in the way we discussed there (after finally understanding, I think, the transposition of x/y indices that you mentioned).

I've compiled this with and without MPI and verified that the interpolation produces good agreement for a simple 3D test case (see below) run with mpb vs. mpb-mpi:

Value type	mpb (res=32)	mpb-mpi (res=32)	mpb (res=64)	mpb-mpi (res=64)
H-norm	0.78812	0.79027	0.83725	0.83490
H-energy	0.65716	0.66058	0.71445	0.71049
D-norm	6.40814	6.40644	6.47873	6.48145
D-energy	6.70220	6.69825	6.81435	6.82025

evaluated at an arbitrary point. The results for mpb-mpi remain unchanged (up to tiny convergence variations) when run with different number processes.

3D test example (click to expand)

(set! resolution 64)
(set! num-bands 3)
(set! output-epsilon (lambda () (print "skipping output-epsilon\n")))

; define a simple geometry with inversion
(set! geometry-lattice (make lattice (size 1 1 1) 
                                     (basis1 1 0 0) (basis2 0 1 0) (basis3 0 0 1)))
(set! geometry (list (make sphere
                        (center 0 0 0) (radius 0.25) 
                        (material (make dielectric (epsilon 13) )))))
; set a k-point
(set! k-points (list (vector3 0.3532 0.1512 0.175)))


(run) ; run the calculation

; compute some interpolating values at an arbitrary real-space point
(define r (vector3 0.0 0.12 -0.22))
(define band-idx 3)

(print "\n\n")
(get-hfield band-idx)
(print "|H|:      " (vector3-norm (get-field-point r)) "\n")
(compute-field-energy)
(print "H-energy: " (get-energy-point r) "\n")

(get-dfield band-idx)
(print "|D|:      " (vector3-norm (get-field-point r)) "\n")
(compute-field-energy)
(print "D-energy: " (get-energy-point r) "\n")

Does this seem like a correct implementation? One concern I have is I didn't really think about the case of HAVE_MPI and no SCALAR_COMPLEX: would that need additional special-casing?

I also snuck in a tiny addition to .gitignore for working in VSCode (can be dropped, of course).

[thchr]

The gentlest of bumps @stevengj: thoughts on this?

[stevengj]

We could actually call this in interp_val on the result of the bilinear interpolation (this is the only function that calls get_val), rather than in get_val, to save a factor of 8 in the number of allreduce calls.

[thchr]

Ah, nice. Just to check my understanding, you mean we could instead reduce over the operations here,

mpb/mpb/fields.c

Lines 628 to 631 in 0d886cf

	return(((D(x,y,z)*(1.0-dx) + D(x2,y,z)*dx) * (1.0-dy) +
	(D(x,y2,z)*(1.0-dx) + D(x2,y2,z)*dx) * dy) * (1.0-dz) +
	((D(x,y,z2)*(1.0-dx) + D(x2,y,z2)*dx) * (1.0-dy) +
	(D(x,y2,z2)*(1.0-dx) + D(x2,y2,z2)*dx) * dy) * dz);

and then return the reduced value from interp_val? I.e. get_val would then be returning a process-specific value when run with MPI?

[stevengj]

Yes. Do the reduce on the interpolated value.

This gives the same result because the interpolation is a linear operation.

[thchr]

Okay, thanks! (my uncertainty was just whether MPI allows a function to return a process-specific value or if it has to be reduced before returning; makes sense that it doesn't need that)

[stevengj]

Functions can do whatever they want. MPI is just a library providing subroutines for a bunch of running processes to communicate with one another — each process is an independent program that can do independent computations.

[stevengj]

The implementation looks correct to me. Would be nice to test a case with mpbi-mpi just to be sure (especially in 2d vs 3d), but on a quick skim I think you have that right. Nice job!

[thchr]

I implemented the suggestion about moving the mpi_allreduce_1 call to interp_val (thanks for explaining the mental mode of MPI to me: now I have a clearer picture of it).

I also compiled and checked the mpbi and mpbi-mpi versions in both 2D and 3D cases (for resolution = 64 and with 4 processes in MPI mode), see below:

2D example: (cylinders on a square lattice; .ctl file below):

Value type	mpb	mpb-mpi	mpbi	mpbi-mpi
H-norm	1.3853	1.3855	1.3856	1.3859
H-energy	1.9375	1.9380	1.9385	1.9391
D-norm	8.7010	8.7002	8.6987	8.6973
D-energy	5.8359	5.8348	5.8328	5.8310

3D example: (spheres on in a cubic lattice; .ctl file from original post)

Value type	mpb	mpb-mpi	mpbi	mpbi-mpi
H-norm	0.8381	0.8384	0.8360	0.8469
H-energy	0.7158	0.7164	0.7123	0.7308
D-norm	6.4777	6.4773	6.4799	6.4682
D-energy	6.8146	6.8119	6.8176	6.7925

The agreement seems to be pretty good across the board, so I reckon the combination of HAS_MPI and no SCALAR_COMPLEX seems to be working as well.

2D-example .ctl file

(set! resolution 64)
(set! num-bands 3)
(set! output-epsilon (lambda () (print "skipping output-epsilon\n")))

; define a simple geometry
(set! geometry-lattice (make lattice (size 1 1 no-size) 
                                     (basis1 1 0) (basis2 0 1)))
(set! geometry (list (make cylinder
                        (center 0 0 0) (radius 0.25) (height infinity)
                        (material (make dielectric (epsilon 13) )))))
; set a k-point
(set! k-points (list (vector3 0.3532 0.1512 0)))

(run) ; run the calculation

; compute some interpolating values at an arbitrary real-space point
(define r (vector3 0.13 0.12 0))
(define band-idx 3)

(print "\n\n")
(get-hfield band-idx)
(print "|H|:      " (vector3-norm (get-field-point r)) "\n")
(compute-field-energy)
(print "H-energy: " (get-energy-point r) "\n")

(get-dfield band-idx)
(print "|D|:      " (vector3-norm (get-field-point r)) "\n")
(compute-field-energy)
(print "D-energy: " (get-energy-point r) "\n")

[thchr]

Seems CI failed on building libctl: not sure why? Should be unrelated to the latest changes and this PR?

[thchr]

CI is good now after the libctl fix.

[stevengj]

Thanks!