← Back to branch summary

~emmanuel-lambert/python-meep/intec

« back to all changes in this revision

Viewing changes to tests-intec/two_dimensional.py

  • Committer: emmanuel.lambert at ugent
  • Date: 2009-11-23 10:23:17 UTC
  • Revision ID: emmanuel.lambert@intec.ugent.be-20091123102317-m00vrwyw0m6kpcmm
further merging with NS / replacing tests dir with default tests and creating new 'tests-intec' subdir

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests-intec/two_dimensional.py
 
1
#!/usr/bin/env python 
 
2
#passed
 
3
from meep import *
 
4
from math import fabs,sqrt
 
5
 
 
6
class epsilon(Callback):
 
7
    def __init__(self):
 
8
        Callback.__init__(self)
 
9
    def double_vec(self,vec):
 
10
        return self.eps(vec)
 
11
    def eps(self,vec):
 
12
        v = vec 
 
13
        r = v.x()*v.x() + v.y()*v.y()
 
14
        dr = sqrt(r)
 
15
        while dr>1:
 
16
            dr-=1
 
17
        if dr > 0.7001:
 
18
            return 12.0
 
19
        return 1.0
 
20
        
 
21
 
 
22
def compare(a, b, n):
 
23
    if (fabs(a-b) > fabs(b)*1.3e-14): 
 
24
        master_printf("%s differs by\t%g out of\t%g\n", n, a-b, b)
 
25
        master_printf("This gives a fractional error of %g\n", fabs(a-b)/fabs(b))
 
26
        return False
 
27
    else: 
 
28
        return True 
 
29
 
 
30
 
 
31
def compare_point(f1, f2, p):
 
32
    m1=monitor_point()
 
33
    m_test=monitor_point()
 
34
 
 
35
    f1.get_point(m_test, p)
 
36
    f2.get_point(m1, p)
 
37
 
 
38
    for i in range(0,10):
 
39
        c = i
 
40
        if (f1.gv.has_field(c)):
 
41
            v1 = m_test.get_component(c)
 
42
            v2 = m1.get_component(c)
 
43
 
 
44
            if (abs(v1 - v2) > 0.0*2e-15*abs(v2)): 
 
45
                master_printf("%s differs:  %g %g out of %g %g\n", component_name(c), (v2-v1).real, (v2-v1).imag, v2.real, v2.imag)
 
46
                master_printf("This comes out to a fractional error of %g\n", abs(v1 - v2)/abs(v2))
 
47
                master_printf("Right now I'm looking at %g %g, time %g\n", p.x(), p.y(), f1.time())
 
48
                return False
 
49
    return True
 
50
 
 
51
 
 
52
 
 
53
def test_pml(epsilon,splitting):
 
54
    a = 10.0
 
55
    v = voltwo(3.0, 2.0, a)
 
56
 
 
57
    s1=structure(v, epsilon, pml(1.0, X) + pml(1.0, Y, High))
 
58
    s=structure(v, epsilon, pml(1.0, X) + pml(1.0, Y, High), identity(), splitting)
 
59
 
 
60
    master_printf( "Testing pml while splitting into %d chunks...\n", splitting)
 
61
 
 
62
    f=fields(s)
 
63
    f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,0.5), 1.0);
 
64
    f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
65
    f1=fields(s1);
 
66
    f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,0.5), 1.0);
 
67
    f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,0.401), 1.0);
 
68
 
 
69
    deltaT = 100.0;
 
70
    ttot = 3.1*deltaT;
 
71
    total_energy_check_time = deltaT;
 
72
 
 
73
    while (f.time() < f.last_source_time()):
 
74
        f.step()
 
75
    while (f1.time() < f1.last_source_time()):
 
76
        f1.step()
 
77
 
 
78
    last_energy = f.total_energy()
 
79
 
 
80
    while (f.time() < ttot):
 
81
        f.step()
 
82
        f1.step()
 
83
        
 
84
        if (f.time() >= total_energy_check_time):
 
85
            if not compare_point(f, f1, vec(0.5  , 0.01)):
 
86
                return False
 
