← Back to branch summary

~ubuntu-branches/debian/sid/meep-mpich2/sid

« back to all changes in this revision

Viewing changes to tests/two_dimensional.cpp

  • Committer: Package Import Robot
  • Author(s): Thorsten Alteholz
  • Date: 2012-06-06 18:00:00 UTC
  • Revision ID: package-import@ubuntu.com-20120606180000-037wajw9zkgziel5
Tags: upstream-1.1.1
Import upstream version 1.1.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/two_dimensional.cpp
 
1
/* Copyright (C) 2005-2009 Massachusetts Institute of Technology  
 
2
%
 
3
%  This program is free software; you can redistribute it and/or modify
 
4
%  it under the terms of the GNU General Public License as published by
 
5
%  the Free Software Foundation; either version 2, or (at your option)
 
6
%  any later version.
 
7
%
 
8
%  This program is distributed in the hope that it will be useful,
 
9
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
%  GNU General Public License for more details.
 
12
%
 
13
%  You should have received a copy of the GNU General Public License
 
14
%  along with this program; if not, write to the Free Software Foundation,
 
15
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 
16
*/
 
17
 
 
18
#include <stdio.h>
 
19
#include <stdlib.h>
 
20
#include <signal.h>
 
21
 
 
22
#include <meep.hpp>
 
23
using namespace meep;
 
24
 
 
25
double one(const vec &) { return 1.0; }
 
26
double targets(const vec &pt) {
 
27
  const double r = sqrt(pt.x()*pt.x() + pt.y()*pt.y());
 
28
  double dr = r;
 
29
  while (dr > 1) dr -= 1;
 
30
  if (dr > 0.7001) return 12.0;
 
31
  return 1.0;
 
32
}
 
33
 
 
34
#if MEEP_SINGLE
 
35
static const double tol = 1e-3, thresh = 1e-5;
 
36
#else
 
37
static const double tol = 1e-11, thresh = 1e-5;
 
38
#endif
 
39
 
 
40
int compare(double a, double b, const char *n) {
 
41
  if (fabs(a-b) > fabs(b)*tol && fabs(b) > thresh) {
 
42
    master_printf("%s differs by\t%g out of\t%g\n", n, a-b, b);
 
43
    master_printf("This gives a fractional error of %g\n", fabs(a-b)/fabs(b));
 
44
    return 0;
 
45
  } else {
 
46
    return 1;
 
47
  }
 
48
}
 
49
 
 
50
int compare_point(fields &f1, fields &f2, const vec &p) {
 
51
  monitor_point m1, m_test;
 
52
  f1.get_point(&m_test, p);
 
53
  f2.get_point(&m1, p);
 
54
  for (int i=0;i<10;i++) {
 
55
    component c = (component) i;
 
56
    if (f1.gv.has_field(c)) {
 
57
      complex<double> v1 = m_test.get_component(c), v2 = m1.get_component(c);
 
58
      if (abs(v1 - v2) > tol * abs(v2) && abs(v2) > thresh) {
 
59
        master_printf("%s differs:  %g %g out of %g %g\n",
 
60
               component_name(c), real(v2-v1), imag(v2-v1), real(v2), imag(v2));
 
61
        master_printf("This comes out to a fractional error of %g\n",
 
62
               abs(v1 - v2)/abs(v2));
 
63
        master_printf("Right now I'm looking at %g %g, time %g\n",
 
64
                      p.x(), p.y(), f1.time());
 
65
        return 0;
 
66
      }
 
67
    }
 
68
  }
 
69
  return 1;
 
70
}
 
71
 
 
72
int test_metal(double eps(const vec &), int splitting, const char *mydirname) {
 
73
  double a = 10.0;
 
74
  double ttot = 17.0;
 
75
 
 
76
  grid_volume gv = voltwo(3.0, 2.0, a);
 
77
  structure s1(gv, eps);
 
78
  structure s(gv, eps, no_pml(), identity(), splitting);
 
79
  s.set_output_directory(mydirname);
 
80
  s1.set_output_directory(mydirname);
 
81
 
 
82
  master_printf("Metal test using %d chunks...\n", splitting);
 
83
  fields f(&s);
 
84
  f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0);
 
85
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
86
  fields f1(&s1);
 
87
  f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0);
 
88
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
89
  double total_energy_check_time = 8.0;
 
