3D half-space simulation in C++ scuff-em returns NAN output

[amrit-poudel]

I am trying to compute (electric) Greens function (say Gzz with r= r') for the half-space geometry using scuff-em C++ library. However, I keep getting NaN output.

For your convenience I have attached my C++, file, .scuffgeo file, and .geo file

I am not sure what I am doing wrong. I understand that there is already an application that computes LDOS and/or greens function, but I would really like to understand why my approach does not seem to work.

I will appreciate any help/suggestion.

Thank you.

***** Halfspace.geo ***

// Mesh.ElementOrder = 2 ;

// Mesh.SecondOrderLinear = 1 ;

Mesh.Algorithm3D = 4 ; // 3D mesh algorithm (1=Delaunay(+Tetgen), 4=Frontal)

// Mesh.Optimize = 1 ; // Optimize the mesh to improve the quality of tetrahedral elements ( for Tetgen)

Mesh.OptimizeNetgen = 1; // Optimize the mesh using Netgen to improve the quality of tetrahedral elements

// Half space parameters

xlength = 1.0; ylength = 1.0;

x1 = 0.0; y1 = 0.0;

x2 = x1 + xlength; y2 = y1;

x3 = x2; y3 = y2 + ylength;

x4 = x1; y4 = y3;

h = 20.0; // height // This is the thickness of the half-space. It should be infinite, but I am using a finite but a large value.

// Mesh parameter lc = 0.4;

// Geometry Point(1) = {x1, y1, -h, lc}; Point(2) = {x2, y2, -h, lc}; Point(3) = {y3, y3, -h, lc}; Point(4) = {x4, y4, -h, lc};

Line(1) = {1, 2}; Line(2) = {2, 3}; Line(3) = {3, 4}; Line(4) = {4, 1};

Line Loop(5) = {1, 2, 3, 4};

// Surface Plane Surface(6) = {5};

// Volume vol[] = Extrude{0.0, 0.0, h}{Surface{6};};

/* vol contains in the following order: [0] - front surface (opposite to source surface) [1] - extruded volume [2], [3] , ... and so on */

// Save surfaces Physical Surface(1001) = {6}; Physical Surface(1002) = {vol[0]}; Physical Surface(1003) = {vol[2]}; Physical Surface(1004) = {vol[3]}; Physical Surface(1005) = {vol[4]}; Physical Surface(1006) = {vol[5]};

// Save volumes // Physical Volume(2001) = {vol[1]}; // We do not need volume discretization for BEM calculation

// Color mesh based on physical entity Mesh.ColorCarousel = 2; // (0=by element type, 1=by elementary entity, 2=by physical entity, 3=by partition)

/ If physical line surface or volume is not specificed then gmsh will save all 1D, 2D and 3D mesh elements. /

**** Halfspace.scuffgeo

Material specification

MATERIAL SILVER cond = 6.3e7 ; epsilon_rr = 3.9 ; epsilon_0 = 8.85e-12 ; Eps(w) = epsilon_rr - i_cond/(w_epsilon_0) ; ENDMATERIAL

Lattice

LATTICE VECTOR 1.0 0 VECTOR 0.0 1.0 ENDLATTICE

Object specification

OBJECT TheHalfspace MESHFILE Halfspace.msh MATERIAL SILVER ENDOBJECT

[amrit-poudel]

****** Script_halfspace_bem.cpp ***

// include all libraries here

define pi 3.14159265358979323

define MASTER 0

// Global variables

// ***** // SI units // ****

// Constants extern const double epsilon_0 = 8.85e-12; extern const double c = 3e8;

// Source parameters extern const double lambda_SI = 1e-6; extern const double freq_SI = c/lambda_SI; extern const double omega_SI = 2.0_pi_freq_SI; extern const double a = 1e-6;

// Material parameters extern const double freq_plasma_SI = 9.01_2.42e14; // 1eV = 2.417e14 Hz (linear frequency), \omegap = sqrt{\sigma/(\tau\epsilon_0)} extern const double sigma_SI = 6.3e7; extern const double epsilon_rr = 3.9;

// **** // Scuffem units // ****

// Source parameters extern const double freq = freq_SI_a/c; extern const double omega = omega_SI_a/c;

// Source and observation locations extern const double sx = 0.0; extern const double sy = 0.0; extern const double sz =0.5; extern const double szstart = sz; extern const double szend = 1.5; extern const double szstep = 0.5; extern const int sznum = floor((szend - szstart)/szstep);

// Main int main(int argc, char **argv) {

// Plot number
int plnum = atoi(argv[1]); // cannot use stoi because argv contains array of pointers to char (and not string type of C++)

char pstr[20]; sprintf(pstr, "%d", (int)plnum);

// Define geometry
RWGGeometry *Geom=new RWGGeometry("Halfspace.scuffgeo", 1.0e-6); //NOTE: I modified the source code such that RWGGeometry object  can take an extra input parameter which specifies the  characteristic length scale of the system. This parameter was hard coded to 1 micron (1e-6 m) in the original scuff-em.

// Preallocate BEM matrix and RHS vector
HMatrix *M  = Geom->AllocateBEMMatrix();
HVector *KN = Geom->AllocateRHSVector();

// Define BZ vectors
// NOTE: integration of any function f(k_x, k_y) over BZ is given by:
// 1/V_BZ\int f(k_x, k_y) dk_x dk_y, where 0<k_x<2*pi/L_x, 0<k_y<2*pi/L_y, V_BZ = 4*pi^2/(L_x*L_y)
double Lx = 1;
double Ly = 1;
double kxstart = 0.0; double kxend = 2.0*pi/Lx; cout << kxend << endl;
double kystart = 0.0; double kyend = 2.0*pi/Ly; cout << kyend << endl;
double kxstep = 0.1*2*pi/Lx; cout << kxstep << endl;
double kystep = 0.5*2*pi/Ly; cout << kystep << endl;
int kxnum, kynum;
kxnum = floor((kxend - kxstart)/kxstep);
kynum = floor((kyend - kystart)/kystep);

cout << "kxnum:" << kxnum << endl;
cout << "kynum:" << kynum << endl;

double kBloch[2] = {kxstart, kystart};

// Source
double Loc[3] = {sx, sy, sz};  // point source location 
cdouble Pol[3] = {0.0, 0.0, 1.0};  // point source polarization and strength
PointSource PS(Loc, Pol, 0); // use 1 for magnetic dipole type

// Store output field in a matrix
HMatrix FxMatrix(sznum, kxnum*kynum, LHM_COMPLEX, LHM_NORMAL); //  complex, non-symmetric matrix
HMatrix FyMatrix(sznum, kxnum*kynum, LHM_COMPLEX, LHM_NORMAL); //  complex, non-symmetric matrix
HMatrix FzMatrix(sznum, kxnum*kynum, LHM_COMPLEX, LHM_NORMAL); //  complex, non-symmetric matrix

cdouble EH[6];

int i, m, n;
for (m=0; m<kxnum; m++){

    for (n=0; n<kynum; n++) {

        // Assemble BEM matrix
        Geom->AssembleBEMMatrix(static_cast<cdouble>(omega), kBloch, M);

        // Factorize BEM matrix
        M->LUFactorize();

        for (i=0; i<sznum; i++) {

            // cout << " loop index:" << i << endl;

            // Assemble RHS vector
            Geom->AssembleRHSVector(static_cast<cdouble>(omega), kBloch, &PS, KN);
            // Geom->AssembleRHSVector(omega, &PW, KN);

            // Solve BEM equation
            M->LUSolve(KN);

            // Obtain total field at multiple points
            Geom->GetFields(0, KN, static_cast<cdouble>(omega), kBloch, Loc, EH);
            // Using &PS gives total field (scattered + incident). Replace &PS by 0 to get scattered field only

            // Store fields in HMatrix
            FxMatrix.SetEntry(i, m*kynum+n, EH[0]);
            FyMatrix.SetEntry(i, m*kynum+n, EH[1]);
            FzMatrix.SetEntry(i, m*kynum+n, EH[2]);

            // Update source location
            Loc[2] = Loc[2] + szstep;

            // Update source
            PS.SetX0(Loc);
        }
        // Restart Loc
        Loc[2] = szstart;

        // Update ky 
        kBloch[1] = kBloch[1] + kystep;
    }
    // Restart ky
    kBloch[1] = kystart;

    // Update kx 
    kBloch[0] = kBloch[0] + kxstep;
}
//Restart kx
kBloch[0] = kxstart;

// Save 
// Create file name 
char fname_Ex[255] = "./plots/Ex_" ; strcat(fname_Ex, pstr);  strcat(fname_Ex, ".hdf5");
char fname_Ey[255] = "./plots/Ey_" ; strcat(fname_Ey, pstr);  strcat(fname_Ey, ".hdf5");
char fname_Ez[255] = "./plots/Ez_" ; strcat(fname_Ez, pstr);  strcat(fname_Ez, ".hdf5");

// Open hdf5 file
void *pFile_Ex = FxMatrix.OpenHDF5Context(fname_Ex);
void *pFile_Ey = FyMatrix.OpenHDF5Context(fname_Ey);
void *pFile_Ez = FzMatrix.OpenHDF5Context(fname_Ez);

// Write to hdf5 file
FxMatrix.ExportToHDF5(pFile_Ex, "Ex");
FyMatrix.ExportToHDF5(pFile_Ey, "Ey");
FzMatrix.ExportToHDF5(pFile_Ez, "Ez");

// Close hdf5 file
FxMatrix.CloseHDF5Context(pFile_Ex);
FyMatrix.CloseHDF5Context(pFile_Ey);
FzMatrix.CloseHDF5Context(pFile_Ez);

return 0;

}