87
            if not compare_point(f, f1, vec(0.46 , 0.33)):
 
88
                return False
 
89
            if not compare_point(f, f1, vec(1.0  , 1.0 )):
 
90
                return False
 
91
            new_energy = f.total_energy()
 
92
 
 
93
            if not compare(new_energy, f1.total_energy(),"   total energy"):
 
94
                return False
 
95
 
 
96
            if new_energy > last_energy*1e-6:
 
97
                master_printf("Energy decaying too slowly: from %g to %g (%g)\n", last_energy, new_energy, new_energy/last_energy)
 
98
                return False
 
99
            else:
 
100
                #print 'Energies:', new_energy, last_energy, new_energy/last_energy
 
101
                master_printf("Got newE/oldE of %g\n", new_energy/last_energy);
 
102
     
 
103
            total_energy_check_time += deltaT;
 
104
  
 
105
    return True 
 
106
 
 
107
 
 
108
def test_pml_tm(epsilon,splitting):
 
109
    a = 10.0
 
110
    v = voltwo(3.0, 3.0, a)
 
111
 
 
112
    s1=structure(v, epsilon, pml(1.0))
 
113
    s=structure(v, epsilon, pml(1.0), identity(), splitting)
 
114
    master_printf("Testing TM pml while splitting into %d chunks...\n", splitting);
 
115
 
 
116
    f=fields(s);
 
117
    f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,1.401), 1.0)
 
118
    f1=fields(s1);
 
119
    f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299,1.401), 1.0)
 
120
 
 
121
    deltaT = 100.0;
 
122
    ttot = 3.1*deltaT;
 
123
    total_energy_check_time = deltaT;
 
124
 
 
125
    while (f.time() < f.last_source_time()):
 
126
        f.step()
 
127
    while (f1.time() < f1.last_source_time()):
 
128
        f1.step()
 
129
 
 
130
    last_energy = f.total_energy()
 
131
 
 
132
    while (f.time() < ttot):
 
133
        f.step()
 
134
        f1.step()
 
135
        
 
136
        if (f.time() >= total_energy_check_time):
 
137
            if not compare_point(f, f1, vec(0.5  , 0.01)):
 
138
                return False
 
139
            if not compare_point(f, f1, vec(0.46 , 0.33)):
 
140
                return False
 
141
            if not compare_point(f, f1, vec(1.0  , 1.0 )):
 
142
                return False
 
143
            new_energy = f.total_energy()
 
144
 
 
145
            if not compare(new_energy, f1.total_energy(),"   total energy"):
 
146
                return False
 
147
 
 
148
            if new_energy > last_energy*4e-6:
 
149
                master_printf("Energy decaying too slowly: from %g to %g (%g)\n", last_energy, new_energy, new_energy/last_energy)
 
150
                return False
 
151
            else:
 
152
                #print 'Energies:', new_energy, last_energy, new_energy/last_energy
 
153
                master_printf("Got newE/oldE of %g\n", new_energy/last_energy);
 
154
    
 
155
            total_energy_check_time += deltaT;
 
156
  
 
157
    return True
 
158
 
 
159
 
 
160
 
 
161
def test_pml_te(epsilon,splitting):
 
162
    a = 10.0
 
163
    v = voltwo(3.0, 3.0, a)
 
164
 
 
165
    s1=structure(v, epsilon, pml(1.0))
 
166
    s=structure(v, epsilon, pml(1.0), identity(), splitting)
 
167
 
 
168
    master_printf("Testing TE pml while splitting into %d chunks...\n", splitting);
 
169
 
 
170
    f=fields(s);
 
171
    f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,1.5), 1.0)
 
172
    f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.37,1.27), 1.0)
 
173
 
 
174
    f1=fields(s1);
 
175
    f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5,1.5), 1.0)
 
176
    f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.37,1.27), 1.0)
 
177
 
 
178
    deltaT = 100.0;
 
179
    ttot = 3.1*deltaT;
 
180
    total_energy_check_time = deltaT;
 
181
 
 
182
    while (f.time() < f.last_source_time()):
 
183
        f.step()
 