90
  while (f.time() < ttot) {
 
91
    f.step();
 
92
    f1.step();
 
93
    if (!compare_point(f, f1, vec(0.5  , 0.01))) return 0;
 
94
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
 
95
    if (!compare_point(f, f1, vec(1.0  , 1.0 ))) return 0;
 
96
    if (f.time() >= total_energy_check_time) {
 
97
      if (!compare(f.total_energy(), f1.total_energy(),
 
98
                   "   total energy")) return 0;
 
99
      if (!compare(f.electric_energy_in_box(gv.surroundings()),
 
100
                   f1.electric_energy_in_box(gv.surroundings()),
 
101
                   "electric energy")) return 0;
 
102
      if (!compare(f.magnetic_energy_in_box(gv.surroundings()),
 
103
                   f1.magnetic_energy_in_box(gv.surroundings()),
 
104
                   "magnetic energy")) return 0;
 
105
      total_energy_check_time += 5.0;
 
106
    }
 
107
  }
 
108
  return 1;
 
109
}
 
110
 
 
111
int test_periodic(double eps(const vec &), int splitting, const char *mydirname) {
 
112
  double a = 10.0;
 
113
  double ttot = 17.0;
 
114
 
 
115
  grid_volume gv = voltwo(3.0, 2.0, a);
 
116
  structure s1(gv, eps);
 
117
  structure s(gv, eps, no_pml(), identity(), splitting);
 
118
  s.set_output_directory(mydirname);
 
119
  s1.set_output_directory(mydirname);
 
120
 
 
121
  master_printf("Periodic test using %d chunks...\n", splitting);
 
122
  fields f(&s);
 
123
  f.use_bloch(vec(0.1,0.7));
 
124
  f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0);
 
125
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
126
  fields f1(&s1);
 
127
  f1.use_bloch(vec(0.1,0.7));
 
128
  f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0);
 
129
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
130
  double total_energy_check_time = 8.0;
 
131
  while (f.time() < ttot) {
 
132
    f.step();
 
133
    f1.step();
 
134
    if (!compare_point(f, f1, vec(0.5  , 0.01))) return 0;
 
135
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
 
136
    if (!compare_point(f, f1, vec(1.0  , 1.0 ))) return 0;
 
137
    if (f.time() >= total_energy_check_time) {
 
138
      if (!compare(f.total_energy(), f1.total_energy(),
 
139
                   "   total energy")) return 0;
 
140
      if (!compare(f.electric_energy_in_box(gv.surroundings()),
 
141
                   f1.electric_energy_in_box(gv.surroundings()),
 
142
                   "electric energy")) return 0;
 
143
      if (!compare(f.magnetic_energy_in_box(gv.surroundings()),
 
144
                   f1.magnetic_energy_in_box(gv.surroundings()),
 
145
                   "magnetic energy")) return 0;
 
146
      total_energy_check_time += 5.0;
 
147
    }
 
148
  }
 
149
  return 1;
 
150
}
 
151
 
 
152
int test_periodic_tm(double eps(const vec &), int splitting, const char *mydirname) {
 
153
  double a = 10.0;
 
154
  double ttot = 17.0;
 
155
 
 
156
  grid_volume gv = voltwo(3.0, 2.0, a);
 
157
  structure s1(gv, eps);
 
158
  structure s(gv, eps, no_pml(), identity(), splitting);
 
159
  s.set_output_directory(mydirname);
 
160
  s1.set_output_directory(mydirname);
 
161
 
 
162
  master_printf("Periodic 2D TM test using %d chunks...\n", splitting);
 
163
  fields f(&s);
 
164
  f.use_bloch(vec(0.1,0.7));
 
165
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
166
  fields f1(&s1);
 
167
  f1.use_bloch(vec(0.1,0.7));
 
168
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
169
  double total_energy_check_time = 8.0;
 
170
  while (f.time() < ttot) {
 
171
    f.step();
 
172
    f1.step();
 
173
    if (!compare_point(f, f1, vec(0.5  , 0.01))) return 0;
 
174
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
 
175
    if (!compare_point(f, f1, vec(1.0  , 1.0 ))) return 0;
 
176
    if (f.time() >= total_energy_check_time) {
 
177
      if (!compare(f.total_energy(), f1.total_energy(),
 
178
                   "   total energy")) return 0;
 
179
      if (!compare(f.electric_energy_in_box(gv.surroundings()),
 
180
                   f1.electric_energy_in_box(gv.surroundings()),
 
181
                   "electric energy")) return 0;
 
182
      if (!compare(f.magnetic_energy_in_box(gv.surroundings()),
 
183
                   f1.magnetic_energy_in_box(gv.surroundings()),
 
184
                   "magnetic energy")) return 0;
 
185
      total_energy_check_time += 5.0;
 
186
    }
 
187
  }
 
188
  return 1;
 
189
}
 
