/
absorber-1d-ll.cpp
128 lines (110 loc) · 4.38 KB
/
absorber-1d-ll.cpp
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
/***************************************************************/
/* test for absorber functionality in libmeepgeom */
/* modeled after absorber_1d.ctl */
/***************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <complex>
#include "meep.hpp"
#include "ctl-math.h"
#include "ctlgeom.h"
#include "meepgeom.hpp"
#ifndef DATADIR
#define DATADIR "./"
#endif
using namespace meep;
typedef std::complex<double> cdouble;
/***************************************************************/
/* dummy material function needed to pass to structure( ) */
/* constructor as a placeholder before we can call */
/* set_materials_from_geometry */
/***************************************************************/
double dummy_eps(const vec &) { return 1.0; }
/***************************************************************/
/* usage: absorber-1d [ --pml ] */
/***************************************************************/
int main(int argc, char *argv[]) {
initialize mpi(argc, argv);
// simple argument parsing
bool use_pml = false;
bool verbose = false;
for (int narg = 1; narg < argc; narg++) {
if (argv[narg] && !strcmp(argv[narg], "--pml"))
use_pml = true;
else if (argv[narg] && !strcmp(argv[narg], "--verbose"))
verbose = true;
else
meep::abort("unrecognized command-line option %s", argv[narg]);
};
if (verbose) master_printf("Using %s.\n", use_pml ? "pml" : "absorber");
double resolution = 20.0;
geometry_lattice.size.x = meep_geom::TINY;
geometry_lattice.size.y = meep_geom::TINY;
geometry_lattice.size.z = 10.0;
grid_volume gv = volone(10.0, resolution);
gv.center_origin();
boundary_region br = use_pml ? pml(1.0) : boundary_region();
structure the_structure(gv, dummy_eps, br);
meep_geom::absorber_list alist = 0;
if (!use_pml) {
alist = meep_geom::create_absorber_list();
meep_geom::add_absorbing_layer(alist, 1.0, Z);
};
geometric_object_list g = {0, 0};
vector3 center = {0, 0, 0};
meep_geom::set_materials_from_geometry(&the_structure, g, center,
true, // use_anisotropic_averaging,
DEFAULT_SUBPIXEL_TOL, // tol
DEFAULT_SUBPIXEL_MAXEVAL, // maxeval
false, // ensure_periodicity
meep_geom::vacuum, alist);
if (alist) meep_geom::destroy_absorber_list(alist);
// (set! sources (list (make source (src (make gaussian-src (frequency (/ 0.803)) (fwidth 0.1)))
// (center 0 0 0) (component Ex))))
fields f(&the_structure);
double fcen = 1.0 / 0.803;
double df = 0.1;
gaussian_src_time src(fcen, df);
vec x0(0.0);
f.add_point_source(Ex, src, x0);
double min_time = 50.0;
double dt_output = 10.0;
double next_output = dt_output;
double start_time = f.time();
double max_field = real(f.get_field(Ex, x0));
double field = max_field;
double f50 = 0.0;
bool done = false;
while (!done) {
f.step();
double t = f.time();
field = real(f.get_field(Ex, x0));
if (t >= 50.0 && f50 == 0.0) f50 = field;
max_field = fmax(fabs(field), max_field);
if (t >= next_output) {
next_output += dt_output;
if (verbose) master_printf("%.6e: %+.12e \n", t, field);
};
done = ((t - start_time) > min_time && (fabs(field) <= 1.0e-6 * fabs(max_field)));
};
double tFinal = f.time();
double fFinal = real(f.get_field(Ex, x0));
// test: compare {f50, tFinal, fFinal} to
// their reference values
double f50_ref = -5.090114e-01;
double tFinal_ref = 9.195000e+01;
double fFinal_ref = sizeof(realnum) == sizeof(float) ? 1.413336e-07 : 1.624782e-07;
if (fabs(f50 - f50_ref) > 1.0e-6 * fabs(f50_ref) ||
fabs(tFinal - tFinal_ref) > 1.0e-6 * fabs(tFinal_ref) ||
fabs(fFinal - fFinal_ref) > 1.0e-6 * fabs(fFinal_ref)) {
master_printf("{f50, tFinal, fFinal}={%e,%e,%e}\n", f50, tFinal, fFinal);
master_printf(" should be:\n");
master_printf("{f50, tFinal, fFinal}={%e,%e,%e}\n", f50_ref, tFinal_ref, fFinal_ref);
meep::abort("Test failed.");
}
else if (verbose)
master_printf("Test successful.\n");
meep_geom::unset_default_material();
return 0;
}