29import matplotlib.pyplot
as plt
32''' Testing routine for different dipole formulations
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.
37 If arg2 is present, the code also calculates verification of derivative terms.
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.
46 testset = int(sys.argv[1])
57plt.switch_backend(
'Agg')
59outfilename =
"./vecpotdip_verify2_"+str(testset)+
".png"
70 line_start = np.array([0,0,0])
75 line_start = np.array([0,0,0])
80 line_start = np.array([-3,-3,-3])
85 line_start = np.array([-3,-3,-3])
90 line_start = np.array([0,0,0])
93 tilt_angle_theta = 45.
95 line_start = np.array([0,0,0])
100 line_start = np.array([0,0,0])
105 tilt_angle_theta = 0.
107 line_start = np.array([0,0,0])
109print(
"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))
111line_phi = np.array([0,45,80,90,110,135])*math.pi/180.
112line_theta = np.zeros(len(line_phi)) + line_theta * math.pi/180.
127fig.set_size_inches(20,30)
128nsubplots=len(line_theta)
129for i
in range(nsubplots):
130 fig.add_subplot(nsubplots,1,i+1)
133radii = np.arange(0.1,100,step)*RE
137fig.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)
139for i
in range(nsubplots):
143 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)
145 xv = line_start[0]*RE + radii*np.sin(line_phi[i])*np.cos(line_theta[i])
146 yv = line_start[1]*RE + radii*np.sin(line_phi[i])*np.sin(line_theta[i])
147 zv = line_start[2]*RE + radii*np.cos(line_phi[i])
149 B1 = np.zeros([nr,4])
150 B2 = np.zeros([nr,4])
151 B3 = np.zeros([nr,4])
152 B4 = np.zeros([nr,4])
156 B1[j,k] = dip.getX(xv[j],yv[j],zv[j],0,k,0) + imfpot.get(xv[j],yv[j],zv[j],0,k,0)
157 B2[j,k] = dip.get_old(xv[j],yv[j],zv[j],0,k,0)
158 B3[j,k] = B2[j,k] + mdip.get_old(xv[j],yv[j],zv[j],0,k,0)
159 B4[j,k] = dip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
160 B4[j,k] = B4[j,k] + mdip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
161 if plotmagnitude
is True:
162 B1[j,3] = np.linalg.norm(B1[j,0:3])
163 B2[j,3] = np.linalg.norm(B2[j,0:3])
164 B3[j,3] = np.linalg.norm(B3[j,0:3])
165 B4[j,3] = np.linalg.norm(B4[j,0:3])
167 colors=[
'r',
'k',
'b',
'magenta']
168 coords = [
'x',
'y',
'z',
'mag']
171 if plotmagnitude
is True:
174 ax.plot(radiiRE, B1[:,k], c=colors[k], linestyle=
'-', linewidth=linewidth, label=
'vectorpot B'+coords[k], zorder=-10)
175 ax.plot(radiiRE, B2[:,k], c=colors[k], linestyle=
'--', linewidth=linewidth, label=
'regular B'+coords[k])
176 ax.plot(radiiRE, B3[:,k], c=colors[k], linestyle=
':', linewidth=linewidth, label=
'reg+mirror B'+coords[k])
177 if tilt_angle_phi<epsilon:
178 ax.plot(radiiRE, B4[:,k], c=colors[k], linestyle=
'-.', linewidth=linewidth, label=
'line+mirror B'+coords[k])
180 ax.set_xlabel(
r"$r$ [$r_\mathrm{E}$]", fontsize=fontsize)
183 ax.set_yscale(
'symlog', linthreshy=linthresh)
184 for item
in ax.get_xticklabels():
185 item.set_fontsize(fontsize)
186 for item
in ax.get_yticklabels():
187 item.set_fontsize(fontsize)
189 ylims = np.array(ax.get_ylim())
196handles, labels = axes[-1].get_legend_handles_labels()
197axes[-1].legend(handles, labels, fontsize=fontsize)
199fig.savefig(outfilename)
208 print(
"derivatives d"+coords[kkk])
209 dB1 = np.zeros([nr,3,3])
210 dB2 = np.zeros([nr,3,3])
211 dB3 = np.zeros([nr,3,3])
212 dB4 = np.zeros([nr,3,3])
216 fig.set_size_inches(20,30)
217 for i
in range(nsubplots):
218 fig.add_subplot(nsubplots,1,i+1)
219 axes = fig.get_axes()
221 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)
223 for i
in range(nsubplots):
224 print(
"derivatives subplot ",i)
227 xv = line_start[0]*RE + radii*np.sin(line_phi[i])*np.cos(line_theta[i])
228 yv = line_start[1]*RE + radii*np.sin(line_phi[i])*np.sin(line_theta[i])
229 zv = line_start[2]*RE + radii*np.cos(line_phi[i])
233 B1[j,k] = dip.getX(xv[j],yv[j],zv[j],0,k,0) + imfpot.get(xv[j],yv[j],zv[j],0,k,0)
234 B2[j,k] = dip.get_old(xv[j],yv[j],zv[j],0,k,0)
235 B3[j,k] = B2[j,k] + mdip.get_old(xv[j],yv[j],zv[j],0,k,0)
236 B4[j,k] = dip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
237 B4[j,k] = B4[j,k] + mdip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
240 dB1[j,k,kk] = dip.getX(xv[j],yv[j],zv[j],1,k,kk) + imfpot.get(xv[j],yv[j],zv[j],1,k,kk)
241 dB2[j,k,kk] = dip.get_old(xv[j],yv[j],zv[j],1,k,kk)
242 dB3[j,k,kk] = dB2[j,k,kk] + mdip.get_old(xv[j],yv[j],zv[j],1,k,kk)
243 dB4[j,k,kk] = dip.get_ldp(xv[j],yv[j],zv[j],1,k,kk)
244 dB4[j,k,kk] = dB4[j,k,kk] + mdip.get_ldp(xv[j],yv[j],zv[j],1,k,kk)
247 for j
in np.arange(1,nr-1):
253 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) + imfpot.get(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - imfpot.get(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
254 if abs(cdbx) > epsilon*B1[j,k]:
255 dB1[j,k,0] = dB1[j,k,0]/cdbx
256 elif (
abs(cdbx)<epsilon*B1[j,k])
and (
abs(dB1[j,k,0])<epsilon*B1[j,k]):
260 if abs(B1[j,k])<epsilon:
263 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)
264 if abs(cdbx) > epsilon*B2[j,k]:
265 dB2[j,k,0] = dB2[j,k,0]/cdbx
266 elif (
abs(cdbx)<epsilon*B2[j,k])
and (
abs(dB2[j,k,0])<epsilon*B2[j,k]):
270 if abs(B2[j,k])<epsilon:
273 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)
274 if abs(cdbx) > epsilon*B3[j,k]:
275 dB3[j,k,0] = dB3[j,k,0]/cdbx
276 elif (
abs(cdbx)<epsilon*B3[j,k])
and (
abs(dB3[j,k,0])<epsilon*B3[j,k]):
280 if abs(B3[j,k])<epsilon:
283 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)
284 if abs(cdbx) > epsilon*B4[j,k]:
285 dB4[j,k,0] = dB4[j,k,0]/cdbx
286 elif (
abs(cdbx)<epsilon*B4[j,k])
and (
abs(dB4[j,k,0])<epsilon*B4[j,k]):
290 if abs(B4[j,k])<epsilon:
295 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)
296 if abs(cdby) > epsilon*B1[j,k]:
297 dB1[j,k,1] = dB1[j,k,1]/cdby
298 elif (
abs(cdby)<epsilon*B1[j,k])
and (
abs(dB1[j,k,1])<epsilon*B1[j,k]):
302 if abs(B1[j,k])<epsilon:
305 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)
306 if abs(cdby) > epsilon*B2[j,k]:
307 dB2[j,k,1] = dB2[j,k,1]/cdby
308 elif (
abs(cdby)<epsilon*B2[j,k])
and (
abs(dB2[j,k,1])<epsilon*B2[j,k]):
312 if abs(B2[j,k])<epsilon:
315 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)
316 if abs(cdby) > epsilon*B3[j,k]:
317 dB3[j,k,1] = dB3[j,k,1]/cdby
318 elif (
abs(cdby)<epsilon*B3[j,k])
and (
abs(dB3[j,k,1])<epsilon*B3[j,k]):
322 if abs(B3[j,k])<epsilon:
325 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)
326 if abs(cdby) > epsilon*B4[j,k]:
327 dB4[j,k,1] = dB4[j,k,1]/cdby
328 elif (
abs(cdby)<epsilon*B4[j,k])
and (
abs(dB4[j,k,1])<epsilon*B4[j,k]):
332 if abs(B4[j,k])<epsilon:
337 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)
338 if abs(cdbz) > epsilon*B1[j,k]:
339 dB1[j,k,2] = dB1[j,k,2]/cdbz
340 elif (
abs(cdbz)<epsilon*B1[j,k])
and (
abs(dB1[j,k,2])<epsilon*B1[j,k]):
344 if abs(B1[j,k])<epsilon:
347 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)
348 if abs(cdbz) > epsilon*B2[j,k]:
349 dB2[j,k,2] = dB2[j,k,2]/cdbz
350 elif (
abs(cdbz)<epsilon*B2[j,k])
and (
abs(dB2[j,k,2])<epsilon*B2[j,k]):
354 if abs(B2[j,k])<epsilon:
357 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)
358 if abs(cdbz) > epsilon*B3[j,k]:
359 dB3[j,k,2] = dB3[j,k,2]/cdbz
360 elif (
abs(cdbz)<epsilon*B3[j,k])
and (
abs(dB3[j,k,2])<epsilon*B3[j,k]):
364 if abs(B3[j,k])<epsilon:
367 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)
368 if abs(cdbz) > epsilon*B4[j,k]:
369 dB4[j,k,2] = dB4[j,k,2]/cdbz
370 elif (
abs(cdbz)<epsilon*B4[j,k])
and (
abs(dB4[j,k,2])<epsilon*B4[j,k]):
374 if abs(B4[j,k])<epsilon:
377 print(np.ma.amin(np.ma.masked_invalid(dB1)),np.ma.amax(np.ma.masked_invalid(dB1)))
378 print(np.ma.amin(np.ma.masked_invalid(dB2)),np.ma.amax(np.ma.masked_invalid(dB2)))
379 print(np.ma.amin(np.ma.masked_invalid(dB3)),np.ma.amax(np.ma.masked_invalid(dB3)))
380 print(np.ma.amin(np.ma.masked_invalid(dB4)),np.ma.amax(np.ma.masked_invalid(dB4)))
388 coords = [
'x',
'y',
'z']
403 ax.plot(radiiRE, dB1[:,k,kkk], c=colors[k], linestyle=
'-', linewidth=linewidth, label=
'vectorpot dB'+coords[k]+
'/d'+coords[kkk])
404 ax.plot(radiiRE, dB2[:,k,kkk], c=colors[k], linestyle=
'--', linewidth=linewidth, label=
'regular dB'+coords[k]+
'/d'+coords[kkk])
405 ax.plot(radiiRE, dB3[:,k,kkk], c=colors[k], linestyle=
':', linewidth=linewidth, label=
'reg+mirror dB'+coords[k]+
'/d'+coords[kkk])
406 if tilt_angle_phi<epsilon:
407 ax.plot(radiiRE, dB4[:,k,kkk], c=colors[k], linestyle=
'-.', linewidth=linewidth, label=
'line+mirror dB'+coords[k]+
'/d'+coords[kkk])
409 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)
412 ax.set_xlabel(
r"$r$ [$r_\mathrm{E}$]", fontsize=fontsize)
417 ax.set_ylim([-.1,1.1])
419 for item
in ax.get_xticklabels():
420 item.set_fontsize(fontsize)
421 for item
in ax.get_yticklabels():
422 item.set_fontsize(fontsize)
424 handles, labels = axes[-1].get_legend_handles_labels()
425 axes[-1].legend(handles, labels, fontsize=fontsize).set_zorder(10)
427 fig.savefig(outfilename[:-4]+
"_d"+coords[kkk]+outfilename[-4:])
static ARCH_HOSTDEV VecSimple< T > abs(const VecSimple< T > &l)