190
 
 
191
int test_pml(double eps(const vec &), int splitting, const char *mydirname) {
 
192
  double a = 10.0;
 
193
 
 
194
  grid_volume gv = voltwo(3.0, 2.0, a);
 
195
  structure s1(gv, eps, pml(1.0, X) + pml(1.0, Y, High));
 
196
  structure s(gv, eps, pml(1.0, X) + pml(1.0, Y, High), identity(), splitting);
 
197
  s.set_output_directory(mydirname);
 
198
  s1.set_output_directory(mydirname);
 
199
 
 
200
  master_printf("Testing pml while splitting into %d chunks...\n", splitting);
 
201
  fields f(&s);
 
202
  f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,0.5), 1.0);
 
203
  f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
204
  fields f1(&s1);
 
205
  f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,0.5), 1.0);
 
206
  f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
207
  const double deltaT = 100.0;
 
208
  const double ttot = 3.1*deltaT;
 
209
  double total_energy_check_time = deltaT;
 
210
 
 
211
  while (f.time() < f.last_source_time()) f.step();
 
212
  while (f1.time() < f1.last_source_time()) f1.step();
 
213
 
 
214
  double last_energy = f.total_energy();
 
215
  while (f.time() < ttot) {
 
216
    f.step();
 
217
    f1.step();
 
218
    if (f.time() >= total_energy_check_time) {
 
219
      if (!compare_point(f, f1, vec(0.5  , 0.01))) return 0;
 
220
      if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
 
221
      if (!compare_point(f, f1, vec(1.0  , 1.0 ))) return 0;
 
222
      const double new_energy = f.total_energy();
 
223
      if (!compare(new_energy, f1.total_energy(),
 
224
                   "   total energy")) return 0;
 
225
      if (new_energy > last_energy*1e-6) {
 
226
        master_printf("Energy decaying too slowly: from %g to %g (%g)\n",
 
227
                      last_energy, new_energy, new_energy/last_energy);
 
228
        return 0;
 
229
      } else {
 
230
        master_printf("Got newE/oldE of %g\n", new_energy/last_energy);
 
231
      }
 
232
      total_energy_check_time += deltaT;
 
233
    }
 
234
  }
 
235
  return 1;
 
236
}
 
237
 
 
238
int test_pml_tm(double eps(const vec &), int splitting, const char *mydirname) {
 
239
  double a = 10.0;
 
240
 
 
241
  grid_volume gv = voltwo(3.0, 3.0, a);
 
242
  structure s1(gv, eps, pml(1.0));
 
243
  structure s(gv, eps, pml(1.0), identity(), splitting);
 
244
  s.set_output_directory(mydirname);
 
245
  s1.set_output_directory(mydirname);
 
246
 
 
247
  master_printf("Testing TM pml while splitting into %d chunks...\n", splitting);
 
248
  fields f(&s);
 
249
  f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,1.401), 1.0);
 
250
  fields f1(&s1);
 
251
  f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,1.401), 1.0);
 
252
  const double deltaT = 100.0;
 
253
  const double ttot = 3.1*deltaT;
 
254
  double total_energy_check_time = deltaT;
 
255
 
 
256
  while (f.time() < f.last_source_time()) f.step();
 
257
  while (f1.time() < f1.last_source_time()) f1.step();
 
258
 
 
259
  double last_energy = f.total_energy();
 
260
  while (f.time() < ttot) {
 
261
    f.step();
 
262
    f1.step();
 
263
    if (f.time() >= total_energy_check_time) {
 
264
      if (!compare_point(f, f1, vec(0.5  , 0.01))) return 0;
 
265
      if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
 
266
      if (!compare_point(f, f1, vec(1.0  , 1.0 ))) return 0;
 
267
      const double new_energy = f.total_energy();
 
268
      if (!compare(new_energy, f1.total_energy(),
 
269
                   "   total energy")) return 0;
 
270
      if (new_energy > last_energy*4e-6) {
 
271
        master_printf("Energy decaying too slowly: from %g to %g (%g)\n",
 
272
                      last_energy, new_energy, new_energy/last_energy);
 
273
        return 0;
 
274
      } else {
 
275
        master_printf("Got newE/oldE of %g\n", new_energy/last_energy);
 
276
      }
 
277
      total_energy_check_time += deltaT;
 
278
    }
 
279
  }
 
