File: stress_tensor.cpp

package info (click to toggle)
meep-mpich2 1.7.0-3
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 25,824 kB
  • sloc: cpp: 27,370; python: 10,574; lisp: 1,213; makefile: 440; sh: 28
file content (87 lines) | stat: -rw-r--r-- 2,739 bytes parent folder | download | duplicates (10)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#include <meep.hpp>
using namespace meep;
using namespace std;

const double sx = 5.0;
const double sy = 3.0;
const double dpml = 1.0;
const double d = 0.35;
const double sw = 1.0;
const double res = 20;

double two_waveguides(const vec &p) {
  if ((fabs(p.x()) >= 0.5*d) && (fabs(p.x()) <= 0.5*d+sw) && (p.y() <= 0.5*sw) && (p.y() >= -0.5*sw))
    return 11.9;
  else
    return 1.0;
}

int main(int argc, char **argv) {
  initialize mpi(argc, argv);
  quiet = true;
  grid_volume gv = vol3d(sx+2*dpml,sy+2*dpml,0,res);
  gv.center_origin();
  const symmetry S = mirror(X,gv) - mirror(Y,gv);
  structure s(gv, two_waveguides, pml(dpml,X)+pml(dpml,Y), S);
  s.set_epsilon(two_waveguides);
  fields f(&s);
  f.use_bloch(vec(0,0,0.5));
  f.add_point_source(Ey, 0.22, 0.06, 0.0, 4.0, vec(0.5*(d+sw),0,0));
  f.add_point_source(Ey, 0.22, 0.06, 0.0, 4.0, vec(-0.5*(d+sw),0,0));

#if 0
  double T = f.last_source_time();
  int iT = T / f.dt;
  while (f.t < iT) {
    if (f.t % (iT / 10) == 0)
      master_printf("%g%% done with source\n", f.time()/T * 100);
    f.step();
  }
  double T2 = 300;
  int iT2 = T2 / f.dt;
  complex<double> *vals = new complex<double>[iT2];
  while (f.t - iT < iT2) {
    if ((f.t - iT) % (iT2 / 10) == 0)
      master_printf("%g%% done with harminv\n", (f.t - iT) * 100.0 / iT2);
    f.step();
    vals[f.t - iT - 1] = f.get_field(Ey, vec(0.5*(d+sw),0.,0.));
  }
  complex<double> amps[8];
  double freqs_re[8], freqs_im[8];

  master_printf("done with timestepping, running harminv...\n");
  int num = do_harminv(vals, iT2, f.dt, 0.19, 0.25, 8, 
  		       amps, freqs_re, freqs_im);
  master_printf("harminv found %d modes\n",num);
  for (int i=0;i<num;i++)
    master_printf("freq %d is %0.6g, %0.6g\n",i,freqs_re[i],freqs_im[i]);
  double freq = freqs_re[0];
  double T3 = T + T2;
#else
  double freq = 0.220847; // hard-code since we are testing force, not harminv
  double T3 = f.last_source_time() + 300;
#endif

  double dpad = 0.1;
  volume box(vec(0.5*d-dpad,-0.5*sw-dpad,0),vec(0.5*d+sw+dpad,0.5*sw+dpad,0));
  dft_flux fluxR = f.add_dft_flux_plane(box, freq, freq, 1);
  dpad = 0.02;
  volume_list line(volume(vec(0.5*d-dpad,-0.5*sy,0),vec(0.5*d-dpad,0.5*sy,0)),
		   Sx);
  dft_force forceR = f.add_dft_force(&line, freq, freq, 1);
  f.zero_fields();
  f.t = 0;
  while (f.time() < T3) {
    if (f.t % 2000 == 0)
      master_printf("%g%% done with flux and force\n", f.time()/T3 * 100);
    f.step();
  }
  double *fl  = fluxR.flux();
  double *fr = forceR.force();
  master_printf("flux is %0.8g, force is %0.8g, F/P = %0.8g\n",
		fl[0], fr[0], -0.5*fr[0]/fl[0]);
  double FoverP = -0.5*fr[0]/fl[0];
  delete fl; delete fr;
  return fabs(FoverP+0.33628872)/0.33628872 > 0.1;
  // MPB: -0.33628872
}