[HomerReid]

I looked through this and don't see anything immediately jumping out at me that would cause a NAN. Where is the NAN appearing?

A quick way to find NANs is to set the environment variable SCUFF_ABORT_ON_FPE=1. (FPE stands for "floating-point exception.") Then your code will abort as soon as any NAN is found anywhere. You can trace through in e.g. the debugger to find where this happened. It's easier to identify NANs at their source rather than trying to trace their propagation through the code.

One problem with your example, which is probably not related to the cause of the NAN, is that half-space geometries cannot be described using the simple OBJECT...ENDOBJECT syntax in the .scuffgeo file. The reason is that these geometries do not involve compact homogeneous regions bounded by closed surfaces. Instead, you have to use REGION and SURFACE statements. See the .scuffgeo files in this example:

http://homerreid.dyndns.org/scuff-em/doc/examples/HalfSpaceLDOS/HalfSpaceLDOS/

[amrit-poudel]

Thank you very much for your quick reply. I ran the simulation all the way to end and when I was trying to read hdf5 files, all the values were NaN.

I will try to put the environment variable in my C++ file, like you suggested, and run it again.

As for the .scuffgeo file for half-space geometry, I was simply following the following example from your webpage. I thought one could approximate half-space geometry by a periodic structure of a rectangular cuboid whose "surfaces straddle the unit-cell boundary" (using your own words).

