← Back to branch summary

~emmanuel-lambert/python-meep/intec

« back to all changes in this revision

Viewing changes to tests/three_d.py

  • Committer: emmanuel.lambert at ugent
  • Date: 2009-12-01 13:23:31 UTC
  • Revision ID: emmanuel.lambert@intec.ugent.be-20091201132331-4y3ad7wvx1fjnzog
merging

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/three_d.py
1
 
#!/usr/bin/env python 
 
1
#!/usr/bin/env python
2
2
#passed
3
3
from meep import *
4
4
from math import fabs,sqrt
5
5
 
 
6
quiet(True)
 
7
 
6
8
class epsilon(Callback):
7
9
    def __init__(self):
8
10
        Callback.__init__(self)
9
 
        
 
11
 
10
12
    def double_vec(self,vec):
11
13
        return self.eps(vec)
12
14
 
13
15
    def eps(self,vec):
14
 
        v = vec 
 
16
        v = vec
15
17
        r = v.x()*v.x() + v.y()*v.y()
16
18
        dr = sqrt(r)
17
19
        while dr>1:
19
21
        if dr > 0.7001:
20
22
            return 12.0
21
23
        return 1.0
22
 
        
 
24
 
23
25
 
24
26
def compare(a, b, n):
25
 
    if (fabs(a-b) > fabs(b)*2.0e-12): 
 
27
    if (fabs(a-b) > fabs(b)*2.0e-12):
26
28
        master_printf("%s differs by\t%g out of\t%g\n", n, a-b, b)
27
29
        master_printf("This gives a fractional error of %g\n", fabs(a-b)/fabs(b))
28
30
        return False
29
 
    else: 
30
 
        return True 
 
31
    else:
 
32
        return True
31
33
 
32
34
 
33
35
def compare_point(f1, f2, p):
43
45
            v1 = m_test.get_component(c)
44
46
            v2 = m1.get_component(c)
45
47
 
46
 
            if (abs(v1 - v2) > 0.0*2e-15*abs(v2)): 
 
48
            if (abs(v1 - v2) > 0.0*2e-15*abs(v2)):
47
49
                master_printf("%s differs:  %g %g out of %g %g\n", component_name(c), real(v2-v1), imag(v2-v1), real(v2), imag(v2))
48
50
                master_printf("This comes out to a fractional error of %g\n", abs(v1 - v2)/abs(v2))
49
51
                master_printf("Right now I'm looking at %g %g, time %g\n", p.x(), p.y(), p.z(), f1.time())
63
65
            v1 = m_test.get_component(c)
64
66
            v2 = m1.get_component(c)
65
67
 
66
 
            if (abs(v1 - v2) > 1e-14 *  (5.0 + abs(v2))): 
 
68
            if (abs(v1 - v2) > 1e-14 *  (5.0 + abs(v2))):
67
69
                master_printf("%s differs:  %g %g out of %g %g\n", component_name(c), real(v2-v1), imag(v2-v1), real(v2), imag(v2))
68
70
                master_printf("This comes out to a fractional error of %g\n", abs(v1 - v2)/abs(v2))
69
71
                master_printf("Right now I'm looking at %g %g, time %g\n", p.x(), p.y(), p.z(), f1.time())
76
78
    v = vol3d(1.5, 1.0, 1.2, a)
77
79
    s=structure(v, eps, pml(0.401))
78
80
    master_printf("Testing pml quality...\n");
79
 
  
 
81
 
80
82
    f=fields(s)
81
83
    f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.751,0.5,0.601), 1.0);
82
84
    deltaT = 10.0
90
92
 
91
93
    while (f.time() < ttot):
92
94
        f.step()
93
 
        
 
95
 
94
96
        if (f.time() >= total_energy_check_time):
95
97
            new_energy = f.total_energy()
96
98
 
101
103
            else:
102
104
                master_printf("Got newE/oldE of %g\n", new_energy/last_energy);
103
105
                #print "Got newE/oldE of ", new_energy/last_energy, " \n"
104
 
      
 
106
 
105
107
            total_energy_check_time += deltaT;
106
 
  
107
 
    return True 
 
108
 
 
109
    return True
108
110
 
109
111
 
110
112
def test_pml_splitting(eps, splitting):
119
121
    f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.099,0.499,0.501), 1.0)
120
122
    ttot = 31.0
121
123
    next_energy_time = 10.0
122
 
    while (f.time() < ttot): 
 
124
    while (f.time() < ttot):
123
125
        f.step()
124
126
        f1.step()
125
127
        if not approx_point(f, f1, vec(0.5  , 0.01 , 1.0 )):
126
128
            return False
127
 
        if not approx_point(f, f1, vec(0.46 , 0.33 , 0.33)): 
 
129
        if not approx_point(f, f1, vec(0.46 , 0.33 , 0.33)):
128
130
            return False
129
131
        if not approx_point(f, f1, vec(1.0  , 1.0  , 0.33)):
130
132
            return False
134
136
            if not compare(f.total_energy(), f1.total_energy(), "   total energy"):
135
137
                return False
136
138
            next_energy_time += 10.0;
137
 
     
 
139
 
138
140
    return 1
139
141
 
140
142
 
151
153
    f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.299,0.401), 1.0)
152
154
    f1=fields(s1)
153
155
    f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299,0.299,0.401), 1.0)
154
 
    
 
156
 
155
157
    total_energy_check_time = 8.0
156
 
 
 
158
 
157
159
    while (f.time() < ttot):
158
160
        f.step()
159
161
        f1.step()
163
165
                return False
164
166
        if not compare_point(f, f1, vec(1.301  , 0.301 ,  0.399 )):
165
167
                return False
166
 
        
 
168
 
167
169
        if (f.time() >= total_energy_check_time):
168
170
            if not compare(f.total_energy(), f1.total_energy(),"   total energy"):
169
171
                return False
171
173
                return False
172
174
            if not compare(f.magnetic_energy_in_box(v.surroundings()), f1.magnetic_energy_in_box(v.surroundings()), "magnetic energy"):
173
175
                return False
174
 
 
 
176
 
175
177
            total_energy_check_time += 5.0
176
 
  
 
178
 
177
179
    return True
178
180
 
179
181
 
192
194
    f1.use_bloch(vec(0.1,0.7,0.3))
193
195
    f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.3,0.25,0.5), 1.0)
194
196
    total_energy_check_time = 8.0
195
 
 
 
197
 
196
198
    while (f.time() < ttot):
197
199
        f.step()
198
200
        f1.step()
202
204
                return False
203
205
        if not compare_point(f, f1, vec(1.0  , 0.25 ,  0.301 )):
204
206
                return False
205
 
        
 
207
 
206
208
        if (f.time() >= total_energy_check_time):
207
209
            if not compare(f.total_energy(), f1.total_energy(),"   total energy"):
208
210
                return False
210
212
                return False
211
213
            if not compare(f.magnetic_energy_in_box(v.surroundings()), f1.magnetic_energy_in_box(v.surroundings()), "magnetic energy"):
212
214
                return False
213
 
 
 
215
 
214
216
            total_energy_check_time += 5.0
215
 
  
 
217
 
216
218
    return True
217
219
 
218
220