Vlasiator ebf0dd394 on dev (v5.4.0 + 1054 commits)
Loading...
Searching...
No Matches
vectorpotentialdipole_verify.py
Go to the documentation of this file.
1#!/usr/bin/env python import matplotlib.pyplot as plt
2
3# /*
4# * This file is part of Vlasiator.
5# * Copyright 2010-2016 Finnish Meteorological Institute
6# * Copyright 2017-2019 University of Helsinki
7# *
8# * For details of usage, see the COPYING file and read the "Rules of the Road"
9# * at http://www.physics.helsinki.fi/vlasiator/
10# *
11# * This program is free software; you can redistribute it and/or modify
12# * it under the terms of the GNU General Public License as published by
13# * the Free Software Foundation; either version 2 of the License, or
14# * (at your option) any later version.
15# *
16# * This program is distributed in the hope that it will be useful,
17# * but WITHOUT ANY WARRANTY; without even the implied warranty of
18# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19# * GNU General Public License for more details.
20# *
21# * You should have received a copy of the GNU General Public License along
22# * with this program; if not, write to the Free Software Foundation, Inc.,
23# * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24# */
25import numpy as np
26import math
27import sys,os
28import pytools as pt
29import matplotlib.pyplot as plt
30import fieldmodels
31
32''' Testing routine for different dipole formulations
33
34 run using "python vectorpotentialdipole_verify.py [arg1] [arg2]" where arg1 is a number from 0 to 4
35 for different test profile starting positions, directions, and dipole tilts.
36
37 If arg2 is present, the code also calculates verification of derivative terms.
38
39 Generates plots of magnetic field components for different dipole models along cuts through the
40 simulation domain. For derivative analysis, calculates the derivatives along said cuts analytically
41 and numerically and outputs the ratio.
42
43'''
44
45if len(sys.argv)!=1:
46 testset = int(sys.argv[1])
47else:
48 testset = 0
49
50if len(sys.argv)!=2:
51 calcderivatives=True
52else:
53 calcderivatives=False
54
55plotmagnitude=False
56
57plt.switch_backend('Agg')
58
59outfilename = "./vecpotdip_verify_"+str(testset)+".png"
60
61RE=6371000.
62#RE=1
63epsilon=1.e-15
64
65if testset==0:
66 tilt_angle_phi = 0.
67 tilt_angle_theta = 0.
68 line_theta = 0.
69 line_start = np.array([0,0,0])
70elif testset==1:
71 tilt_angle_phi = 0.
72 tilt_angle_theta = 0.
73 line_theta = 45.
74 line_start = np.array([0,0,0])
75elif testset==2:
76 tilt_angle_phi = 0.
77 tilt_angle_theta = 0.
78 line_theta = 0.
79 line_start = np.array([-3,-3,-3])
80elif testset==3:
81 tilt_angle_phi = 0.
82 tilt_angle_theta = 0.
83 line_theta = 45.
84 line_start = np.array([-3,-3,-3])
85elif testset==4:
86 tilt_angle_phi = 10
87 tilt_angle_theta = 0.
88 line_theta = 0.
89 line_start = np.array([0,0,0])
90elif testset==5:
91 tilt_angle_phi = 10
92 tilt_angle_theta = 45.
93 line_theta = 0.
94 line_start = np.array([0,0,0])
95else: # Same as 0
96 print("Default")
97 tilt_angle_phi = 0.
98 tilt_angle_theta = 0.
99 line_theta = 0.
100 line_start = np.array([0,0,0])
101
102print("Test set "+str(testset)+" line start "+str(line_start)+" tilt phi "+str(tilt_angle_phi)+" tilt theta "+str(tilt_angle_theta)+" line theta "+str(line_theta))
103
104line_phi = np.array([0,45,80,90,110,135])*math.pi/180.
105line_theta = np.zeros(len(line_phi)) + line_theta * math.pi/180.
106#line_start = np.array([-5,-5,-5])
107step = 0.1
108
109linewidth=2
110linthresh=1.e-10
111fontsize=20
112
113#fieldmodels.dipole.set_dipole(centerx, centery, centerz, tilt_phi, tilt_theta, mult=1.0, radius_f=None, radius_z=None):
114dip = fieldmodels.dipole(0,0,0,tilt_angle_phi,tilt_angle_theta)
115mdip = fieldmodels.dipole(80*RE,0,0,tilt_angle_phi,180.-tilt_angle_theta)
116
117# Create figure
118fig = plt.figure()
119fig.set_size_inches(20,30)
120nsubplots=len(line_theta)
121for i in range(nsubplots):
122 fig.add_subplot(nsubplots,1,i+1)
123axes = fig.get_axes()
124
125radii = np.arange(0.1,100,step)*RE
126nr=len(radii)
127radiiRE = radii/RE
128
129fig.suptitle(r"Profiles starting from ("+str(line_start[0])+","+str(line_start[1])+","+str(line_start[2])+") [RE] with dipole tilt $\Phi="+str(int(tilt_angle_phi))+"$, $\Theta="+str(int(tilt_angle_theta))+"$", fontsize=fontsize)
130
131for i in range(nsubplots):
132 print("subplot ",i)
133 ax = axes[i]
134
135 ax.text(0.2,0.08,r"profile with $\theta="+str(int(line_theta[i]*180./math.pi))+"$, $\phi="+str(int(line_phi[i]*180./math.pi))+"$",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
136
137 xv = line_start[0]*RE + radii*np.sin(line_phi[i])*np.cos(line_theta[i])
138 yv = line_start[1]*RE + radii*np.sin(line_phi[i])*np.sin(line_theta[i])
139 zv = line_start[2]*RE + radii*np.cos(line_phi[i])
140
141 B1 = np.zeros([nr,4]) # X-scaled vector dipole
142 B2 = np.zeros([nr,4]) # regular dipole
143 B3 = np.zeros([nr,4]) # regular dipole + mirror dipole
144 B4 = np.zeros([nr,4]) # line dipole + mirror dipole
145
146 for j in range(nr):
147 for k in range(3):
148 B1[j,k] = dip.getX(xv[j],yv[j],zv[j],0,k,0)
149 B2[j,k] = dip.get_old(xv[j],yv[j],zv[j],0,k,0)
150 B3[j,k] = B2[j,k] + mdip.get_old(xv[j],yv[j],zv[j],0,k,0)
151 B4[j,k] = dip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
152 B4[j,k] = B4[j,k] + mdip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
153 if plotmagnitude is True:
154 B1[j,3] = np.linalg.norm(B1[j,0:3])
155 B2[j,3] = np.linalg.norm(B2[j,0:3])
156 B3[j,3] = np.linalg.norm(B3[j,0:3])
157 B4[j,3] = np.linalg.norm(B4[j,0:3])
158
159 colors=['r','k','b','magenta']
160 coords = ['x','y','z','mag']
161
162 plotrange = range(3)
163 if plotmagnitude is True:
164 plotrange = range(4)
165 for k in plotrange:
166 ax.plot(radiiRE, B1[:,k], c=colors[k], linestyle='-', linewidth=linewidth, label='vectorpot B'+coords[k], zorder=-10)
167 ax.plot(radiiRE, B2[:,k], c=colors[k], linestyle='--', linewidth=linewidth, label='regular B'+coords[k])
168 ax.plot(radiiRE, B3[:,k], c=colors[k], linestyle=':', linewidth=linewidth, label='reg+mirror B'+coords[k])
169 if tilt_angle_phi<epsilon:
170 ax.plot(radiiRE, B4[:,k], c=colors[k], linestyle='-.', linewidth=linewidth, label='line+mirror B'+coords[k])
171
172 ax.set_xlabel(r"$r$ [$r_\mathrm{E}$]", fontsize=fontsize)
173 ax.set_xlim([0,70])
174 #ax.set_yscale('log', nonposy='clip')
175 ax.set_yscale('symlog', linthreshy=linthresh)
176 for item in ax.get_xticklabels():
177 item.set_fontsize(fontsize)
178 for item in ax.get_yticklabels():
179 item.set_fontsize(fontsize)
180
181 ylims = np.array(ax.get_ylim())
182 if ylims[0] < -1e-4:
183 ylims[0] = -1e-4
184 if ylims[1] > 1e-4:
185 ylims[1] = 1e-4
186 ax.set_ylim(ylims)
187
188handles, labels = axes[-1].get_legend_handles_labels()
189axes[-1].legend(handles, labels, fontsize=fontsize)
190
191fig.savefig(outfilename)
192plt.close()
193
194
195
196if calcderivatives:
197 # Derivatives
198 step2=0.00001 # distance in each direction for calculating numerical derivative
199 for kkk in range(3):
200 print("derivatives d"+coords[kkk])
201 dB1 = np.zeros([nr,3,3])
202 dB2 = np.zeros([nr,3,3])
203 dB3 = np.zeros([nr,3,3])
204 dB4 = np.zeros([nr,3,3])
205
206 # Create figure
207 fig = plt.figure()
208 fig.set_size_inches(20,30)
209 for i in range(nsubplots):
210 fig.add_subplot(nsubplots,1,i+1)
211 axes = fig.get_axes()
212
213 fig.suptitle(r"Numerical and analytical derivative ratios, profiles starting from ("+str(line_start[0])+","+str(line_start[1])+","+str(line_start[2])+") [RE] with dipole tilt $\Phi="+str(int(tilt_angle_phi))+"$, $\Theta="+str(int(tilt_angle_theta))+"$", fontsize=fontsize)
214
215 for i in range(nsubplots):
216 print("derivatives subplot ",i)
217 ax = axes[i]
218
219 xv = line_start[0]*RE + radii*np.sin(line_phi[i])*np.cos(line_theta[i])
220 yv = line_start[1]*RE + radii*np.sin(line_phi[i])*np.sin(line_theta[i])
221 zv = line_start[2]*RE + radii*np.cos(line_phi[i])
222
223 for j in range(nr):
224 for k in range(3):
225 B1[j,k] = dip.getX(xv[j],yv[j],zv[j],0,k,0)
226 B2[j,k] = dip.get_old(xv[j],yv[j],zv[j],0,k,0)
227 B3[j,k] = B2[j,k] + mdip.get_old(xv[j],yv[j],zv[j],0,k,0)
228 B4[j,k] = dip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
229 B4[j,k] = B4[j,k] + mdip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
230 #for kk in range(3):
231 kk=kkk
232 dB1[j,k,kk] = dip.getX(xv[j],yv[j],zv[j],1,k,kk)
233 dB2[j,k,kk] = dip.get_old(xv[j],yv[j],zv[j],1,k,kk)
234 dB3[j,k,kk] = dB2[j,k,kk] + mdip.get_old(xv[j],yv[j],zv[j],1,k,kk)
235 dB4[j,k,kk] = dip.get_ldp(xv[j],yv[j],zv[j],1,k,kk)
236 dB4[j,k,kk] = dB4[j,k,kk] + mdip.get_ldp(xv[j],yv[j],zv[j],1,k,kk)
237
238 # analytical derivative vs numerical derivative
239 for j in np.arange(1,nr-1):
240 for k in range(3):
241
242 # d/dx
243 if kkk==0:
244 cdbx=(dip.getX(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.getX(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
245 if abs(cdbx) > epsilon*B1[j,k]:
246 dB1[j,k,0] = dB1[j,k,0]/cdbx
247 elif (abs(cdbx)<epsilon*B1[j,k]) and (abs(dB1[j,k,0])<epsilon*B1[j,k]):
248 dB1[j,k,0] = None
249 else:
250 dB1[j,k,0] = 0
251 if abs(B1[j,k])<epsilon:
252 dB1[j,k,0] = None
253
254 cdbx=(dip.get_old(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.get_old(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
255 if abs(cdbx) > epsilon*B2[j,k]:
256 dB2[j,k,0] = dB2[j,k,0]/cdbx
257 elif (abs(cdbx)<epsilon*B2[j,k]) and (abs(dB2[j,k,0])<epsilon*B2[j,k]):
258 dB2[j,k,0] = None
259 else:
260 dB2[j,k,0] = 0
261 if abs(B2[j,k])<epsilon:
262 dB2[j,k,0] = None
263
264 cdbx=(dip.get_old(xv[j]+step2*RE,yv[j],zv[j],0,k,0)+mdip.get_old(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.get_old(xv[j]-step2*RE,yv[j],zv[j],0,k,0) - mdip.get_old(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
265 if abs(cdbx) > epsilon*B3[j,k]:
266 dB3[j,k,0] = dB3[j,k,0]/cdbx
267 elif (abs(cdbx)<epsilon*B3[j,k]) and (abs(dB3[j,k,0])<epsilon*B3[j,k]):
268 dB3[j,k,0] = None
269 else:
270 dB3[j,k,0] = 0
271 if abs(B3[j,k])<epsilon:
272 dB3[j,k,0] = None
273
274 cdbx=(dip.get_ldp(xv[j]+step2*RE,yv[j],zv[j],0,k,0)+mdip.get_ldp(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.get_ldp(xv[j]-step2*RE,yv[j],zv[j],0,k,0) - mdip.get_ldp(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
275 if abs(cdbx) > epsilon*B4[j,k]:
276 dB4[j,k,0] = dB4[j,k,0]/cdbx
277 elif (abs(cdbx)<epsilon*B4[j,k]) and (abs(dB4[j,k,0])<epsilon*B4[j,k]):
278 dB4[j,k,0] = None
279 else:
280 dB4[j,k,0] = 0
281 if abs(B4[j,k])<epsilon:
282 dB4[j,k,0] = None
283
284 # d/dy
285 if kkk==1:
286 cdby=(dip.getX(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.getX(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
287 if abs(cdby) > epsilon*B1[j,k]:
288 dB1[j,k,1] = dB1[j,k,1]/cdby
289 elif (abs(cdby)<epsilon*B1[j,k]) and (abs(dB1[j,k,1])<epsilon*B1[j,k]):
290 dB1[j,k,1] = None
291 else:
292 dB1[j,k,1] = 0
293 if abs(B1[j,k])<epsilon:
294 dB1[j,k,1] = None
295
296 cdby=(dip.get_old(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.get_old(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
297 if abs(cdby) > epsilon*B2[j,k]:
298 dB2[j,k,1] = dB2[j,k,1]/cdby
299 elif (abs(cdby)<epsilon*B2[j,k]) and (abs(dB2[j,k,1])<epsilon*B2[j,k]):
300 dB2[j,k,1] = None
301 else:
302 dB2[j,k,1] = 0
303 if abs(B2[j,k])<epsilon:
304 dB2[j,k,1] = None
305
306 cdby=(dip.get_old(xv[j],yv[j]+step2*RE,zv[j],0,k,0)+mdip.get_old(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.get_old(xv[j],yv[j]-step2*RE,zv[j],0,k,0) - mdip.get_old(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
307 if abs(cdby) > epsilon*B3[j,k]:
308 dB3[j,k,1] = dB3[j,k,1]/cdby
309 elif (abs(cdby)<epsilon*B3[j,k]) and (abs(dB3[j,k,1])<epsilon*B3[j,k]):
310 dB3[j,k,1] = None
311 else:
312 dB3[j,k,1] = 0
313 if abs(B3[j,k])<epsilon:
314 dB3[j,k,1] = None
315
316 cdby=(dip.get_ldp(xv[j],yv[j]+step2*RE,zv[j],0,k,0)+mdip.get_ldp(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.get_ldp(xv[j],yv[j]-step2*RE,zv[j],0,k,0) - mdip.get_ldp(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
317 if abs(cdby) > epsilon*B4[j,k]:
318 dB4[j,k,1] = dB4[j,k,1]/cdby
319 elif (abs(cdby)<epsilon*B4[j,k]) and (abs(dB4[j,k,1])<epsilon*B4[j,k]):
320 dB4[j,k,1] = None
321 else:
322 dB4[j,k,1] = 0
323 if abs(B4[j,k])<epsilon:
324 dB4[j,k,1] = None
325
326 # d/dz
327 if kkk==2:
328 cdbz=(dip.getX(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.getX(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
329 if abs(cdbz) > epsilon*B1[j,k]:
330 dB1[j,k,2] = dB1[j,k,2]/cdbz
331 elif (abs(cdbz)<epsilon*B1[j,k]) and (abs(dB1[j,k,2])<epsilon*B1[j,k]):
332 dB1[j,k,2] = None
333 else:
334 dB1[j,k,2] = 0
335 if abs(B1[j,k])<epsilon:
336 dB1[j,k,2] = None
337
338 cdbz=(dip.get_old(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.get_old(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
339 if abs(cdbz) > epsilon*B2[j,k]:
340 dB2[j,k,2] = dB2[j,k,2]/cdbz
341 elif (abs(cdbz)<epsilon*B2[j,k]) and (abs(dB2[j,k,2])<epsilon*B2[j,k]):
342 dB2[j,k,2] = None
343 else:
344 dB2[j,k,2] = 0
345 if abs(B2[j,k])<epsilon:
346 dB2[j,k,2] = None
347
348 cdbz=(dip.get_old(xv[j],yv[j],zv[j]+step2*RE,0,k,0)+mdip.get_old(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.get_old(xv[j],yv[j],zv[j]-step2*RE,0,k,0) - mdip.get_old(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
349 if abs(cdbz) > epsilon*B3[j,k]:
350 dB3[j,k,2] = dB3[j,k,2]/cdbz
351 elif (abs(cdbz)<epsilon*B3[j,k]) and (abs(dB3[j,k,2])<epsilon*B3[j,k]):
352 dB3[j,k,2] = None
353 else:
354 dB3[j,k,2] = 0
355 if abs(B3[j,k])<epsilon:
356 dB3[j,k,2] = None
357
358 cdbz=(dip.get_ldp(xv[j],yv[j],zv[j]+step2*RE,0,k,0)+mdip.get_ldp(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.get_ldp(xv[j],yv[j],zv[j]-step2*RE,0,k,0) - mdip.get_ldp(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
359 if abs(cdbz) > epsilon*B4[j,k]:
360 dB4[j,k,2] = dB4[j,k,2]/cdbz
361 elif (abs(cdbz)<epsilon*B4[j,k]) and (abs(dB4[j,k,2])<epsilon*B4[j,k]):
362 dB4[j,k,2] = None
363 else:
364 dB4[j,k,2] = 0
365 if abs(B4[j,k])<epsilon:
366 dB4[j,k,2] = None
367
368 print(np.ma.amin(np.ma.masked_invalid(dB1)),np.ma.amax(np.ma.masked_invalid(dB1)))
369 print(np.ma.amin(np.ma.masked_invalid(dB2)),np.ma.amax(np.ma.masked_invalid(dB2)))
370 print(np.ma.amin(np.ma.masked_invalid(dB3)),np.ma.amax(np.ma.masked_invalid(dB3)))
371 print(np.ma.amin(np.ma.masked_invalid(dB4)),np.ma.amax(np.ma.masked_invalid(dB4)))
372
373 # print(np.amin(B1),np.amax(B1))
374 # print(np.amin(B2),np.amax(B2))
375 # print(np.amin(B3),np.amax(B3))
376 # print(np.amin(B4),np.amax(B4))
377
378 colors=['r','k','b']
379 coords = ['x','y','z']
380
381 # for k in range(3):
382 # B1[:,k] = abs(B1[:,k])
383 # B2[:,k] = abs(B2[:,k])
384 # B3[:,k] = abs(B3[:,k])
385 # B4[:,k] = abs(B4[:,k])
386 # for kk in range(3):
387 # dB1[:,k,kk] = abs(dB1[:,k,kk])
388 # dB2[:,k,kk] = abs(dB2[:,k,kk])
389 # dB3[:,k,kk] = abs(dB3[:,k,kk])
390 # dB4[:,k,kk] = abs(dB4[:,k,kk])
391
392
393 for k in range(3):
394 ax.plot(radiiRE, dB1[:,k,kkk], c=colors[k], linestyle='-', linewidth=linewidth, label='vectorpot dB'+coords[k]+'/d'+coords[kkk])
395 ax.plot(radiiRE, dB2[:,k,kkk], c=colors[k], linestyle='--', linewidth=linewidth, label='regular dB'+coords[k]+'/d'+coords[kkk])
396 ax.plot(radiiRE, dB3[:,k,kkk], c=colors[k], linestyle=':', linewidth=linewidth, label='reg+mirror dB'+coords[k]+'/d'+coords[kkk])
397 if tilt_angle_phi<epsilon:
398 ax.plot(radiiRE, dB4[:,k,kkk], c=colors[k], linestyle='-.', linewidth=linewidth, label='line+mirror dB'+coords[k]+'/d'+coords[kkk])
399
400 ax.text(0.2,0.08,r"profile with $\theta="+str(int(line_theta[i]*180./math.pi))+"$, $\phi="+str(int(line_phi[i]*180./math.pi))+"$",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
401
402
403 ax.set_xlabel(r"$r$ [$r_\mathrm{E}$]", fontsize=fontsize)
404 ax.set_xlim([1,70])
405 #ax.set_yscale('log', nonposy='clip')
406 #ax.set_ylim([1.e-6,1.e-1])
407 #ax.set_ylim([1.e-26,1.e-21])
408 ax.set_ylim([-.1,1.1])
409
410 for item in ax.get_xticklabels():
411 item.set_fontsize(fontsize)
412 for item in ax.get_yticklabels():
413 item.set_fontsize(fontsize)
414
415 handles, labels = axes[-1].get_legend_handles_labels()
416 axes[-1].legend(handles, labels, fontsize=fontsize).set_zorder(10)
417
418 fig.savefig(outfilename[:-4]+"_d"+coords[kkk]+outfilename[-4:])
419 plt.close()
420
421
abs
static ARCH_HOSTDEV VecSimple< T > abs(const VecSimple< T > &l)
Definition vectorclass_fallback.h:81