HomerReid / scuff-em

A comprehensive and full-featured computational physics suite for boundary-element analysis of electromagnetic scattering, fluctuation-induced phenomena (Casimir forces and radiative heat transfer), nanophotonics, RF device engineering, electrostatics, and more. Includes a core library with C++ and python APIs as well as many command-line applications.
http://www.homerreid.com/scuff-em
GNU General Public License v2.0
125 stars 50 forks source link

Different results from python GetPFT and scuff-scatter #195

Closed jmllorens closed 5 years ago

jmllorens commented 5 years ago

Hi Homer,

I found the new interface in libscuffSolver very useful and easy to handle. I am starting to migrate some of my models from command line to python, and I found a problem with the output of GetPFT which differs from that of scuff-scatter. The scattered power is the same, but the absorbed power and rest of magnitudes differ. I have been browsing the source code, but I haven't found any remarkable difference between the different calls.

In the attached zip file, there is a small example with a simple mesh file, the scuffgeo file, Args and a short python script.

Any help is welcome.

José.

test.zip

jmllorens commented 5 years ago

I think I have finally found the problems. The first one is that CachedIF in scuffSolver.cc always gets a copy of a PointSource. The second one is that the copy method of PlaneWave and PointSource does not copy the attributes of the base class IncField.

  1. In scuffSolver.cc:684 the copy of the source to CachedIF takes place:

    else if (IF)
    { G->AssembleRHSVector(Omega, IF, KN);
     if (CachedIF){
       delete CachedIF;
     }
     PlaneWave *PW = (PlaneWave *) IF;
     if (PW) CachedIF=new PlaneWave(*PW);
     PointSource *PS = (PointSource *) IF;
     if (PS) CachedIF=new PointSource(*PS);
    
    }

    The two if clauses are always true and CachedIF ends up with a copy of PS.

  2. The copy method of PlaneWave in PlaneWave.cc:50:

    PlaneWave::PlaneWave(const PlaneWave &PW)
    { 
    memcpy(E0, PW.E0, 3*sizeof(cdouble));
    memcpy(nHat, PW.nHat, 3*sizeof(double));
    SetRegionLabel(PW.RegionLabel);
    }

    does not include a copy of the IncField attributes. Therefore, Omega is equal to -1 +0j. The same problem happens for PointSource.

It would be nice to find a nice solution for both problems. I can only make a dirty patch at the moment.

José.

jmllorens commented 5 years ago

The pull-request #202 can fix this, although some refinements are welcome.