184
    while (f1.time() < f1.last_source_time()):
 
185
        f1.step()
 
186
 
 
187
    last_energy = f.total_energy()
 
188
 
 
189
    while (f.time() < ttot):
 
190
        f.step()
 
191
        f1.step()
 
192
        
 
193
        if (f.time() >= total_energy_check_time):
 
194
            if not compare_point(f, f1, vec(0.5  , 0.01)):
 
195
                return False
 
196
            if not compare_point(f, f1, vec(0.46 , 0.33)):
 
197
                return False
 
198
            if not compare_point(f, f1, vec(1.0  , 1.0 )):
 
199
                return False
 
200
            new_energy = f.total_energy()
 
201
 
 
202
            if not compare(new_energy, f1.total_energy(),"   total energy"):
 
203
                return False
 
204
 
 
205
            if new_energy > last_energy*1.1e-6:
 
206
                master_printf("Energy decaying too slowly: from %g to %g (%g)\n", last_energy, new_energy, new_energy/last_energy)
 
207
                return False
 
208
            else:
 
209
                master_printf("Got newE/oldE of %g\n", new_energy/last_energy);
 
210
                      
 
211
            total_energy_check_time += deltaT;
 
212
    
 
213
    return True
 
214
 
 
215
def test_metal(epsilon,splitting):
 
216
    a = 10.0
 
217
    ttot = 17.0
 
218
    v = voltwo(3.0, 3.0, a)
 
219
    #no_pml = boundary_region()
 
220
 
 
221
    s1=structure(v, epsilon)
 
222
    s=structure(v, epsilon, no_pml(), identity(), splitting)
 
223
    master_printf("Metal test using %d chunks...\n", splitting)
 
224
 
 
225
    f=fields(s)
 
226
    f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0)
 
227
    f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0)
 
228
    f1=fields(s1)
 
229
    f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0)
 
230
    f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0)
 
231
    
 
232
    total_energy_check_time = 8.0
 
233
 
 
234
    while (f.time() < ttot):
 
235
        f.step()
 
236
        f1.step()
 
237
        if not compare_point(f, f1, vec(0.5  , 0.01)):
 
238
                return False
 
239
        if not compare_point(f, f1, vec(0.46 , 0.33)):
 
240
                return False
 
241
        if not compare_point(f, f1, vec(1.0  , 1.0 )):
 
242
                return False
 
243
        
 
244
        if (f.time() >= total_energy_check_time):
 
245
            if not compare(f.total_energy(), f1.total_energy(),"   total energy"):
 
246
                return False
 
247
            if not compare(f.electric_energy_in_box(v.surroundings()), f1.electric_energy_in_box(v.surroundings()), "electric energy"):
 
248
                return False
 
249
            if not compare(f.magnetic_energy_in_box(v.surroundings()), f1.magnetic_energy_in_box(v.surroundings()), "magnetic energy"):
 
250
                return False
 
251
 
 
252
            total_energy_check_time += 5.0
 
253
  
 
254
    return True
 
255
 
 
256
def test_periodic(epsilon,splitting):
 
257
    a = 10.0
 
258
    ttot = 17.0
 
259
    v = voltwo(3.0, 2.0, a)
 
260
    #no_pml = boundary_region()
 
261
 
 
262
    s1=structure(v, epsilon)
 
263
    s=structure(v, epsilon, no_pml(), identity(), splitting)
 
264
    master_printf("Periodic test using %d chunks...\n", splitting);
 
265
 
 
266
    f=fields_complex(s)
 
267
    f.use_bloch(vec(0.1,0.7))
 
268
    f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0)
 
269
    f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0)
 
270
    f1=fields_complex(s1)
 
271
    f1.use_bloch(vec(0.1,0.7))
 
272
    f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.5), 1.0)
 
273
    f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0)
 
274
    
 
275
    total_energy_check_time = 8.0
 
276
 
 
277
    while (f.time() < ttot):
 
278
        f.step()
 
279
        f1.step()
 
280
        if not compare_point(f, f1, vec(0.5  , 0.01)):
 
281
                return False
 
282
        if not compare_point(f, f1, vec(0.46 , 0.33)):
 
