66import numpy as np
77import h5py
88from tqdm import tqdm
9- from scipy .constants import c as c_light
109from wakis .sources import Beam
11- from wakis .field_monitors import FieldMonitor
1210
1311class RoutinesMixin ():
1412
1513 def emsolve (self , Nt , source = None , callback = None ,
16- save = False , fields = ['E' ], components = ['Abs' ], save_every = 1 , subdomain = None ,
17- plot = False , plot_every = 1 , use_etd = False ,
14+ save = False , fields = ['E' ], components = ['Abs' ], save_every = 1 , subdomain = None ,
15+ plot = False , plot_every = 1 , use_etd = False ,
1816 plot3d = False , ** kwargs ):
1917 '''
2018 Run the simulation and save the selected field components in HDF5 files
@@ -26,13 +24,13 @@ def emsolve(self, Nt, source=None, callback=None,
2624 Nt: int
2725 Number of timesteps to run
2826 source: source object
29- source object from `sources.py` defining the time-dependednt source.
27+ source object from `sources.py` defining the time-dependednt source.
3028 It should have an update function `source.update(solver, t)`
3129 save: bool
3230 Flag to enable saving the field in HDF5 format
3331 fields: list, default ['E']
3432 3D field magnitude ('E', 'H', or 'J') to save
35- 'Ex', 'Hy', etc., is also accepted and will override
33+ 'Ex', 'Hy', etc., is also accepted and will override
3634 the `components` parameter.
3735 components: list, default ['z']
3836 Field compoonent ('x', 'y', 'z', 'Abs') to save. It will be overriden
@@ -49,10 +47,10 @@ def emsolve(self, Nt, source=None, callback=None,
4947 Number of timesteps between consecutive plots
5048 **kwargs:
5149 Keyword arguments to be passed to the Plot2D function.
52- * Default kwargs used for 2D plotting:
50+ * Default kwargs used for 2D plotting:
5351 {'field':'E', 'component':'z',
54- 'plane':'ZY', 'pos':0.5, 'title':'Ez',
55- 'cmap':'rainbow', 'patch_reverse':True,
52+ 'plane':'ZY', 'pos':0.5, 'title':'Ez',
53+ 'cmap':'rainbow', 'patch_reverse':True,
5654 'off_screen': True, 'interpolation':'spline36'}
5755 * Default kwargs used for 3D plotting:
5856 {'field':'E', 'component':'z',
@@ -71,8 +69,9 @@ def emsolve(self, Nt, source=None, callback=None,
7169 h5py
7270 '''
7371 self .Nt = Nt
74- if source is not None : self .source = source
75-
72+ if source is not None :
73+ self .source = source
74+
7675 if save :
7776
7877 hfs = {}
@@ -94,8 +93,8 @@ def emsolve(self, Nt, source=None, callback=None,
9493
9594 if plot :
9695 plotkw = {'field' :'E' , 'component' :'z' ,
97- 'plane' :'ZY' , 'pos' :0.5 , 'cmap' :'rainbow' ,
98- 'patch_reverse' :True , 'title' :'Ez' ,
96+ 'plane' :'ZY' , 'pos' :0.5 , 'cmap' :'rainbow' ,
97+ 'patch_reverse' :True , 'title' :'Ez' ,
9998 'off_screen' : True , 'interpolation' :'spline36' }
10099 plotkw .update (kwargs )
101100
@@ -105,28 +104,28 @@ def emsolve(self, Nt, source=None, callback=None,
105104 'title' :'Ez' , 'cmap' :'jet' , 'clip_volume' :False , 'clip_normal' :'-y' ,
106105 'field_on_stl' :True , 'field_opacity' :1.0 ,
107106 'off_screen' :True , 'zoom' :1.0 , 'nan_opacity' :1.0 }
108-
107+
109108 plotkw .update (kwargs )
110109
111110 # get ABC values
112111 if self .activate_abc :
113112 E_abc_2 , H_abc_2 = self .get_abc ()
114113 E_abc_1 , H_abc_1 = self .get_abc ()
115114
116- # Time loop
115+ # Time loop
117116 for n in tqdm (range (Nt )):
118117
119- if source is not None :
118+ if source is not None :
120119 source .update (self , n * self .dt )
121120
122121 if save and n % save_every == 0 :
123122 for field in hfs .keys ():
124123 try :
125124 d = getattr (self , field [0 ])[xx ,yy ,zz ,field [1 :]]
126- except :
125+ except AttributeError :
127126 raise (f'Component { field } not valid. Input must have a \
128127 field ["E", "H", "J"] and a component ["x", "y", "z", "Abs"]' )
129-
128+
130129 # Save timestep in HDF5
131130 hfs [field ]['#' + str (n ).zfill (5 )] = d
132131
@@ -155,38 +154,38 @@ def emsolve(self, Nt, source=None, callback=None,
155154 for hf in hfs :
156155 hf .close ()
157156
158- def wakesolve (self , wakelength ,
159- wake = None ,
157+ def wakesolve (self , wakelength ,
158+ wake = None ,
160159 callback = None ,
161- compute_plane = 'both' ,
160+ compute_plane = 'both' ,
162161 plot = False , plot_from = None , plot_every = 1 , plot_until = None ,
163- save_J = False ,
162+ save_J = False ,
164163 add_space = None , #for legacy
165164 use_field_monitor = False ,
166165 field_monitor = None ,
167- use_edt = None , #deprecated
166+ use_edt = None , #deprecated
168167 ** kwargs ):
169168 '''
170169 Run the EM simulation and compute the longitudinal (z) and transverse (x,y)
171- wake potential WP(s) and impedance Z(s).
172-
173- The `Ez` field is saved every timestep in a subdomain (xtest, ytest, z) around
170+ wake potential WP(s) and impedance Z(s).
171+
172+ The `Ez` field is saved every timestep in a subdomain (xtest, ytest, z) around
174173 the beam trajectory in HDF5 format file `Ez.h5`.
175174
176- The computed results are available as Solver class attributes:
175+ The computed results are available as Solver class attributes:
177176 - wake potential: WP (longitudinal), WPx, WPy (transverse) [V/pC]
178177 - impedance: Z (longitudinal), Zx, Zy (transverse) [Ohm]
179178 - beam charge distribution: lambdas (distance) [C/m] lambdaf (spectrum) [C]
180179
181180 Parameters:
182181 -----------
183182 wakelength: float
184- Desired length of the wake in [m] to be computed
185-
183+ Desired length of the wake in [m] to be computed
184+
186185 Maximum simulation time in [s] can be computed from the wakelength parameter as:
187- .. math:: t_{max} = t_{inj} + (wakelength + (z_{max}-z_{min}))/c
186+ .. math:: t_{max} = t_{inj} + (wakelength + (z_{max}-z_{min}))/c
188187 wake: Wake obj, default None
189- `Wake()` object containing the information needed to run
188+ `Wake()` object containing the information needed to run
190189 the wake solver calculation. See Wake() docstring for more information.
191190 Can be passed at `Solver()` instantiation as parameter too.
192191 save_J: bool, default False
@@ -201,11 +200,11 @@ def wakesolve(self, wakelength,
201200 Number of timesteps between consecutive plots
202201 **kwargs:
203202 Keyword arguments to be passed to the Plot2D function.
204- Default kwargs used:
205- {'plane':'ZY', 'pos':0.5, 'title':'Ez',
206- 'cmap':'rainbow', 'patch_reverse':True,
203+ Default kwargs used:
204+ {'plane':'ZY', 'pos':0.5, 'title':'Ez',
205+ 'cmap':'rainbow', 'patch_reverse':True,
207206 'off_screen': True, 'interpolation':'spline36'}
208-
207+
209208 Raises:
210209 -------
211210 AttributeError:
@@ -218,7 +217,8 @@ def wakesolve(self, wakelength,
218217 h5py
219218 '''
220219
221- if wake is not None : self .wake = wake
220+ if wake is not None :
221+ self .wake = wake
222222 if self .wake is None :
223223 raise ('Wake solver information not passed to the solver instantiation' )
224224
@@ -228,26 +228,26 @@ def wakesolve(self, wakelength,
228228 # plot params defaults
229229 if plot :
230230 plotkw = {'plane' :'ZY' , 'pos' :0.5 , 'title' :'Ez' ,
231- 'cmap' :'rainbow' , 'patch_reverse' :True ,
231+ 'cmap' :'rainbow' , 'patch_reverse' :True ,
232232 'off_screen' : True , 'interpolation' :'spline36' }
233233 plotkw .update (kwargs )
234-
234+
235235 # integration path (test position)
236236 self .xtest , self .ytest = self .wake .xtest , self .wake .ytest
237237 self .ixt , self .iyt = np .abs (self .x - self .xtest ).argmin (), np .abs (self .y - self .ytest ).argmin ()
238238 if compute_plane .lower () == 'longitudinal' :
239239 xx , yy = self .ixt , self .iyt
240240 else :
241- xx , yy = slice (self .ixt - 1 , self .ixt + 2 ), slice (self .iyt - 1 , self .iyt + 2 )
241+ xx , yy = slice (self .ixt - 1 , self .ixt + 2 ), slice (self .iyt - 1 , self .iyt + 2 )
242242
243- # Compute simulation time
243+ # Compute simulation time
244244 self .wake .wakelength = wakelength
245245 self .ti = self .wake .ti
246246 self .v = self .wake .v
247247 if self .use_mpi : #E- should it be zmin, zmax instead?
248248 z = self .Z # use global coords
249249 zz = slice (0 , self .NZ )
250- else :
250+ else :
251251 z = self .z
252252 zz = slice (0 , self .Nz )
253253
@@ -256,13 +256,13 @@ def wakesolve(self, wakelength,
256256 self .tmax , self .Nt = tmax , Nt
257257
258258 # Add beam source
259- beam = Beam (q = self .wake .q ,
260- sigmaz = self .wake .sigmaz ,
259+ beam = Beam (q = self .wake .q ,
260+ sigmaz = self .wake .sigmaz ,
261261 beta = self .wake .beta ,
262- xsource = self .wake .xsource , ysource = self .wake .ysource ,
262+ xsource = self .wake .xsource , ysource = self .wake .ysource ,
263263 )
264-
265- # hdf5
264+
265+ # hdf5
266266 self .Ez_file = self .wake .Ez_file
267267 hf = None # needed for MPI
268268 if self .use_mpi :
@@ -291,10 +291,12 @@ def save_to_h5(self, hf, field, x, y, z, comp, n):
291291 if self .rank == 0 :
292292 hf ['#' + str (n ).zfill (5 )] = _field
293293 else :
294- hf ['#' + str (n ).zfill (5 )] = getattr (self , field )[x , y , z , comp ]
295-
296- if plot_until is None : plot_until = Nt
297- if plot_from is None : plot_from = int (self .ti / self .dt )
294+ hf ['#' + str (n ).zfill (5 )] = getattr (self , field )[x , y , z , comp ]
295+
296+ if plot_until is None :
297+ plot_until = Nt
298+ if plot_from is None :
299+ plot_from = int (self .ti / self .dt )
298300
299301 print ('Running electromagnetic time-domain simulation...' )
300302 for n in tqdm (range (Nt )):
@@ -303,10 +305,10 @@ def save_to_h5(self, hf, field, x, y, z, comp, n):
303305 beam .update (self , n * self .dt )
304306
305307 # Save
306- save_to_h5 (self , hf , 'E' , xx , yy , zz , 'z' , n )
308+ save_to_h5 (self , hf , 'E' , xx , yy , zz , 'z' , n )
307309 if save_J :
308- save_to_h5 (self , hfJ , 'J' , xx , yy , zz , 'z' , n )
309-
310+ save_to_h5 (self , hfJ , 'J' , xx , yy , zz , 'z' , n )
311+
310312 # Advance
311313 self .one_step ()
312314
@@ -324,7 +326,7 @@ def save_to_h5(self, hf, field, x, y, z, comp, n):
324326 if callback is not None :
325327 callback (self , n * self .dt )
326328
327- # End of time loop
329+ # End of time loop
328330 if self .use_mpi :
329331 if self .rank == 0 :
330332 hf .close ()
@@ -340,4 +342,4 @@ def save_to_h5(self, hf, field, x, y, z, comp, n):
340342
341343 # Compute wakefield magnitudes is done inside WakeSolver
342344 self .wake .solve (compute_plane = compute_plane )
343-
345+
0 commit comments