280
  return 1;
 
281
}
 
282
 
 
283
int test_pml_te(double eps(const vec &), int splitting, const char *mydirname) {
 
284
  double a = 10.0;
 
285
 
 
286
  grid_volume gv = voltwo(3.0, 3.0, a);
 
287
  structure s1(gv, eps, pml(1.0));
 
288
  structure s(gv, eps, pml(1.0), identity(), splitting);
 
289
  s.set_output_directory(mydirname);
 
290
  s1.set_output_directory(mydirname);
 
291
 
 
292
  master_printf("Testing TE pml while splitting into %d chunks...\n", splitting);
 
293
  fields f(&s);
 
294
  f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,1.5), 1.0);
 
295
  f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.37,1.27), 1.0);
 
296
  fields f1(&s1);
 
297
  f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,1.5), 1.0);
 
298
  f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.37,1.27), 1.0);
 
299
  const double deltaT = 100.0;
 
300
  const double ttot = 3.1*deltaT;
 
301
  double total_energy_check_time = deltaT;
 
302
 
 
303
  while (f.time() < f.last_source_time()) f.step();
 
304
  while (f1.time() < f1.last_source_time()) f1.step();
 
305
 
 
306
  double last_energy = f.total_energy();
 
307
  while (f.time() < ttot) {
 
308
    f.step();
 
309
    f1.step();
 
310
    if (f.time() >= total_energy_check_time) {
 
311
      if (!compare_point(f, f1, vec(0.5  , 0.01))) return 0;
 
312
      if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
 
313
      if (!compare_point(f, f1, vec(1.0  , 1.0 ))) return 0;
 
314
      const double new_energy = f.total_energy();
 
315
      if (!compare(new_energy, f1.total_energy(),
 
316
                   "   total energy")) return 0;
 
317
      if (new_energy > last_energy*1.1e-6) {
 
318
        master_printf("Energy decaying too slowly: from %g to %g (%g)\n",
 
319
                      last_energy, new_energy, new_energy/last_energy);
 
320
        return 0;
 
321
      } else {
 
322
        master_printf("Got newE/oldE of %g\n", new_energy/last_energy);
 
323
      }
 
324
      total_energy_check_time += deltaT;
 
325
    }
 
326
  }
 
327
  return 1;
 
328
}
 
329
 
 
330
int main(int argc, char **argv) {
 
331
  initialize mpi(argc, argv);
 
332
  quiet = true;
 
333
  const char *mydirname = "two_dimensional-out";
 
334
  trash_output_directory(mydirname);
 
335
  master_printf("Testing 2D...\n");
 
336
 
 
337
  for (int s=2;s<4;s++)
 
338
    if (!test_pml(one, s, mydirname)) abort("error in test_pml vacuum\n");
 
339
 
 
340
  for (int s=2;s<4;s++)
 
341
    if (!test_pml_tm(one, s, mydirname))
 
342
      abort("error in test_pml_tm vacuum\n");
 
343
 
 
344
  for (int s=2;s<4;s++)
 
345
    if (!test_pml_te(one, s, mydirname))
 
346
      abort("error in test_pml_te vacuum\n");
 
347
 
 
348
  for (int s=2;s<4;s++)
 
349
    if (!test_metal(one, s, mydirname)) abort("error in test_metal vacuum\n");
 
350
  //if (!test_metal(one, 200, mydirname)) abort("error in test_metal vacuum\n");
 
351
 
 
352
  for (int s=2;s<5;s++)
 
353
    if (!test_metal(targets, s, mydirname)) abort("error in test_metal targets\n");
 
354
  //if (!test_metal(targets, 60, mydirname)) abort("error in test_metal targets\n");
 
355
 
 
356
  for (int s=2;s<5;s++)
 
357
    if (!test_periodic(targets, s, mydirname))
 
358
      abort("error in test_periodic targets\n");
 
359
  //if (!test_periodic(one, 200, mydirname))
 
360
  //  abort("error in test_periodic targets\n");
 
361
 
 
362
  for (int s=2;s<4;s++)
 
363
    if (!test_periodic_tm(one, s, mydirname))
 
364
      abort("error in test_periodic_tm vacuum\n");
 
365
 
 
366
  return 0;
 
367
}