" Surfaces that straddle the unit-cell boundary

More generally, your unit-cell geometry may include surfaces that straddle the unit-cell boundary. (This will be the case whenever the infinite surfaces you are describing are connected, as opposed to the isolated arrays of objects considered above). For example, here's .scuffgeo file representing a thin metallic film perforated by an array of holes. The parameters here are chosen to mimic the geometry studied by Martin-Moreno et al, Physical Review Letters 86 1114 (2001). "

http://homerreid.dyndns.org/scuff-EM/reference/scuffEMGeometries.shtml

Could you please clarify why this approximation is not accurate ? May be I am missing something here.

That said, I understand that repeated rectangular cuboids do not quite form the half-space geometry. Instead, they form an infinitely extended metallic film (in x and y directions) with finite thickness in the z-direction. However, one can show analytically, that the greens function for thin film approaches that of the half-space, if the source/evaluation point << film thickness.

This is exactly what I am trying to test with your scuff-em c++ library.

If this is incorrect, I will follow your suggestion and use REGION and SURFACE.

Thanks again for answering my questions. I really appreciate it.

[HomerReid]

The Martin-Moreno geometry is topologically distinct from a half-space geometry. The regions above and below the perforated film belong to the same contiguous region of space. Given points above and below the perforated film, I can find a path connecting the points that lies entirely in the exterior region (in this case the path would go through one of the holes in the film). From the boundary-element perspective, the topology of the exterior medium in this case is equivalent to the simplest geometries involving compact objects in vacuum. The fact that the perforated thin film is not compact, but is instead extended, does not change this topology (although it does require care when meshing the unit cell, as discussed on the web page you referenced above).

