Preparations for SCUFF-EM benchmarks

Former-commit-id: 9ab0dfb39833e0da0db78680a90984160e24e3e9

examples/scuff-em-benchmarks/materials/matprop.dat

@@ -0,0 +1,54 @@
+# Drude-Lorentz models, with the same constants as in qpms.
+# Compared to qpms, sign of the imaginary part is swapped here due
+# to different time-frequency transform convention.
+
+MATERIAL LDSilver
+# Electronvolt over reduced Planck's constant; the omegas and gammas
+# are defined in these units.
+   eh = 1519267460583196.5;
+   omegap = 9.01
+   f0 = 0.84
+   f1 = 0.065
+   f2 = 0.124
+   f3 = 0.111
+   f4 = 0.84
+   f5 = 5.646
+   omega0 = 0.0
+   omega1 = 0.816
+   omega2 = 4.481
+   omega3 = 8.185
+   omega4 = 9.083
+   omega5 = 20.29
+   gamma0 = 0.053
+   gamma1 = 3.886
+   gamma2 = 0.452
+   gamma3 = 0.065
+   gamma4 = 0.916
+   gamma5 = 2.419
+   Eps(w) = f0 * (omegap * eh)^2 / ((omega0 * eh)^2 - w^2 + I * w * gamma0) + f1 * (omegap * eh)^2 / ((omega1 * eh)^2 - w^2 + I * w * gamma1) + f2 * (omegap * eh)^2 / ((omega2 * eh)^2 - w^2 + I * w * gamma2) + f3 * (omegap * eh)^2 / ((omega3 * eh)^2 - w^2 + I * w * gamma3) + f4 * (omegap * eh)^2 / ((omega4 * eh)^2 - w^2 + I * w * gamma4) + f5 * (omegap * eh)^2 / ((omega5 * eh)^2 - w^2 + I * w * gamma5);
+ENDMATERIAL
+
+MATERIAL LDGold
+   eh = 1519267460583196.5;
+   omegap = 9.03
+   f0 = 0.76
+   f1 = 0.024
+   f2 = 0.010
+   f3 = 0.071
+   f4 = 0.601
+   f5 = 4.384
+   omega0 = 0.0
+   omega1 = 0.415
+   omega2 = 0.83
+   omega3 = 2.969
+   omega4 = 4.304
+   omega5 = 13.32
+   gamma0 = 0.053
+   gamma1 = 0.241
+   gamma2 = 0.345
+   gamma3 = 0.87
+   gamma4 = 2.494
+   gamma5 = 2.214
+   Eps(w) = f0 * (omegap * eh)^2 / ((omega0 * eh)^2 - w^2 + I * w * gamma0) + f1 * (omegap * eh)^2 / ((omega1 * eh)^2 - w^2 + I * w * gamma1) + f2 * (omegap * eh)^2 / ((omega2 * eh)^2 - w^2 + I * w * gamma2) + f3 * (omegap * eh)^2 / ((omega3 * eh)^2 - w^2 + I * w * gamma3) + f4 * (omegap * eh)^2 / ((omega4 * eh)^2 - w^2 + I * w * gamma4) + f5 * (omegap * eh)^2 / ((omega5 * eh)^2 - w^2 + I * w * gamma5);
+ENDMATERIAL
+

examples/scuff-em-benchmarks/scatter/20x20/Makefile

@@ -0,0 +1,4 @@
+
+
+freqlist_scuff: freqlist_eV
+	../../../../misc/omega_eV2scuff.py -o freqlist_scuff freqlist_eV 

examples/scuff-em-benchmarks/scatter/20x20/freqlist_eV

@@ -0,0 +1,3 @@
+1.0
+1.5
+2.0

examples/scuff-em-benchmarks/shapes/cylinder.geo.template

@@ -0,0 +1,62 @@
+// radius, height in nanometers
+R=RADIUS/1000.;
+h=HEIGHT/1000.;
+
+
+l=ELEMSIZ/1000.;
+
+// Circle centers
+Point(1) = { 0,  0,  h/2.,l};
+Point(2) = { 0,  0, -h/2.,l};
+
+// Upper circle arc limits
+Point(3) = { R,  0,  h/2.,l};
+Point(4) = { 0,  R,  h/2.,l};
+Point(5) = {-R,  0,  h/2.,l};
+Point(6) = { 0, -R,  h/2.,l};
+
+// Upper circle
+Circle(1) = {3, 1, 4};
+Circle(2) = {4, 1, 5};
+Circle(3) = {5, 1, 6};
+Circle(4) = {6, 1, 3};
+
+// Lower circle arc limits
+Point(7)  = { R,  0, -h/2.,l};
+Point(8)  = { 0,  R, -h/2.,l};
+Point(9)  = {-R,  0, -h/2.,l};
+Point(10) = { 0, -R, -h/2.,l};
+
+// Lower circle
+Circle(5) = { 7, 2,  8};
+Circle(6) = { 8, 2,  9};
+Circle(7) = { 9, 2, 10};
+Circle(8) = {10, 2,  7};
+
+Line( 9) = {3,7};
+Line(10) = {4,8};
+Line(11) = {5,9};
+Line(12) = {6,10};
+
+Line Loop(13) = {1, 2, 3, 4};
+Line Loop(14) = {5, 6, 7, 8};
+
+Line Loop(15) = {1, 10, -5,  -9};
+Line Loop(16) = {2, 11, -6, -10};
+Line Loop(17) = {3, 12, -7, -11};
+Line Loop(18) = {4,  9, -8, -12};
+
+Ruled Surface(19) = {13};
+Ruled Surface(20) = {14};
+Ruled Surface(21) = {15};
+Ruled Surface(22) = {16};
+Ruled Surface(23) = {17};
+Ruled Surface(24) = {18};
+
+
+
+Physical Surface(1) = {19,20,21,22,23,24};
+
+
+
+