283
                return False
 
284
        if not compare_point(f, f1, vec(1.0  , 1.0 )):
 
285
                return False
 
286
        
 
287
        if (f.time() >= total_energy_check_time):
 
288
            if not compare(f.total_energy(), f1.total_energy(),"   total energy"):
 
289
                return False
 
290
            if not compare(f.electric_energy_in_box(v.surroundings()), f1.electric_energy_in_box(v.surroundings()), "electric energy"):
 
291
                return False
 
292
            if not compare(f.magnetic_energy_in_box(v.surroundings()), f1.magnetic_energy_in_box(v.surroundings()), "magnetic energy"):
 
293
                return False
 
294
 
 
295
            total_energy_check_time += 5.0
 
296
  
 
297
    return True
 
298
 
 
299
def test_periodic_tm(epsilon,splitting):
 
300
    a = 10.0
 
301
    ttot = 17.0
 
302
    v = voltwo(3.0, 2.0, a)
 
303
    #no_pml = boundary_region()
 
304
 
 
305
    s1=structure(v, epsilon)
 
306
    s=structure(v, epsilon, no_pml(), identity(), splitting)
 
307
    master_printf("Periodic 2D TM test using %d chunks...\n", splitting);
 
308
 
 
309
    f=fields_complex(s)
 
310
    f.use_bloch(vec(0.1,0.7))
 
311
    f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0)
 
312
    f1=fields_complex(s1)
 
313
    f1.use_bloch(vec(0.1,0.7))
 
314
    f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.401), 1.0)
 
315
    
 
316
    total_energy_check_time = 8.0
 
317
 
 
318
    while (f.time() < ttot):
 
319
        f.step()
 
320
        f1.step()
 
321
        if not compare_point(f, f1, vec(0.5  , 0.01)):
 
322
                return False
 
323
        if not compare_point(f, f1, vec(0.46 , 0.33)):
 
324
                return False
 
325
        if not compare_point(f, f1, vec(1.0  , 1.0 )):
 
326
                return False
 
327
        
 
328
        if (f.time() >= total_energy_check_time):
 
329
            if not compare(f.total_energy(), f1.total_energy(),"   total energy"):
 
330
                return False
 
331
            if not compare(f.electric_energy_in_box(v.surroundings()), f1.electric_energy_in_box(v.surroundings()), "electric energy"):
 
332
                return False
 
333
            if not compare(f.magnetic_energy_in_box(v.surroundings()), f1.magnetic_energy_in_box(v.surroundings()), "magnetic energy"):
 
334
                return False
 
335
 
 
336
            total_energy_check_time += 5.0
 
337
  
 
338
    return True
 
339
 
 
340
 
 
341
 
 
342
master_printf('Testing 2D ....\n')
 
343
 
 
344
 
 
345
for s in range(2,4):
 
346
    if not test_pml(EPS,s):
 
347
        abort("error in test_pml vacuum\n")
 
348
 
 
349
for s in range(2,4):
 
350
    if not test_pml_tm(EPS,s):
 
351
        abort("error in test_pml TM vacuum\n")
 
352
 
 
353
for s in range(2,4):
 
354
    if not test_pml_te(EPS,s):
 
355
        abort("error in test_pml TE vacuum\n")
 
356
 
 
357
for s in range(2,4):
 
358
    if not test_metal(EPS,s):
 
359
        abort("error in test_metal vacuum\n")
 
360
 
 
361
trg = epsilon()
 
362
set_EPS_Callback(trg.__disown__())
 
363
 
 
364
for s in range(2,5):
 
365
    if not test_metal(EPS,s):
 
366
        abort("error in test_metal targets \n")
 
367
 
 
368
for s in range(2,5):
 
369
    if not test_periodic(EPS,s):
 
370
        abort("error in test_periodic targets \n")
 
371
 
 
372
del_EPS_Callback()
 
373
 
 
374
for s in range(2,4):
 
375
    if not test_periodic_tm(EPS,s):
 
376
        abort("error in test_periodic TM targets \n")
 
377
 
 
378
 
 
379