In contrast, for a half-space geometry the interface surface splits the exterior medium into two infinite disconnected regions. There is no way to describe such a configuration as one or more connected homogeneous regions embedded in a connected exterior region.

The same would be true for a non-perforated thin film, like in this example:

http://homerreid.dyndns.org/scuff-EM/scuff-transmission/#ThinFilm

Again, there is no way to describe this geometry in terms of homogeneous regions embedded in a connected exterior region.

[amrit-poudel]

Thanks for your reply. This makes sense.

I modified my .scuffgeo file, following your suggestion.

Material specification

MATERIAL SILVER cond = 6.3e7 ; epsilon_rr = 3.9 ; epsilon_0 = 8.85e-12 ; Eps(w) = epsilon_rr - i_cond/(w_epsilon_0) ; ENDMATERIAL

Lattice

LATTICE VECTOR 1.0 0 VECTOR 0.0 1.0 ENDLATTICE

**

Half space geometry with top metallic surface at z=0

REGION Upperemptyspace MATERIAL Vacuum REGION Halfspace MATERIAL SILVER

SURFACE Halfspacesurface MESHFILE Halfspacesurface.msh REGIONS Upperemptyspace Halfspace ENDSURFACE

One last thing I wanted to confirm was units of frequency. You have already answered this question before, so i am sorry for repeating here, but just wanted to make sure.

In my .scuffgeo file above, I have units of material parameters expressed in SI units. As such, eps(\omega) expects \omega in SI units (without scaling by any factor).

However, in my C++ file (above), I pass dimensionless \omega, that is \omega = \omega_SI*a/c, where \omega_SI = angular frequency in Hz, c= 3e8 m/s, and a = characteristic length (default of 1 micron, 1e-6 m)

// Define geometry double a = 1.0e-6. RWGGeometry *Geom=new RWGGeometry("Halfspace.scuffgeo", a); // modified source code

// Assemble BEM matrix Geom->AssembleBEMMatrix(static_cast(omega), kBloch, M); // dimensionless omega

Will .scuffgeo file pick up the correct units of \omega?

[amrit-poudel]

My another question is related to your .geo file for half space geometry.

In your .geo file, you have used Extrude { 0, L, 0 } { Line{1}; Layers{N}; }

Is it really important to extrude both the line and the mesh together for periodic geometry calculation in scuffem?

Does the mesh must be structured for periodic geometry? What about for non-periodic geometry?

Will it be in accurate to define the "surface" for half-space geometry in the following way:

Line(1) = {1, 2}; Line(2) = {2, 3}; Line(3) = {3, 4}; Line(4) = {4, 1};

Line Loop(5) = {1, 2, 3, 4};

// Surface Plane Surface(6) = {5};

// Save surfaces Physical Surface(7) = {6}; // gmsh will save the mesh for surface 6 only!

Thank you for your time and help.

[HomerReid]

Not sure I understand your question. As discussed in the online documentation, the material-property description expects the w variable to be specified in units of rad/sec, as is customary for frequency-dependent dielectric-function data as tabulated in standard references, while the command line --omega argument is interpreted in units of 3e14 rad/sec, so that typical values of --omega for nanophotonics problems are on the order of 1.0 as opposed to 1.0e14. My assumption is that the 14 orders of magnitude separating these two unit systems should make it easy to figure out when there is a mismatch.

Here is the dielectric model I use for silver:

# this model taken from:
# http://www.opticsinfobase.org/abstract.cfm?URI=ao-24-24-4493
MATERIAL Silver
  wp    = 1.37e16
  gamma = 2.73e13
  Eps(w) = 1 + wp^2 / (i*w*(i*w+gamma));
ENDMATERIAL

Regarding the question about meshing, I have found the "extrusion" feature in GMSH to be useful for ensuring that unit-cell meshes satisfy the criterion required for SCUFF-EM surface meshes that straddle unit-cell boundaries. (As discussed in the online documentation, this criterion is that each triangle edge that lies on one segment of a unit-cell boundary must have an image edge on the opposite segment of the unit-cell boundary.) But this is certainly not the only way to ensure this constraint, and indeed you could always design and describe your own surface mesh that satisfies this constraint, bypassing GMSH entirely if you like.

[amrit-poudel]