examples/scuff-em-benchmarks/shapes/extrude.sh

@@ -0,0 +1,28 @@
+#!/bin/sh
+
+sed -e "s/RADIUS/30./g" -e "s/HEIGHT/30./g" -e "s/ELEMSIZ/4./g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r30_h30_fine.msh -2 
+
+sed -e "s/RADIUS/30./g" -e "s/HEIGHT/30./g" -e "s/ELEMSIZ/7./g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r30_h30.msh -2 
+
+sed -e "s/RADIUS/30./g" -e "s/HEIGHT/30./g" -e "s/ELEMSIZ/15./g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r30_h30_rough.msh -2 
+
+sed -e "s/RADIUS/30./g" -e "s/HEIGHT/30./g" -e "s/ELEMSIZ/25./g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r30_h30_veryrough.msh -2 
+
+sed -e "s/RADIUS/100./g" -e "s/HEIGHT/50./g" -e "s/ELEMSIZ/13.3/g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r100_h50_fine.msh -2 
+
+sed -e "s/RADIUS/100./g" -e "s/HEIGHT/50./g" -e "s/ELEMSIZ/23.3/g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r100_h50.msh -2 
+
+sed -e "s/RADIUS/100./g" -e "s/HEIGHT/50./g" -e "s/ELEMSIZ/50./g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r100_h50_rough.msh -2 
+
+sed -e "s/RADIUS/100./g" -e "s/HEIGHT/50./g" -e "s/ELEMSIZ/83.3/g" cylinder.geo.template >tmp.geo \
+        && gmsh tmp.geo -o cylinder_r100_h50_veryrough.msh -2 
+
+rm tmp.geo
+

examples/scuff-em-benchmarks/shapes/sphere.geo.template

@@ -0,0 +1,60 @@
+//
+// gmsh geometry specification for a sphere of radius R = RADIUS nm
+// 
+
+//************************************************************
+//* input parameters      
+//************************************************************
+R = RADIUS/1000.;    // radius
+
+//************************************************************
+//* meshing finenesses ***************************************
+//************************************************************
+l3 = ELEMSIZ/1000.;  // fineness at north pole
+l2 = ELEMSIZ/1000.;  // fineness at equator
+l1 = ELEMSIZ/1000.;  // fineness at south pole
+
+//************************************************************
+//* upper sphere *********************************************
+//************************************************************
+Point(1) = {  0 ,    0,  0.0,  l2};
+Point(2) = {  R,    0,  0.0,  l2};
+Point(3) = {  0 ,   R,  0.0,  l2};
+Circle(1) = {2,1,3};
+Point(4) = { -R,    0,  0.0,  l2};
+Point(5) = {   0,  -R,  0.0,  l2};
+Circle(2) = {3,1,4};
+Circle(3) = {4,1,5};
+Circle(4) = {5,1,2};
+Point(6) = {   0,    0,  0.0+R, l3};
+Point(7) = {   0,    0,  0.0-R, l1};
+Circle(5) = {3,1,6};
+Circle(6) = {6,1,5};
+Circle(7) = {5,1,7};
+Circle(8) = {7,1,3};
+Circle(9) = {2,1,7};
+Circle(10) = {7,1,4};
+Circle(11) = {4,1,6};
+Circle(12) = {6,1,2};
+Line Loop(13) = {2,8,-10};
+Ruled Surface(14) = {13};
+Line Loop(15) = {10,3,7};
+Ruled Surface(16) = {15};
+Line Loop(17) = {-8,-9,1};
+Ruled Surface(18) = {17};
+Line Loop(19) = {-11,-2,5};
+Ruled Surface(20) = {19};
+Line Loop(21) = {-5,-12,-1};
+Ruled Surface(22) = {21};
+Line Loop(23) = {-3,11,6};
+Ruled Surface(24) = {23};
+Line Loop(25) = {-7,4,9};
+Ruled Surface(26) = {25};
+Line Loop(27) = {-4,12,-6};
+Ruled Surface(28) = {27};
+Physical Surface(1) = {28,26,16,14,20,24,22,18};
+
+//************************************************************
+//* reference point to get outward-pointing surface normals right
+//************************************************************
+Physical Point(1) = {1};

misc/omega_eV2scuff.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+import sys
+import argparse
+
+ap = argparse.ArgumentParser()
+ap.add_argument("-o", type=str, help='Output file (if not specified, standard output).')
+ap.add_argument("input_file", type=str, help="Input file (if not specified, standard input).")
+
+a = ap.parse_args()
+
+if a.o:
+    output = open(a.o, 'w')
+else:
+    output = sys.stdout
+if a.input_file:
+    input = open(a.input_file, 'r')
+else:
+    input = sys.stdin
+
+from qpms.constants import eV2SU
+
+for l in input:
+    print(eV2SU(float(l)), file=output)
+
+output.close()
+input.close()
+