Sorry for not being very clear, but I understood what you said. I just need to pass dimensionless omega in C++ file, and this gets converted to appropriate units internally.

As for the mesh, I assume this image edge constraint is limited to periodic geometry that straddle unit cell boundaries.

For all other periodic or non periodic geometries, I assume one can use structured or unstructured meshes without any constraints.

Thank you for all your help.

[amrit-poudel]

Could you please clarify why you do not use the keyword "Physical" in your .geo file? Not specifying "Physical" in .geo file will result in a mesh file with all kinds of elements saved.

I used your .geo file (Square_N.geo from http://homerreid.dyndns.org/scuff-em/doc/examples/HalfSpaceLDOS/HalfSpaceLDOS/) to generate a 2D surface for half-space geometry calculation, but the mesh file that GMSH generated when I clicked Mesh->Define->2D (in gmsh GUI), there well all kinds of elements saved, with some element type with only two nodes (these are the elements on a line where you can have only two nodes per element. However, we do not need these elements for BEM calculation).

However, if you specify Physical keyword, then GMSH will only save elements of that surface of interest. This surface will then have triangular elements with 3 nodes.

[HomerReid]

The .msh file reader in SCUFF-EM ignores elements that do not have exactly 3 nodes.

[amrit-poudel]

Thanks for clarifying.

I have been trying to run a 3D half-space geometry calculation with

export SCUFF_ABORT_ON_FPE=1, to determine where the program comes across NaN.

However, the calculation seems to get stuck (but does not abort) at matrix assembly:

Geom->AssembleBEMMatrix(omega, kBloch, M);

---

Here is my scuffgeo file.

\ Material specification MATERIAL SILVER cond = 6.3e7 ; epsilon_rr = 3.9 ; epsilon_0 = 8.85e-12 ; Eps(w) = epsilon_rr - i_cond/(w_epsilon_0) ; ENDMATERIAL

\ Lattice LATTICE VECTOR 1.0 0 VECTOR 0.0 1.0 ENDLATTICE

\ ** \ Half space geometry with top metallic surface at z=0

\ This is default: \ REGION Exterior MATERIAL VACUUM

REGION Upperemptyspace MATERIAL VACUUM REGION Lowermetalspace MATERIAL SILVER

SURFACE Halfspacesurface MESHFILE Halfspacesurface.msh REGIONS Upperemptyspace Lowermetalspace ENDSURFACE

I used your .geo file to specify half-space surface geometry, except that I use N =10.

My program has been running for past 38 hrs using 4 cores (intel ivybridge), but has not finished assembling a single BEM matrix (for a single bloch vector) yet.

Is it suppose to take this long (which looks unreasonable to me) ? Thanks.

Update: Finally, the program aborted due to FPE error during BEM assembly.

I am not sure why this happens. I am using your example .geo file for 3D half space calculation.

[amrit-poudel]

Do you have any insight why the program returns NaN when it tries to assemble BEM matrix?

I am using your .geo file to generate the mesh file and .scuffgeo file is attached above.

The program runs fine, but returns NaN.

I did not compile the library with -g option, so I am unable to use any debugger.

I will appreciate any other suggestions.

[HomerReid]

Looking at your source code, I don't see anything obviously wrong, but I also did not 100% follow the series of numerical constant definitions at the beginning. It's possible that somehow there is a funny number getting passed to AssembleBEMMatrix---for example, if the lengthscale or frequency scale is getting configured incorrectly, then the code might attempt to evaluate the dielectric function of silver at a frequency of something like 3.0 rad/sec or even 3.0e^{-14} rad/sec instead of something on the order of 3e14 rad/sec, in which case the code might be erroneously trying to do calculations at a wavelength that is something like 1000x (or even 1e14 x) smaller than the triangle size. This would explain why the BEM matrix assembly bogged down before it crashed---when the wavelength is very small the code switches over to a more expensive technique for computing BEM matrix elements.

Really the only way to debug these things is to run the code in the debugger, break on entry to AssembleBEMMatrix, and print out the numerical values of various parameters such as Omega and later Eps and Mu.

Also, try running your code with a much coarser mesh (try N=3 in the .geo file) to see if the problem is still there. If so, you can run in the debugger, and the FPE abort will happen much sooner, so you can report which line is causing the crash.

[amrit-poudel]

I reduced the mesh size to N =3, and ran the program in debug mode. It takes more than 30 hrs to get to the point where it aborts. I do not have sufficient computing resources right now to run this on debug mode for that long (I printed out value of \omega, and that looks fine), so I could not identify which part of AssembleBEMMatrix is causing the problem.

I would be interested in running your C++ script for half-space geometry on my cluster, if that is okay.

log file looks like this. I do not understand why there are so many grid points (with N=3 surface of the half-space)!

08/31/15::15:22:43: Assembling BEM matrix block (0,0) 08/31/15::15:22:43: Step 1: Contributions of innermost grid cells... 08/31/15::15:22:43: ...(1,1) block... 08/31/15::15:22:44: ...(1,0) block... 08/31/15::15:22:47: ...(1,-1) block... 08/31/15::15:22:48: ...(0,1) block... 08/31/15::15:22:50: ...(0,0) block... 08/31/15::15:22:53: Step 2: Contributions of outer grid cells... 08/31/15::15:22:53: Creating GBar accelerator for region Upperemptyspace, surfaces Halfspacesurface,Halfspacesurface 08/31/15::15:22:54: Creating interpolation grid with (95x95x2)=18050 grid points... 08/31/15::15:23:00: ...interpolation table constructed! 08/31/15::15:23:00: Creating GBar accelerator for region Lowermetalspace, surfaces Halfspacesurface,Halfspacesurface 08/31/15::15:23:00: Cutting off interpolation table at Lx=LMax=5.128483e-02. 08/31/15::15:23:00: Cutting off interpolation table at Ly=LMax=5.128483e-02. 08/31/15::15:23:03: Creating interpolation grid with (1839x1839x2)=6763842 grid points...

It gets stuck here.

[HomerReid]

This can happen in cases where the wavelength in the material is very short (in metals, for example) so the code thinks it has to sample the periodic Green's function with extremely high precision.

This behavior can be controlled by setting the environment variable SCUFF_INTERPOLATION_TOLERANCE. Try saying

% export SCUFF_INTERPOLATION_TOLERANCE=1e-2

to see if that improves the situation.

[amrit-poudel]

I used that option and the program reaches to the point where it encounters FPE much faster than before.

From the debugger , it looks like AddGFullTerm function in GBarAccelerator.cc is the source of error. (the debugger fails to display the source code, so I am not sure which line in AddGFullTerm exactly, but it looks like exp is causing the error.)

I will appreciate any other suggestions.

[HomerReid]

OK, so, here's how to pinpoint what is going wrong:

• Run the code in GDB until it crashes.
• At the GDB prompt, type where. This will tell you the precise line of C++ code on which the code crashed, for example GBarAccelerator.cc:123 or something (where 123 is the line number).
• Set a breakpoint at that line (break GBarAccelerator.cc:123) and run the code again from the beginning.
• Now GDB will stop at the line that causes the core dump, and you can inspect the various variables (type things like print Omega, etc) in an attempt to deduce what is provoking the FPE.

If you can't see what's going wrong, post a screenshot or console dump of your GDB session and we'll take a look.

As I said before, I suspect the issue is related to the change you made regarding lengthscales and frequency units. I think the dielectric function of silver is getting evaluated at a frequency that is something like 3.0 instead of something like 3.0e14, and thus the value of Epsilon is getting set to something like -1e15 + 1e23i instead of something more reasonable like -1e2 + 1e2i and this is causing numerical problems.

BTW, what happens if you run with the same geometry and input files using my scuff-ldos code as discussed in the following example?

http://homerreid.dyndns.org/scuff-em/doc/examples/HalfSpaceLDOS/HalfSpaceLDOS/

[amrit-poudel]

Here is the output from where:

It appears to me that the error is from exp function of the math library, whose source code is not accessible to me since I did not compile that library. This is why gdb does not point to a specific line of the source code.

Does this have to do with the intel's mkl library?

(gdb) where

0 0x000000354060ebb0 in __ieee754_exp () from /lib64/libm.so.6

1 0x0000003540626558 in cexp () from /lib64/libm.so.6

2 0x00002aaaaaaee22b in scuff::AddGFullTerm(double, std::complex, double, double, double, std::complex_) ()

from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

3 0x00002aaaaaaf01ef in scuff::GetGBar2D(double, scuff::GBarAccelerator, std::complex, std::complex_) ()

from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

4 0x00002aaaaaaf051a in scuff::GetGBar(double, scuff::GBarAccelerator, std::complex, std::complex, bool) ()

from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

5 0x00002aaaaaad9241 in scuff::AssembleInnerPPIIntegrand(double, double, double, double, double, std::complex, scuff::GBarAccelerator, bool, int, int, double, std::complex, std::complex, std::complex_) ()

from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

6 0x00002aaaaaad9f9c in scuff::GetPPIsCubature(scuff::GetPPIArgStruct, int, int, double_, double, double**, double_) ()

from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

7 0x00002aaaaaada4f6 in scuff::GetPanelPanelInteractions(scuff::GetPPIArgStruct_) () from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

8 0x00002aaaaaadaf4b in scuff::GetPanelPanelInteractions(scuff::GetPPIArgStruct, std::complex, std::complex, std::complex) () from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

9 0x00002aaaaaad85ed in scuff::GetEdgeEdgeInteractions(scuff::GetEEIArgStruct_) () from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

10 0x00002aaaaaad6a95 in scuff::GSSIThread(void_) () from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

11 0x00002aaaaaad77dd in scuff::GetSurfaceSurfaceInteractions () from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

12 0x00002aaaaaad7d2f in scuff::GetSurfaceSurfaceInteractions(scuff::GetSSIArgStruct_) ()

from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

13 0x00002aaaaaad58a3 in scuff::RWGGeometry::AssembleBEMMatrixBlock(int, int, std::complex, double, HMatrix*, HMatrix**, int, int, void, bool, int, HMatrix_, double) () from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

14 0x00002aaaaaad5ee4 in scuff::RWGGeometry::AssembleBEMMatrix(std::complex, double, HMatrix) ()

from /home/api852/HPC_LIB/INTEL/DEBUG_SCUFFEM/lib/libscuff.so.0

15 0x0000000000401996 in main ()

[HomerReid]

In GDB type up and then up again. This will take you up to the line in AddGFullTerm that makes the call to cexp. Then use print to print out the value of the argument that is being passed to cexp. This will be some crazily huge number. Then figure out why this number is getting set to a crazily huge value. It is probably the product of several factors, and one of the factors is some crazily huge number, and the reason for this is almost certainly that your silver dielectric function is getting evaluated at a frequency like 1.0 instead of a frequency like 1.0e14.

[amrit-poudel]

I got the following value for omega

(gdb) print omega $1 = 6.2831853071795853

Looks like R is invalid/huge.

(gdb) print R $2 = {< invalid float value >, < invalid float value >, -inf, 0, < invalid float value >, < invalid float value >}

[HomerReid]

No, that number looks fine. The value of Omega should be on the order of 1. The difficulty would only be present inside the code that evaluates the material properties of Silver. That code is supposed to be multipying Omega by the appropriate conversion factor (in this case, 3e14) to obtain a number that can be plugged into your definition of the dielectric function of silver. My hypothesis is that the conversion factor inside that code is getting set incorrectly. The symptom of this would be values of Epsilon that are much too large.

At any rate, a value of Omega like the one you are seeing would not cause the cexp routine to have floating-point errors. What is the argument of cexp when it crashes?

[amrit-poudel]

I am unable to use a print statement in gdb on any other arguments (other than R) of exp [ AddGFullTerm() function does not directly call cexp]

// source code

void AddGFullTerm(double R[3], cdouble k, double kBloch[2], double Lx, double Ly, cdouble T[10]) { double RmL[3]; double r2, r; cdouble PhaseFactor, IKR, Phi, Psi, Zeta, Upsilon;

PhaseFactor=exp( II_(Lx_kBloch[0] + Ly*kBloch[1]) );

RmL[0]=R[0] - Lx; RmL[1]=R[1] - Ly; RmL[2]=R[2];

r2=RmL[0]_RmL[0] + RmL[1]_RmL[1] + RmL[2]_RmL[2]; r=sqrt(r2); if ( r < 1.0e-8 ) return; IKR=II_k_r; Phi=exp(IKR)/(4.0_M_PIr); T[0] += PhaseFactor \ Phi;

Psi=(IKR-1.0)Phi/r2; Zeta=(3.0 + IKR(-3.0 + IKR))_Phi/(r2r2); Upsilon=(-15.0 + IKR(15.0 + IKR_(-6.0 + IKR)))_Phi/(r2_r2*r2);

T[1] += PhaseFactor * RmL[0] * Psi; T[2] += PhaseFactor * RmL[1] * Psi; T[3] += PhaseFactor * RmL[2] * Psi; T[4] += PhaseFactor * (Psi + RmL[0] * RmL[0] * Zeta); T[5] += PhaseFactor * ( RmL[0] * RmL[1] * Zeta); T[6] += PhaseFactor * ( RmL[0] * RmL[2] * Zeta); T[7] += PhaseFactor * (Psi + RmL[1] * RmL[1] * Zeta); T[8] += PhaseFactor * ( RmL[1] * RmL[2] * Zeta); T[9] += PhaseFactor * (Psi + RmL[2] * RmL[2] * Zeta);

}

[HomerReid]

The problem is almost certainly the value of k, which elsewhere gets set to sqrt(Epsilon*Mu)*Omega. Please check the value of k. If it is a crazily large value, please check the value of Epsilon that the code is obtaining from calls to e.g. GetEps().

[amrit-poudel]

I printed out the value of k in

void AddGFullTerm(double R[3], cdouble k, double kBloch[2], double Lx, double Ly, cdouble T[10])

The real part of k = omega = 6.2831853071795853, and the imaginary part of k =0.

[HomerReid]

There are two calls to exp in this routine. Do you know which one of them is causing the FPE? If you are having trouble getting GDB to print the value of the argument that is being passed to this function, that is probably because GDB doesn't know how to multiply complex numbers. So, just print the values of each of the various factors individually. For example, if the problematic call to exp is the first one, whose argument is `II*(Lx*kBloch[0] + Ly*kBloch[1]), then print LX, print kBloch[0], etc.

[HomerReid]

What value of kBloch is your code passing to AssembleBEMMatrix? The components of this vector should be numbers on the order of 1.0.

[amrit-poudel]

From my code:

double Lx = 1; double Ly = 1; double kxstart = 0.0; double kxend = 2.0*pi/Lx;

double kystart = 0.0; double kyend = 2.0*pi/Ly;

double kBloch[2] = {kxstart, kystart};

kx and ky are all of the order of 1.0

---

The problem with gdb is that it does not recognize Lx, Ly or kBloch etc. I am just using Log(" ") to print out the values. (I could have used cout also)

(gdb) print kBloch[0] No symbol "kBloch" in current context.

---

UPDATE: I just found out that the value of k (when FPE occurs) is of the order of 100!

real k:2.732671e+02. and imag k: -2.729852e+02.

I do not understand why the value of k changes. should not it be constant? (of the order of sqrt{eps}*w/c)

Initially, the value of k was

real part of k = omega = 6.2831853071795853, and the imaginary part of k =0.

[HomerReid]

Ah, this explains it. The problem is not that k is too large, but that its imaginary part is negative. This means that the real part of the quantity i*k is positive and large, so your call to exp is attempting to evaluate a number like exp(+272) and that is causing a floating-point exception

This is happening because your dielectric function is designed for use with the electrical-engineering sign convention (in which the time dependence of harmonic fields and currents is exp( + j \omega t) ), whereas SCUFF-EM uses the physics convention of exp( - i \omega t). In the physics convention, the imaginary part of the dielectric function needs to be positive, whereas in the EE convention it needs to be negative.

The fix is easy: just change your dielectric function to read

Eps(w) = epsilon_rr + i*cond/(w*epsilon_0) ;

instead of what you have now, which is

Eps(w) = epsilon_rr - i*cond/(w*epsilon_0) ;

[amrit-poudel]

It looks like that was it! Thanks a lot!

I will close this thread after I reproduce your result for half-space geometry. Thank you very much for your time and help. I really appreciate it.

UPDATE:

I just found out that using

cdouble omega(2*pi, 0.0);

gives NaN for first few values of Bloch wave vector.

[HomerReid]

There seem to be several unrelated issues all jumbled together on this issue thread. If the NAN issue is resolved, could you please close that issue and file new issues for the other points? I will be happy to respond to them individually.

[amrit-poudel]

Sorry about that. I moved all my comments to a new thread and will now close this thread. Thanks!