Coverage for kwave/kWaveSimulation.py: 18%
581 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-03-24 12:06 -0700
1import logging
2from dataclasses import dataclass
4import numpy as np
5from deprecated import deprecated
7from kwave.data import Vector
8from kwave.kgrid import kWaveGrid
9from kwave.kmedium import kWaveMedium
10from kwave.ksensor import kSensor
11from kwave.ksource import kSource
12from kwave.ktransducer import NotATransducer
13from kwave.kWaveSimulation_helper import (
14 create_absorption_variables,
15 display_simulation_params,
16 expand_grid_matrices,
17 scale_source_terms_func,
18 set_sound_speed_ref,
19)
20from kwave.options.simulation_options import SimulationOptions, SimulationType
21from kwave.recorder import Recorder
22from kwave.utils.checks import check_stability
23from kwave.utils.colormap import get_color_map
24from kwave.utils.conversion import cart2grid, cast_to_type
25from kwave.utils.data import get_date_string, get_smallest_possible_type
26from kwave.utils.dotdictionary import dotdict
27from kwave.utils.filters import smooth
28from kwave.utils.matlab import matlab_find, matlab_mask
29from kwave.utils.matrix import num_dim2
32@dataclass
33class kWaveSimulation(object):
34 def __init__(
35 self, kgrid: kWaveGrid, source: kSource, sensor: NotATransducer, medium: kWaveMedium, simulation_options: SimulationOptions
36 ):
37 self.precision = None
38 self.kgrid = kgrid
39 self.medium = medium
40 self.source = source
41 self.sensor = sensor
42 self.options = simulation_options
44 # =========================================================================
45 # FLAGS WHICH DEPEND ON USER INPUTS (THESE SHOULD NOT BE MODIFIED)
46 # =========================================================================
47 # flags which control the characteristics of the sensor
48 #: Whether the sensor sensor mask is defined by cuboid corners
49 #: Whether time reversal simulation is enabled
51 # check if the sensor mask is defined as a list of cuboid corners
52 if self.sensor.mask is not None and self.sensor.mask.shape[0] == (2 * self.kgrid.dim): 52 ↛ 53line 52 didn't jump to line 53 because the condition on line 52 was never true
53 self.userarg_cuboid_corners = True
54 else:
55 self.userarg_cuboid_corners = False
57 # check if performing time reversal, and replace inputs to explicitly use a
58 # source with a dirichlet boundary condition
59 if self.sensor.time_reversal_boundary_data is not None: 59 ↛ 77line 59 didn't jump to line 77 because the condition on line 59 was always true
60 # define a new source structure
61 self.source = kSource()
62 self.source.p_mask = self.sensor.mask
63 self.source.p = np.flip(self.sensor.time_reversal_boundary_data, 1)
64 self.source.p_mode = "dirichlet"
66 # define a new sensor structure
67 Nx, Ny, Nz = self.kgrid.Nx, self.kgrid.Ny, self.kgrid.Nz
68 self.sensor = kSensor(mask=np.ones((Nx, Ny, max(1, Nz))), record=["p_final"])
70 # set time reversal flag
71 self.userarg_time_rev = True
73 self.record = Recorder()
74 self.record.p = False
75 else:
76 # set time reversal flag
77 self.userarg_time_rev = False
79 #: Whether sensor.mask should be re-ordered.
80 #: True if sensor.mask is Cartesian with nearest neighbour interpolation which is calculated using a binary mask
81 #: and thus must be re-ordered
82 self.reorder_data = False
84 #: Whether the sensor.mask is binary
85 self.binary_sensor_mask = True
87 #: If tse sensor is an object of the kWaveTransducer class
88 self.transducer_sensor = False
90 self.record = Recorder()
92 # transducer source flags
93 #: transducer is object of kWaveTransducer class
94 self.transducer_source = False
96 #: Apply receive elevation focus on the transducer
97 self.transducer_receive_elevation_focus = False
99 # general
100 self.COLOR_MAP = get_color_map() #: default color map
101 self.ESTIMATE_SIM_TIME_STEPS = 50 #: time steps used to estimate simulation time
102 self.HIGHEST_PRIME_FACTOR_WARNING = 7 #: largest prime factor before warning
103 self.KSPACE_CFL = 0.3 #: default CFL value used if kgrid.t_array is set to 'auto'
104 self.PSTD_CFL = 0.1 #: default CFL value used if kgrid.t_array is set to 'auto'
106 # source types
107 self.SOURCE_S_MODE_DEF = "additive" #: source mode for stress sources
108 self.SOURCE_P_MODE_DEF = "additive" #: source mode for pressure sources
109 self.SOURCE_U_MODE_DEF = "additive" #: source mode for velocity sources
111 # filenames
112 self.STREAM_TO_DISK_FILENAME = "temp_sensor_data.bin" #: default disk stream filename
113 self.LOG_NAME = ["k-Wave-Log-", get_date_string()] #: default log filename
115 self.calling_func_name = None
116 logging.log(logging.INFO, f" start time: {get_date_string()}")
118 self.c_ref, self.c_ref_compression, self.c_ref_shear = [None] * 3
119 self.transducer_input_signal = None
121 #: Indexing variable corresponding to the location of all the pressure source elements
122 self.p_source_pos_index = None
123 #: Indexing variable corresponding to the location of all the velocity source elements
124 self.u_source_pos_index = None
125 #: Indexing variable corresponding to the location of all the stress source elements
126 self.s_source_pos_index = None
128 #: Delay mask that accounts for the beamforming delays and elevation focussing
129 self.delay_mask = None
131 self.absorb_nabla1 = None #: absorbing fractional Laplacian operator
132 self.absorb_tau = None #: absorbing fractional Laplacian coefficient
133 self.absorb_nabla2 = None #: dispersive fractional Laplacian operator
134 self.absorb_eta = None #: dispersive fractional Laplacian coefficient
136 self.dt = None #: Alias to kgrid.dt
137 self.rho0 = None #: Alias to medium.density
138 self.c0 = None #: Alias to medium.sound_speed
139 self.index_data_type = None
141 @property
142 def equation_of_state(self):
143 """
144 Returns:
145 Set equation of state variable
146 """
147 if self.medium.absorbing:
148 if self.medium.stokes:
149 return "stokes"
150 else:
151 return "absorbing"
152 else:
153 return "loseless"
155 @property
156 def use_sensor(self):
157 """
158 Returns:
159 False if no output of any kind is required
161 """
162 return self.sensor is not None
164 @property
165 def blank_sensor(self):
166 """
167 Returns
168 True if sensor.mask is not defined but _max_all or _final variables are still recorded
170 """
171 fields = ["p", "p_max", "p_min", "p_rms", "u", "u_non_staggered", "u_split_field", "u_max", "u_min", "u_rms", "I", "I_avg"]
172 if not (isinstance(self.sensor, NotATransducer) or any(self.record.is_set(fields)) or self.time_rev):
173 return True
174 return False
176 @property
177 def kelvin_voigt_model(self):
178 """
179 Returns:
180 Whether the simulation is elastic with absorption
182 """
183 return False
185 @property
186 def nonuniform_grid(self):
187 """
188 Returns:
189 True if the computational grid is non-uniform
191 """
192 return self.kgrid.nonuniform
194 @property
195 @deprecated(version="0.4.1", reason="Use TimeReversal class instead")
196 def time_rev(self) -> bool:
197 if self.sensor and not isinstance(self.sensor, NotATransducer): 197 ↛ 199line 197 didn't jump to line 199 because the condition on line 197 was always true
198 return not self.options.simulation_type.is_elastic_simulation() and self.sensor.time_reversal_boundary_data is not None
199 return self.userarg_time_rev
201 @property
202 @deprecated(version="0.4.1", reason="Use TimeReversal class instead")
203 def elastic_time_rev(self):
204 """
205 Returns:
206 True if using time reversal with the elastic code
208 """
209 return False
211 @property
212 def compute_directivity(self):
213 """
214 Returns:
215 True if directivity calculations in 2D are used by setting sensor.directivity_angle
217 """
218 if self.sensor is not None and not isinstance(self.sensor, NotATransducer):
219 if self.kgrid.dim == 2:
220 # check for sensor directivity input and set flag
221 directivity = self.sensor.directivity
222 if directivity is not None and directivity.angle is not None:
223 return True
224 return False
226 @property
227 def cuboid_corners(self):
228 """
229 Returns:
230 Whether the sensor.mask is a list of cuboid corners
231 """
232 if self.sensor is not None and not isinstance(self.sensor, NotATransducer):
233 if not self.blank_sensor and self.sensor.mask.shape[0] == 2 * self.kgrid.dim:
234 return True
235 return self.userarg_cuboid_corners
237 ##############
238 # flags which control the types of source used
239 ##############
240 @property
241 def source_p0(self): # initial pressure
242 """
243 Returns:
244 Whether initial pressure source is present (default=False)
246 """
247 flag = False # default
248 if not isinstance(self.source, NotATransducer) and self.source.p0 is not None:
249 # set flag
250 flag = True
251 return flag
253 @property
254 def source_p0_elastic(self): # initial pressure in the elastic code
255 """
256 Returns:
257 Whether initial pressure source is present in the elastic code (default=False)
259 """
260 # Not clear where this flag is set
261 return False
263 @property
264 def source_p(self):
265 """
266 Returns:
267 Whether time-varying pressure source is present (default=False)
269 """
270 flag = False
271 if not isinstance(self.source, NotATransducer) and self.source.p is not None:
272 # set source flag to the length of the source, this allows source.p
273 # to be shorter than kgrid.Nt
274 flag = len(self.source.p[0])
275 return flag
277 @property
278 def source_p_labelled(self): # time-varying pressure with labelled source mask
279 """
280 Returns:
281 True/False if labelled/binary source mask, respectively.
283 """
284 flag = False
285 if not isinstance(self.source, NotATransducer) and self.source.p is not None:
286 # check if the mask is binary or labelled
287 p_unique = np.unique(self.source.p_mask)
288 flag = not (p_unique.size <= 2 and p_unique.sum() == 1)
289 return flag
291 @property
292 def source_ux(self) -> bool:
293 """
294 Returns:
295 Whether time-varying particle velocity source is used in X-direction
297 """
298 flag = False
299 if not isinstance(self.source, NotATransducer) and self.source.ux is not None:
300 # set source flgs to the length of the sources, this allows the
301 # inputs to be defined independently and be of any length
302 flag = len(self.source.ux[0])
303 return flag
305 @property
306 def source_uy(self) -> bool:
307 """
308 Returns:
309 Whether time-varying particle velocity source is used in Y-direction
311 """
312 flag = False
313 if not isinstance(self.source, NotATransducer) and self.source.uy is not None:
314 # set source flgs to the length of the sources, this allows the
315 # inputs to be defined independently and be of any length
316 flag = len(self.source.uy[0])
317 return flag
319 @property
320 def source_uz(self) -> bool:
321 """
322 Returns:
323 Whether time-varying particle velocity source is used in Z-direction
325 """
326 flag = False
327 if not isinstance(self.source, NotATransducer) and self.source.uz is not None:
328 # set source flgs to the length of the sources, this allows the
329 # inputs to be defined independently and be of any length
330 flag = len(self.source.uz[0])
331 return flag
333 @property
334 def source_u_labelled(self):
335 """
336 Returns:
337 Whether time-varying velocity source with labelled source mask is present (default=False)
339 """
340 flag = False
341 if not isinstance(self.source, NotATransducer) and self.source.u_mask is not None:
342 # check if the mask is binary or labelled
343 u_unique = np.unique(self.source.u_mask)
344 if u_unique.size <= 2 and u_unique.sum() == 1:
345 # binary source mask
346 flag = False
347 else:
348 # labelled source mask
349 flag = True
350 return flag
352 @property
353 def source_sxx(self):
354 """
355 Returns:
356 Whether time-varying stress source in X->X direction is present (default=False)
357 """
358 flag = False
359 if not isinstance(self.source, NotATransducer) and self.source.sxx is not None:
360 flag = len(self.source.sxx[0])
361 return flag
363 @property
364 def source_syy(self):
365 """
366 Returns:
367 Whether time-varying stress source in Y->Y direction is present (default=False)
369 """
370 flag = False
371 if not isinstance(self.source, NotATransducer) and self.source.syy is not None:
372 flag = len(self.source.syy[0])
373 return flag
375 @property
376 def source_szz(self):
377 """
378 Returns:
379 Whether time-varying stress source in Z->Z direction is present (default=False)
381 """
382 flag = False
383 if not isinstance(self.source, NotATransducer) and self.source.szz is not None:
384 flag = len(self.source.szz[0])
385 return flag
387 @property
388 def source_sxy(self):
389 """
390 Returns:
391 Whether time-varying stress source in X->Y direction is present (default=False)
393 """
394 flag = False
395 if not isinstance(self.source, NotATransducer) and self.source.sxy is not None:
396 flag = len(self.source.sxy[0])
397 return flag
399 @property
400 def source_sxz(self):
401 """
402 Returns:
403 Whether time-varying stress source in X->Z direction is present (default=False)
405 """
406 flag = False
407 if not isinstance(self.source, NotATransducer) and self.source.sxz is not None:
408 flag = len(self.source.sxz[0])
409 return flag
411 @property
412 def source_syz(self):
413 """
414 Returns:
415 Whether time-varying stress source in Y->Z direction is present (default=False)
417 """
418 flag = False
419 if not isinstance(self.source, NotATransducer) and self.source.syz is not None:
420 flag = len(self.source.syz[0])
421 return flag
423 @property
424 def source_s_labelled(self):
425 """
426 Returns:
427 Whether time-varying stress source with labelled source mask is present (default=False)
429 """
430 flag = False
431 if not isinstance(self.source, NotATransducer) and self.source.s_mask is not None:
432 # check if the mask is binary or labelled
433 s_unique = np.unique(self.source.s_mask)
434 if s_unique.size <= 2 and s_unique.sum() == 1:
435 # binary source mask
436 flag = False
437 else:
438 # labelled source mask
439 flag = True
440 return flag
442 @property
443 def use_w_source_correction_p(self):
444 """
445 Returns:
446 Whether to use the w source correction instead of the k-space source correction for pressure sources
447 """
448 flag = False
449 if not isinstance(self.source, NotATransducer):
450 if self.source.p is not None and self.source.p_frequency_ref is not None:
451 flag = True
452 return flag
454 @property
455 def use_w_source_correction_u(self):
456 """
457 Returns:
458 Whether to use the w source correction instead of the k-space source correction for velocity sources
459 """
460 flag = False
461 if not isinstance(self.source, NotATransducer):
462 if any([(getattr(self.source, k) is not None) for k in ["ux", "uy", "uz", "u_mask"]]):
463 if self.source.u_frequency_ref is not None:
464 flag = True
465 return flag
467 def input_checking(self, calling_func_name) -> None:
468 """
469 Check the input fields for correctness and validness
471 Args:
472 calling_func_name: Name of the script that calls this function
474 Returns:
475 None
476 """
477 self.calling_func_name = calling_func_name
479 k_dim = self.kgrid.dim
481 self.check_calling_func_name_and_dim(calling_func_name, k_dim)
483 # run subscript to check optional inputs
484 self.options = SimulationOptions.option_factory(self.kgrid, self.options)
485 opt = self.options
487 # TODO(Walter): clean this up with getters in simulation options pml size
488 pml_x_size, pml_y_size, pml_z_size = opt.pml_x_size, opt.pml_y_size, opt.pml_z_size
489 pml_size = Vector([pml_x_size, pml_y_size, pml_z_size])
491 is_elastic_code = opt.simulation_type.is_elastic_simulation()
492 self.print_start_status(is_elastic_code=is_elastic_code)
493 self.set_index_data_type()
495 user_medium_density_input = self.check_medium(self.medium, self.kgrid.k, simulation_type=opt.simulation_type)
497 # select the reference sound speed used in the k-space operator
498 self.c_ref, self.c_ref_compression, self.c_ref_shear = set_sound_speed_ref(self.medium, opt.simulation_type)
500 self.check_source(k_dim, self.kgrid.Nt)
501 self.check_sensor(k_dim)
502 self.check_kgrid_time()
503 self.precision = self.select_precision(opt)
504 self.check_input_combinations(opt, user_medium_density_input, k_dim, pml_size, self.kgrid.N)
506 # run subscript to display time step, max supported frequency etc.
507 display_simulation_params(self.kgrid, self.medium, is_elastic_code)
509 self.smooth_and_enlarge(self.source, k_dim, Vector(self.kgrid.N), opt)
510 self.create_sensor_variables()
511 self.create_absorption_vars()
512 self.assign_pseudonyms(self.medium, self.kgrid)
513 self.scale_source_terms(opt.scale_source_terms)
514 self.create_pml_indices(
515 kgrid_dim=self.kgrid.dim,
516 kgrid_N=Vector(self.kgrid.N),
517 pml_size=pml_size,
518 pml_inside=opt.pml_inside,
519 is_axisymmetric=opt.simulation_type.is_axisymmetric(),
520 )
522 @staticmethod
523 def check_calling_func_name_and_dim(calling_func_name, kgrid_dim) -> None:
524 """
525 Check correct function has been called for the dimensionality of kgrid
527 Args:
528 calling_func_name: Name of the script that makes calls to kWaveSimulation
529 kgrid_dim: Dimensionality of the kWaveGrid
531 Returns:
532 None
533 """
534 assert not calling_func_name.startswith(("pstdElastic", "kspaceElastic")), "Elastic simulation is not supported."
536 if calling_func_name == "kspaceFirstOrder1D":
537 assert kgrid_dim == 1, f"kgrid has the wrong dimensionality for {calling_func_name}."
538 elif calling_func_name in ["kspaceFirstOrder2D", "pstdElastic2D", "kspaceElastic2D", "kspaceFirstOrderAS"]:
539 assert kgrid_dim == 2, f"kgrid has the wrong dimensionality for {calling_func_name}."
540 elif calling_func_name in ["kspaceFirstOrder3D", "pstdElastic3D", "kspaceElastic3D"]:
541 assert kgrid_dim == 3, f"kgrid has the wrong dimensionality for {calling_func_name}."
543 @staticmethod
544 def print_start_status(is_elastic_code: bool) -> None:
545 """
546 Update command-line status with the start time
548 Args:
549 is_elastic_code: is the simulation elastic
551 Returns:
552 None
553 """
554 if is_elastic_code: # pragma: no cover
555 logging.log(logging.INFO, "Running k-Wave elastic simulation...")
556 else:
557 logging.log(logging.INFO, "Running k-Wave simulation...")
558 logging.log(logging.INFO, f" start time: {get_date_string()}")
560 def set_index_data_type(self) -> None:
561 """
562 Pre-calculate the data type needed to store the matrix indices given the
563 total number of grid points: indexing variables will be created using this data type to save memory
565 Returns:
566 None
567 """
568 total_grid_points = self.kgrid.total_grid_points
569 self.index_data_type = get_smallest_possible_type(total_grid_points, "uint", default="double")
571 @staticmethod
572 def check_medium(medium, kgrid_k, simulation_type: SimulationType) -> bool:
573 """
574 Check the properties of the medium structure for correctness and validity
576 Args:
577 medium: kWaveMedium instance
578 kgrid_k: kWaveGrid.k matrix
579 is_elastic: Whether the simulation is elastic
580 is_axisymmetric: Whether the simulation is axisymmetric
582 Returns:
583 Medium Density
584 """
586 # if using the fluid code, allow the density field to be blank if the medium is homogeneous
587 if (not simulation_type.is_elastic_simulation()) and medium.density is None and medium.sound_speed.size == 1:
588 user_medium_density_input = False
589 medium.density = 1
590 else:
591 medium.ensure_defined("density")
592 user_medium_density_input = True
594 # check medium absorption inputs for the fluid code
595 is_absorbing = any(medium.is_defined("alpha_coeff", "alpha_power"))
596 is_stokes = simulation_type.is_axisymmetric() or medium.alpha_mode == "stokes"
597 medium.set_absorbing(is_absorbing, is_stokes)
599 if is_absorbing:
600 medium.check_fields(kgrid_k.shape)
601 return user_medium_density_input
603 def check_sensor(self, kgrid_dim) -> None:
604 """
605 Check the Sensor properties for correctness and validity
607 Args:
608 k_dim: kWaveGrid dimensionality
610 Returns:
611 None
612 """
613 # =========================================================================
614 # CHECK SENSOR STRUCTURE INPUTS
615 # =========================================================================
616 # check sensor fields
617 if self.sensor is not None:
618 # check the sensor input is valid
619 # TODO FARID move this check as a type checking
620 assert isinstance(
621 self.sensor, (kSensor, NotATransducer)
622 ), "sensor must be defined as an object of the kSensor or kWaveTransducer class."
624 # check if sensor is a transducer, otherwise check input fields
625 if not isinstance(self.sensor, NotATransducer):
626 if kgrid_dim == 2:
627 # check for sensor directivity input and set flag
628 directivity = self.sensor.directivity
629 if directivity is not None and self.sensor.directivity.angle is not None:
630 # make sure the sensor mask is not blank
631 assert self.sensor.mask is not None, "The mask must be defined for the sensor"
633 # check sensor.directivity.pattern and sensor.mask have the same size
634 assert (
635 directivity.angle.shape == self.sensor.mask.shape
636 ), "sensor.directivity.angle and sensor.mask must be the same size."
638 # check if directivity size input exists, otherwise make it
639 # a constant times kgrid.dx
640 if directivity.size is None:
641 directivity.set_default_size(self.kgrid)
643 # find the unique directivity angles
644 # assign the wavenumber vectors
645 directivity.set_unique_angles(self.sensor.mask)
646 directivity.set_wavenumbers(self.kgrid)
648 # check for time reversal inputs and set flags
649 if not self.options.simulation_type.is_elastic_simulation() and self.sensor.time_reversal_boundary_data is not None:
650 self.record.p = False
652 # check for sensor.record and set usage flgs - if no flgs are
653 # given, the time history of the acoustic pressure is recorded by
654 # default
655 if self.sensor.record is not None:
656 # check for time reversal data
657 if self.time_rev:
658 logging.log(logging.WARN, "sensor.record is not used for time reversal reconstructions")
660 # check the input is a cell array
661 assert isinstance(self.sensor.record, list), 'sensor.record must be given as a list, e.g. ["p", "u"]'
663 # check the sensor record flgs
664 self.record.set_flags_from_list(self.sensor.record, self.options.simulation_type.is_elastic_simulation())
666 # enforce the sensor.mask field unless just recording the max_all
667 # and _final variables
668 fields = ["p", "p_max", "p_min", "p_rms", "u", "u_non_staggered", "u_split_field", "u_max", "u_min", "u_rms", "I", "I_avg"]
669 if any(self.record.is_set(fields)):
670 assert self.sensor.mask is not None
672 # check if sensor mask is a binary grid, a set of cuboid corners,
673 # or a set of Cartesian interpolation points
674 if not self.blank_sensor:
675 if (kgrid_dim == 3 and num_dim2(self.sensor.mask) == 3) or (
676 kgrid_dim != 3 and (self.sensor.mask.shape == self.kgrid.k.shape)
677 ):
678 # check the grid is binary
679 assert self.sensor.mask.sum() == (
680 self.sensor.mask.size - (self.sensor.mask == 0).sum()
681 ), "sensor.mask must be a binary grid (numeric values must be 0 or 1)."
683 # check the grid is not empty
684 assert self.sensor.mask.sum() != 0, "sensor.mask must be a binary grid with at least one element set to 1."
686 elif self.sensor.mask.shape[0] == 2 * kgrid_dim:
687 # make sure the points are integers
688 assert np.all(self.sensor.mask % 1 == 0), "sensor.mask cuboid corner indices must be integers."
690 # store a copy of the cuboid corners
691 self.record.cuboid_corners_list = self.sensor.mask
693 # check the list makes sense
694 if np.any(self.sensor.mask[self.kgrid.dim :, :] - self.sensor.mask[: self.kgrid.dim, :] < 0):
695 if kgrid_dim == 1:
696 raise ValueError("sensor.mask cuboid corners must be defined " "as [x1, x2; ...]." " where x2 => x1, etc.")
697 elif kgrid_dim == 2:
698 raise ValueError(
699 "sensor.mask cuboid corners must be defined " "as [x1, y1, x2, y2; ...]." " where x2 => x1, etc."
700 )
701 elif kgrid_dim == 3:
702 raise ValueError(
703 "sensor.mask cuboid corners must be defined"
704 " as [x1, y1, z1, x2, y2, z2; ...]."
705 " where x2 => x1, etc."
706 )
708 # check the list are within bounds
709 if np.any(self.sensor.mask < 1):
710 raise ValueError("sensor.mask cuboid corners must be within the grid.")
711 else:
712 if kgrid_dim == 1:
713 if np.any(self.sensor.mask > self.kgrid.Nx):
714 raise ValueError("sensor.mask cuboid corners must be within the grid.")
715 elif kgrid_dim == 2:
716 if np.any(self.sensor.mask[[0, 2], :] > self.kgrid.Nx) or np.any(
717 self.sensor.mask[[1, 3], :] > self.kgrid.Ny
718 ):
719 raise ValueError("sensor.mask cuboid corners must be within the grid.")
720 elif kgrid_dim == 3:
721 if (
722 np.any(self.sensor.mask[[0, 3], :] > self.kgrid.Nx)
723 or np.any(self.sensor.mask[[1, 4], :] > self.kgrid.Ny)
724 or np.any(self.sensor.mask[[2, 5], :] > self.kgrid.Nz)
725 ):
726 raise ValueError("sensor.mask cuboid corners must be within the grid.")
728 # create a binary mask for display from the list of corners
729 # TODO FARID mask should be option_factory in sensor not here
730 self.sensor.mask = np.zeros_like(self.kgrid.k, dtype=bool)
731 cuboid_corners_list = self.record.cuboid_corners_list
732 for cuboid_index in range(cuboid_corners_list.shape[1]):
733 if self.kgrid.dim == 1:
734 self.sensor.mask[cuboid_corners_list[0, cuboid_index] : cuboid_corners_list[1, cuboid_index]] = 1
735 if self.kgrid.dim == 2:
736 self.sensor.mask[
737 cuboid_corners_list[0, cuboid_index] : cuboid_corners_list[2, cuboid_index],
738 cuboid_corners_list[1, cuboid_index] : cuboid_corners_list[3, cuboid_index],
739 ] = 1
740 if self.kgrid.dim == 3:
741 self.sensor.mask[
742 cuboid_corners_list[0, cuboid_index] : cuboid_corners_list[3, cuboid_index],
743 cuboid_corners_list[1, cuboid_index] : cuboid_corners_list[4, cuboid_index],
744 cuboid_corners_list[2, cuboid_index] : cuboid_corners_list[5, cuboid_index],
745 ] = 1
746 else:
747 # check the Cartesian sensor mask is the correct size
748 # (1 x N, 2 x N, 3 x N)
749 assert (
750 self.sensor.mask.shape[0] == kgrid_dim and num_dim2(self.sensor.mask) <= 2
751 ), f"Cartesian sensor.mask for a {kgrid_dim}D simulation must be given as a {kgrid_dim} by N array."
753 # set Cartesian mask flag (this is modified in
754 # createStorageVariables if the interpolation setting is
755 # set to nearest)
756 self.binary_sensor_mask = False
758 # extract Cartesian data from sensor mask
759 if kgrid_dim == 1:
760 # align sensor data as a column vector to be the
761 # same as kgrid.x_vec so that calls to interp1
762 # return data in the correct dimension
763 self.sensor_x = np.reshape((self.sensor.mask, (-1, 1)))
765 # add sensor_x to the record structure for use with
766 # the _extractSensorData subfunction
767 self.record.sensor_x = self.sensor_x
768 "record.sensor_x = sensor_x;"
770 elif kgrid_dim == 2:
771 self.sensor_x = self.sensor.mask[0, :]
772 self.sensor_y = self.sensor.mask[1, :]
773 elif kgrid_dim == 3:
774 self.sensor_x = self.sensor.mask[0, :]
775 self.sensor_y = self.sensor.mask[1, :]
776 self.sensor_z = self.sensor.mask[2, :]
778 # compute an equivalent sensor mask using nearest neighbour
779 # interpolation, if flgs.time_rev = false and
780 # cartesian_interp = 'linear' then this is only used for
781 # display, if flgs.time_rev = true or cartesian_interp =
782 # 'nearest' this grid is used as the sensor.mask
783 self.sensor.mask, self.order_index, self.reorder_index = cart2grid(
784 self.kgrid, self.sensor.mask, self.options.simulation_type.is_axisymmetric()
785 )
787 # if in time reversal mode, reorder the p0 input data in
788 # the order of the binary sensor_mask
789 if self.time_rev:
790 raise NotImplementedError
791 """
792 # append the reordering data
793 new_col_pos = length(sensor.time_reversal_boundary_data(1, :)) + 1;
794 sensor.time_reversal_boundary_data(:, new_col_pos) = order_index;
796 # reorder p0 based on the order_index
797 sensor.time_reversal_boundary_data = sort_rows(sensor.time_reversal_boundary_data, new_col_pos);
799 # remove the reordering data
800 sensor.time_reversal_boundary_data = sensor.time_reversal_boundary_data(:, 1:new_col_pos - 1);
801 """
802 else:
803 # set transducer sensor flag
804 self.transducer_sensor = True
805 self.record.p = False
807 # check to see if there is an elevation focus
808 if not np.isinf(self.sensor.elevation_focus_distance):
809 # set flag
810 self.transducer_receive_elevation_focus = True
812 # get the elevation mask that is used to extract the correct values
813 # from the sensor data buffer for averaging
814 self.transducer_receive_mask = self.sensor.elevation_beamforming_mask
816 # check for directivity inputs with time reversal
817 if kgrid_dim == 2 and self.use_sensor and self.compute_directivity and self.time_rev:
818 logging.log(logging.WARN, "sensor directivity fields are not used for time reversal.")
820 def check_source(self, k_dim, k_Nt) -> None:
821 """
822 Check the source properties for correctness and validity
824 Args:
825 kgrid_dim: kWaveGrid dimension
826 k_Nt: Number of time steps in kWaveGrid
828 Returns:
829 None
830 """
831 # =========================================================================
832 # CHECK SOURCE STRUCTURE INPUTS
833 # =========================================================================
835 # check source inputs
836 if not isinstance(self.source, (kSource, NotATransducer)):
837 # allow an invalid or empty source input if computing time reversal,
838 # otherwise return error
839 assert self.time_rev, "source must be defined as an object of the kSource or kWaveTransducer classes."
841 elif not isinstance(self.source, NotATransducer):
842 # --------------------------
843 # SOURCE IS NOT A TRANSDUCER
844 # --------------------------
846 """
847 check allowable source types
849 Depending on the kgrid dimensionality and the simulation type,
850 following fields are allowed & might be use:
852 kgrid.dim == 1:
853 non-elastic code:
854 ['p0', 'p', 'p_mask', 'p_mode', 'p_frequency_ref', 'ux', 'u_mask', 'u_mode', 'u_frequency_ref']
855 kgrid.dim == 2:
856 non-elastic code:
857 ['p0', 'p', 'p_mask', 'p_mode', 'p_frequency_ref', 'ux', 'uy', 'u_mask', 'u_mode', 'u_frequency_ref']
858 elastic code:
859 ['p0', 'sxx', 'syy', 'sxy', 's_mask', 's_mode', 'ux', 'uy', 'u_mask', 'u_mode']
860 kgrid.dim == 3:
861 non-elastic code:
862 ['p0', 'p', 'p_mask', 'p_mode', 'p_frequency_ref', 'ux', 'uy', 'uz', 'u_mask', 'u_mode', 'u_frequency_ref']
863 elastic code:
864 ['p0', 'sxx', 'syy', 'szz', 'sxy', 'sxz', 'syz', 's_mask', 's_mode', 'ux', 'uy', 'uz', 'u_mask', 'u_mode']
865 """
867 self.source.validate(self.kgrid)
869 # check for a time varying pressure source input
870 if self.source.p is not None:
871 # check the source mode input is valid
872 if self.source.p_mode is None:
873 self.source.p_mode = self.SOURCE_P_MODE_DEF
875 if self.source_p > k_Nt:
876 logging.log(logging.WARN, " source.p has more time points than kgrid.Nt, remaining time points will not be used.")
878 # create an indexing variable corresponding to the location of all the source elements
879 self.p_source_pos_index = matlab_find(self.source.p_mask)
881 # check if the mask is binary or labelled
882 p_unique = np.unique(self.source.p_mask)
884 # create a second indexing variable
885 if p_unique.size <= 2 and p_unique.sum() == 1:
886 # set signal index to all elements
887 self.p_source_sig_index = ":"
888 else:
889 # set signal index to the labels (this allows one input signal
890 # to be used for each source label)
891 self.p_source_sig_index = self.source.p_mask(self.source.p_mask != 0)
893 # convert the data type depending on the number of indices
894 self.p_source_pos_index = cast_to_type(self.p_source_pos_index, self.index_data_type)
895 if self.source_p_labelled:
896 self.p_source_sig_index = cast_to_type(self.p_source_sig_index, self.index_data_type)
898 # check for time varying velocity source input and set source flag
899 if any([(getattr(self.source, k) is not None) for k in ["ux", "uy", "uz", "u_mask"]]):
900 # check the source mode input is valid
901 if self.source.u_mode is None:
902 self.source.u_mode = self.SOURCE_U_MODE_DEF
904 # create an indexing variable corresponding to the location of all
905 # the source elements
906 self.u_source_pos_index = matlab_find(self.source.u_mask)
908 # check if the mask is binary or labelled
909 u_unique = np.unique(self.source.u_mask)
911 # create a second indexing variable
912 if u_unique.size <= 2 and u_unique.sum() == 1:
913 # set signal index to all elements
914 self.u_source_sig_index = ":"
915 else:
916 # set signal index to the labels (this allows one input signal
917 # to be used for each source label)
918 self.u_source_sig_index = self.source.u_mask[self.source.u_mask != 0]
920 # convert the data type depending on the number of indices
921 self.u_source_pos_index = cast_to_type(self.u_source_pos_index, self.index_data_type)
922 if self.source_u_labelled:
923 self.u_source_sig_index = cast_to_type(self.u_source_sig_index, self.index_data_type)
925 # check for time varying stress source input and set source flag
926 if any([(getattr(self.source, k) is not None) for k in ["sxx", "syy", "szz", "sxy", "sxz", "syz", "s_mask"]]):
927 # create an indexing variable corresponding to the location of all
928 # the source elements
929 raise NotImplementedError
930 "s_source_pos_index = find(source.s_mask != 0);"
932 # check if the mask is binary or labelled
933 "s_unique = unique(source.s_mask);"
935 # create a second indexing variable
936 if eng.eval("numel(s_unique) <= 2 && sum(s_unique) == 1"): # noqa: F821
937 # set signal index to all elements
938 eng.workspace["s_source_sig_index"] = ":" # noqa: F821
940 else:
941 # set signal index to the labels (this allows one input signal
942 # to be used for each source label)
943 s_source_sig_index = source.s_mask(source.s_mask != 0) # noqa
945 f"s_source_pos_index = {self.index_data_type}(s_source_pos_index);"
946 if self.source_s_labelled:
947 f"s_source_sig_index = {self.index_data_type}(s_source_sig_index);"
949 else:
950 # ----------------------
951 # SOURCE IS A TRANSDUCER
952 # ----------------------
954 # if the sensor is a transducer, check that the simulation is in 3D
955 assert k_dim == 3, "Transducer inputs are only compatible with 3D simulations."
957 # get the input signal - this is appended with zeros if required to
958 # account for the beamforming delays (this will throw an error if the
959 # input signal is not defined)
960 self.transducer_input_signal = self.source.input_signal
962 # get the delay mask that accounts for the beamforming delays and
963 # elevation focussing; this is used so that a single time series can be
964 # applied to the complete transducer mask with different delays
965 delay_mask = self.source.delay_mask()
967 # set source flag - this should be the length of signal minus the
968 # maximum delay
969 self.transducer_source = self.transducer_input_signal.size - delay_mask.max()
971 # get the active elements mask
972 active_elements_mask = self.source.active_elements_mask
974 # get the apodization mask if not set to 'Rectangular' and convert to a
975 # linear array
976 if self.source.transmit_apodization == "Rectangular":
977 self.transducer_transmit_apodization = 1
978 else:
979 self.transducer_transmit_apodization = self.source.transmit_apodization_mask
980 self.transducer_transmit_apodization = self.transducer_transmit_apodization[active_elements_mask != 0]
982 # create indexing variable corresponding to the active elements
983 # and convert the data type depending on the number of indices
984 self.u_source_pos_index = matlab_find(active_elements_mask).astype(self.index_data_type)
986 # convert the delay mask to an indexing variable (this doesn't need to
987 # be modified if the grid is expanded) which tells each point in the
988 # source mask which point in the input_signal should be used
989 delay_mask = matlab_mask(delay_mask, active_elements_mask != 0) # compatibility
991 # convert the data type depending on the maximum value of the delay
992 # mask and the length of the source
993 smallest_type = get_smallest_possible_type(delay_mask.max(), "uint")
994 if smallest_type is not None:
995 delay_mask = delay_mask.astype(smallest_type)
997 # move forward by 1 as a delay of 0 corresponds to the first point in the input signal
998 self.delay_mask = delay_mask + 1
1000 # clean up unused variables
1001 del active_elements_mask
1003 def check_kgrid_time(self) -> None:
1004 """
1005 Check time-related kWaveGrid inputs
1007 Returns:
1008 None
1009 """
1011 # check kgrid for t_array existence, and create if not defined
1012 if isinstance(self.kgrid.t_array, str) and self.kgrid.t_array == "auto":
1013 # check for time reversal mode
1014 if self.time_rev:
1015 raise ValueError("kgrid.t_array (Nt and dt) must be defined explicitly in time reversal mode.")
1017 # check for time varying sources
1018 if (not self.source_p0_elastic) and (
1019 self.source_p
1020 or self.source_ux
1021 or self.source_uy
1022 or self.source_uz
1023 or self.source_sxx
1024 or self.source_syy
1025 or self.source_szz
1026 or self.source_sxy
1027 or self.source_sxz
1028 or self.source_syz
1029 ):
1030 raise ValueError("kgrid.t_array (Nt and dt) must be defined explicitly when using a time-varying source.")
1032 # create the time array using the compressional sound speed
1033 self.kgrid.makeTime(self.medium.sound_speed, self.KSPACE_CFL)
1035 # check kgrid.t_array for stability given medium properties
1036 if not self.options.simulation_type.is_elastic_simulation():
1037 # calculate the largest timestep for which the model is stable
1039 dt_stability_limit = check_stability(self.kgrid, self.medium)
1041 # give a warning if the timestep is larger than stability limit allows
1042 if self.kgrid.dt > dt_stability_limit:
1043 logging.log(logging.WARN, " time step may be too large for a stable simulation.")
1045 @staticmethod
1046 def select_precision(opt: SimulationOptions):
1047 """
1048 Select the minimal precision for storing the data
1050 Args:
1051 opt: SimulationOptions instance
1053 Returns:
1054 Minimal precision for variable allocation
1056 """
1057 # set storage variable type based on data_cast - this enables the
1058 # output variables to be directly created in the data_cast format,
1059 # rather than creating them in double precision and then casting them
1060 if opt.data_cast == "off":
1061 precision = "double"
1062 elif opt.data_cast == "single":
1063 precision = "single"
1064 elif opt.data_cast == "gsingle":
1065 precision = "single"
1066 elif opt.data_cast == "gdouble":
1067 precision = "double"
1068 elif opt.data_cast == "gpuArray":
1069 raise NotImplementedError("gpuArray is not supported in Python-version")
1070 elif opt.data_cast == "kWaveGPUsingle":
1071 precision = "single"
1072 elif opt.data_cast == "kWaveGPUdouble":
1073 precision = "double"
1074 else:
1075 raise ValueError("'Unknown ''DataCast'' option'")
1076 return precision
1078 def check_input_combinations(self, opt: SimulationOptions, user_medium_density_input: bool, k_dim, pml_size, kgrid_N) -> None:
1079 """
1080 Check the input combinations for correctness and validity
1082 Args:
1083 opt: SimulationOptions instance
1084 user_medium_density_input: Medium Density
1085 k_dim: kWaveGrid dimensionality
1086 pml_size: Size of the PML
1087 kgrid_N: kWaveGrid size in each direction
1089 Returns:
1090 None
1091 """
1092 # =========================================================================
1093 # CHECK FOR VALID INPUT COMBINATIONS
1094 # =========================================================================
1096 # enforce density input if velocity sources or output are being used
1097 if not user_medium_density_input and (
1098 self.source_ux or self.source_uy or self.source_uz or self.record.u or self.record.u_max or self.record.u_rms
1099 ):
1100 raise ValueError(
1101 "medium.density must be explicitly defined " "if velocity inputs or outputs are used, even in homogeneous media."
1102 )
1104 # TODO(walter): move to check medium
1105 # enforce density input if nonlinear equations are being used
1106 if not user_medium_density_input and self.medium.is_nonlinear():
1107 raise ValueError("medium.density must be explicitly defined if medium.BonA is specified.")
1109 # check sensor compatibility options for flgs.compute_directivity
1110 if self.use_sensor and k_dim == 2 and self.compute_directivity and not self.binary_sensor_mask and opt.cartesian_interp == "linear":
1111 raise ValueError(
1112 "sensor directivity fields are only compatible " "with binary sensor masks or " "CartInterp" " set to " "nearest" "."
1113 )
1115 # check for split velocity output
1116 if self.record.u_split_field and not self.binary_sensor_mask:
1117 raise ValueError("The option sensor.record = {" "u_split_field" "} is only compatible " "with a binary sensor mask.")
1119 # check input options for data streaming *****
1120 if opt.stream_to_disk:
1121 if not self.use_sensor or self.time_rev:
1122 raise ValueError(
1123 "The optional input "
1124 "StreamToDisk"
1125 " is currently only compatible "
1126 "with forward simulations using a non-zero sensor mask."
1127 )
1128 elif self.sensor.record is not None and self.sensor.record.ismember(self.record.flags[1:]).any():
1129 raise ValueError(
1130 "The optional input " "StreamToDisk" " is currently only compatible " "with sensor.record = {" "p" "} (the default)."
1131 )
1133 is_axisymmetric = self.options.simulation_type.is_axisymmetric()
1134 # make sure the PML size is smaller than the grid if PMLInside is true
1135 if opt.pml_inside and (
1136 (k_dim == 1 and (pml_size.x * 2 > self.kgrid.Nx))
1137 or (k_dim == 2 and not is_axisymmetric and ((pml_size.x * 2 > kgrid_N[0]) or (pml_size.y * 2 > kgrid_N[1])))
1138 or (k_dim == 2 and is_axisymmetric and ((pml_size.x * 2 > kgrid_N[0]) or (pml_size.y > kgrid_N[1])))
1139 or (k_dim == 3 and ((pml_size.x * 2 > kgrid_N[0]) or (pml_size.x * 2 > kgrid_N[1]) or (pml_size.z * 2 > kgrid_N[2])))
1140 ):
1141 raise ValueError("The size of the PML must be smaller than the size of the grid.")
1143 # make sure the PML is inside if using a non-uniform grid
1144 if self.nonuniform_grid and not opt.pml_inside:
1145 raise ValueError("''PMLInside'' must be true for simulations using non-uniform grids.")
1147 # check for compatible input options if saving to disk
1148 # modified by Farid | disabled temporarily!
1149 # if k_dim == 3 and isinstance(self.options.save_to_disk, str) and
1150 # (not self.use_sensor or not self.binary_sensor_mask or self.time_rev):
1151 # raise ValueError('The optional input ''SaveToDisk'' is currently only compatible
1152 # with forward simulations using a non-zero binary sensor mask.')
1154 # check the record start time is within range
1155 record_start_index = self.sensor.record_start_index
1156 if self.use_sensor and ((record_start_index > self.kgrid.Nt) or (record_start_index < 1)):
1157 raise ValueError("sensor.record_start_index must be between 1 and the number of time steps.")
1159 # ensure 'WSWA' symmetry if using axisymmetric code with 'SaveToDisk'
1160 if is_axisymmetric and self.options.radial_symmetry != "WSWA" and isinstance(self.options.save_to_disk, str):
1161 # display a warning only if using WSWS symmetry (not WSWA-FFT)
1162 if self.options.radial_symmetry.startswith("WSWS"):
1163 logging.log(
1164 logging.WARN, " Optional input " "RadialSymmetry" " changed to " "WSWA" " for compatibility with " "SaveToDisk" "."
1165 )
1167 # update setting
1168 self.options.radial_symmetry = "WSWA"
1170 # ensure p0 smoothing is switched off if p0 is empty
1171 if not self.source_p0:
1172 self.options.smooth_p0 = False
1174 # start log if required
1175 if opt.create_log:
1176 raise NotImplementedError(f"diary({self.LOG_NAME}.txt');")
1178 # update command line status
1179 if self.time_rev:
1180 logging.log(logging.INFO, " time reversal mode")
1182 # cleanup unused variables
1183 for k in list(self.__dict__.keys()):
1184 if k.endswith("_DEF"):
1185 delattr(self, k)
1187 def smooth_and_enlarge(self, source, k_dim, kgrid_N, opt: SimulationOptions) -> None:
1188 """
1189 Smooth and enlarge grids
1191 Args:
1192 source: kWaveSource instance
1193 k_dim: kWaveGrid dimensionality
1194 kgrid_N: kWaveGrid size in each direction
1195 opt: SimulationOptions
1197 Returns:
1198 None
1199 """
1201 # smooth the initial pressure distribution p0 if required, and then restore
1202 # the maximum magnitude
1203 # NOTE 1: if p0 has any values at the edge of the domain, the smoothing
1204 # may cause part of p0 to wrap to the other side of the domain
1205 # NOTE 2: p0 is smoothed before the grid is expanded to ensure that p0 is
1206 # exactly zero within the PML
1207 # NOTE 3: for the axisymmetric code, p0 is smoothed assuming WS origin
1208 # symmetry
1209 if self.source_p0 and self.options.smooth_p0:
1210 # update command line status
1211 logging.log(logging.INFO, " smoothing p0 distribution...")
1213 if self.options.simulation_type.is_axisymmetric():
1214 if self.options.radial_symmetry in ["WSWA-FFT", "WSWA"]:
1215 # create a new kWave grid object with expanded radial grid
1216 kgrid_exp = kWaveGrid([kgrid_N.x, kgrid_N.y * 4], [self.kgrid.dx, self.kgrid.dy])
1218 # mirror p0 in radial dimension using WSWA symmetry
1219 self.source.p0 = self.source.p0.astype(float)
1220 p0_exp = np.zeros((kgrid_exp.Nx, kgrid_exp.Ny))
1221 p0_exp[:, kgrid_N.y * 0 + 0 : kgrid_N.y * 1] = self.source.p0
1222 p0_exp[:, kgrid_N.y * 1 + 1 : kgrid_N.y * 2] = -np.fliplr(self.source.p0[:, 1:])
1223 p0_exp[:, kgrid_N.y * 2 + 0 : kgrid_N.y * 3] = -self.source.p0
1224 p0_exp[:, kgrid_N.y * 3 + 1 : kgrid_N.y * 4] = np.fliplr(self.source.p0[:, 1:])
1226 elif self.options.radial_symmetry in ["WSWS-FFT", "WSWS"]:
1227 # create a new kWave grid object with expanded radial grid
1228 kgrid_exp = kWaveGrid([kgrid_N.x, kgrid_N.y * 2 - 2], [self.kgrid.dx, self.kgrid.dy])
1230 # mirror p0 in radial dimension using WSWS symmetry
1231 p0_exp = np.zeros((kgrid_exp.Nx, kgrid_exp.Ny))
1232 p0_exp[:, 1 : kgrid_N.y] = source.p0
1233 p0_exp[:, kgrid_N.y + 0 : kgrid_N.y * 2 - 2] = np.fliplr(source.p0[:, 1:-1])
1235 # smooth p0
1236 p0_exp = smooth(p0_exp, True)
1238 # trim back to original size
1239 source.p0 = p0_exp[:, 0 : self.kgrid.Ny]
1241 # clean up unused variables
1242 del kgrid_exp
1243 del p0_exp
1244 else:
1245 source.p0 = smooth(source.p0, True)
1247 # expand the computational grid if the PML is set to be outside the input
1248 # grid defined by the user
1249 if opt.pml_inside is False:
1250 expand_results = expand_grid_matrices(
1251 self.kgrid,
1252 self.medium,
1253 self.source,
1254 self.sensor,
1255 self.options,
1256 dotdict(
1257 {
1258 "p_source_pos_index": self.p_source_pos_index,
1259 "u_source_pos_index": self.u_source_pos_index,
1260 "s_source_pos_index": self.s_source_pos_index,
1261 }
1262 ),
1263 dotdict(
1264 {
1265 "axisymmetric": self.options.simulation_type.is_axisymmetric(),
1266 "use_sensor": self.use_sensor,
1267 "blank_sensor": self.blank_sensor,
1268 "cuboid_corners": self.cuboid_corners,
1269 "source_p0": self.source_p0,
1270 "source_p": self.source_p,
1271 "source_ux": self.source_ux,
1272 "source_uy": self.source_uy,
1273 "source_uz": self.source_uz,
1274 "transducer_source": self.transducer_source,
1275 "source_sxx": self.source_sxx,
1276 "source_syy": self.source_syy,
1277 "source_szz": self.source_szz,
1278 "source_sxy": self.source_sxy,
1279 "source_sxz": self.source_sxz,
1280 "source_syz": self.source_syz,
1281 }
1282 ),
1283 )
1284 self.kgrid, self.index_data_type, self.p_source_pos_index, self.u_source_pos_index, self.s_source_pos_index = expand_results
1286 # get maximum prime factors
1287 if self.options.simulation_type.is_axisymmetric():
1288 prime_facs = self.kgrid.highest_prime_factors(self.options.radial_symmetry[:4])
1289 else:
1290 prime_facs = self.kgrid.highest_prime_factors()
1292 # give warning for bad dimension sizes
1293 if prime_facs.max() > self.HIGHEST_PRIME_FACTOR_WARNING:
1294 prime_facs = prime_facs[prime_facs != 0]
1295 logging.log(logging.WARN, f"Highest prime factors in each dimension are {prime_facs}")
1296 logging.log(logging.WARN, "Use dimension sizes with lower prime factors to improve speed")
1297 del prime_facs
1299 # smooth the sound speed distribution if required
1300 if opt.smooth_c0 and num_dim2(self.medium.sound_speed) == k_dim and self.medium.sound_speed.size > 1:
1301 logging.log(logging.INFO, " smoothing sound speed distribution...")
1302 self.medium.sound_speed = smooth(self.medium.sound_speed)
1304 # smooth the ambient density distribution if required
1305 if opt.smooth_rho0 and num_dim2(self.medium.density) == k_dim and self.medium.density.size > 1:
1306 logging.log(logging.INFO, "smoothing density distribution...")
1307 self.medium.density = smooth(self.medium.density)
1309 def create_sensor_variables(self) -> None:
1310 """
1311 Create the sensor related variables
1313 Returns:
1314 None
1315 """
1316 # define the output variables and mask indices if using the sensor
1317 if self.use_sensor:
1318 if not self.blank_sensor or self.options.save_to_disk:
1319 if self.cuboid_corners:
1320 # create empty list of sensor indices
1321 self.sensor_mask_index = []
1323 # loop through the list of cuboid corners, and extract the
1324 # sensor mask indices for each cube
1325 for cuboid_index in range(self.record.cuboid_corners_list.shape[1]):
1326 # create empty binary mask
1327 temp_mask = np.zeros_like(self.kgrid.k, dtype=bool)
1329 if self.kgrid.dim == 1:
1330 self.sensor.mask[
1331 self.record.cuboid_corners_list[0, cuboid_index] : self.record.cuboid_corners_list[1, cuboid_index]
1332 ] = 1
1333 if self.kgrid.dim == 2:
1334 self.sensor.mask[
1335 self.record.cuboid_corners_list[0, cuboid_index] : self.record.cuboid_corners_list[2, cuboid_index],
1336 self.record.cuboid_corners_list[1, cuboid_index] : self.record.cuboid_corners_list[3, cuboid_index],
1337 ] = 1
1338 if self.kgrid.dim == 3:
1339 self.sensor.mask[
1340 self.record.cuboid_corners_list[0, cuboid_index] : self.record.cuboid_corners_list[3, cuboid_index],
1341 self.record.cuboid_corners_list[1, cuboid_index] : self.record.cuboid_corners_list[4, cuboid_index],
1342 self.record.cuboid_corners_list[2, cuboid_index] : self.record.cuboid_corners_list[5, cuboid_index],
1343 ] = 1
1345 # extract mask indices
1346 self.sensor_mask_index.append(matlab_find(temp_mask))
1347 self.sensor_mask_index = np.array(self.sensor_mask_index)
1349 # cleanup unused variables
1350 del temp_mask
1352 else:
1353 # create mask indices (this works for both normal sensor and
1354 # transducer inputs)
1355 self.sensor_mask_index = np.where(self.sensor.mask.flatten(order="F") != 0)[0] + 1 # +1 due to matlab indexing
1356 self.sensor_mask_index = np.expand_dims(self.sensor_mask_index, -1) # compatibility, n => [n, 1]
1358 # convert the data type depending on the number of indices (this saves
1359 # memory)
1360 self.sensor_mask_index = cast_to_type(self.sensor_mask_index, self.index_data_type)
1362 else:
1363 # set the sensor mask index variable to be empty
1364 self.sensor_mask_index = []
1366 def create_absorption_vars(self) -> None:
1367 """
1368 Create absorption variables for the fluid code based on
1369 the expanded and smoothed values of the medium parameters (if not saving to disk)
1371 Returns:
1372 None
1373 """
1374 if not self.options.simulation_type.is_elastic_simulation() and not self.options.save_to_disk:
1375 self.absorb_nabla1, self.absorb_nabla2, self.absorb_tau, self.absorb_eta = create_absorption_variables(
1376 self.kgrid, self.medium, self.equation_of_state
1377 )
1379 def assign_pseudonyms(self, medium: kWaveMedium, kgrid: kWaveGrid) -> None:
1380 """
1381 Shorten commonly used field names (these act only as pointers provided that the values aren't modified)
1382 (done after enlarging and smoothing the grids)
1384 Args:
1385 medium: kWaveMedium instance
1386 kgrid: kWaveGrid instance
1388 Returns:
1389 None
1390 """
1391 self.dt = float(kgrid.dt)
1392 self.rho0 = medium.density
1393 self.c0 = medium.sound_speed
1395 def scale_source_terms(self, is_scale_source_terms) -> None:
1396 """
1397 Scale the source terms based on the expanded and smoothed values of the medium parameters
1399 Args:
1400 is_scale_source_terms: Should the source terms be scaled
1402 Returns:
1403 None
1404 """
1405 if not is_scale_source_terms:
1406 return
1407 try:
1408 p_source_pos_index = self.p_source_pos_index
1409 except AttributeError:
1410 p_source_pos_index = None
1412 try:
1413 s_source_pos_index = self.s_source_pos_index
1414 except AttributeError:
1415 s_source_pos_index = None
1417 try:
1418 u_source_pos_index = self.u_source_pos_index
1419 except AttributeError:
1420 u_source_pos_index = None
1422 self.transducer_input_signal = scale_source_terms_func(
1423 self.c0,
1424 self.dt,
1425 self.kgrid,
1426 self.source,
1427 p_source_pos_index,
1428 s_source_pos_index,
1429 u_source_pos_index,
1430 self.transducer_input_signal,
1431 dotdict(
1432 {
1433 "nonuniform_grid": self.nonuniform_grid,
1434 "source_ux": self.source_ux,
1435 "source_uy": self.source_uy,
1436 "source_uz": self.source_uz,
1437 "transducer_source": self.transducer_source,
1438 "source_p": self.source_p,
1439 "source_p0": self.source_p0,
1440 "use_w_source_correction_p": self.use_w_source_correction_p,
1441 "use_w_source_correction_u": self.use_w_source_correction_u,
1442 "source_sxx": self.source_sxx,
1443 "source_syy": self.source_syy,
1444 "source_szz": self.source_szz,
1445 "source_sxy": self.source_sxy,
1446 "source_sxz": self.source_sxz,
1447 "source_syz": self.source_syz,
1448 }
1449 ),
1450 )
1452 def create_pml_indices(self, kgrid_dim, kgrid_N: Vector, pml_size, pml_inside, is_axisymmetric):
1453 """
1454 Define index variables to remove the PML from the display if the optional
1455 input 'PlotPML' is set to false
1457 Args:
1458 kgrid_dim: kWaveGrid dimensionality
1459 kgrid_N: kWaveGrid size in each direction
1460 pml_size: Size of the PML
1461 pml_inside: Whether the PML is inside the grid defined by the user
1462 is_axisymmetric: Whether the simulation is axisymmetric
1464 """
1465 # comment by Farid: PlotPML is always False in Python version,
1466 # therefore if statement removed
1467 if kgrid_dim == 1:
1468 self.x1 = pml_size.x + 1.0
1469 self.x2 = kgrid_N.x - pml_size.x
1470 elif kgrid_dim == 2:
1471 self.x1 = pml_size.x + 1.0
1472 self.x2 = kgrid_N.x - pml_size.x
1473 if is_axisymmetric:
1474 self.y1 = 1.0
1475 else:
1476 self.y1 = pml_size.y + 1.0
1477 self.y2 = kgrid_N.y - pml_size.y
1478 elif kgrid_dim == 3:
1479 self.x1 = pml_size.x + 1.0
1480 self.x2 = kgrid_N.x - pml_size.x
1481 self.y1 = pml_size.y + 1.0
1482 self.y2 = kgrid_N.y - pml_size.y
1483 self.z1 = pml_size.z + 1.0
1484 self.z2 = kgrid_N.z - pml_size.z
1486 # define index variables to allow original grid size to be maintained for
1487 # the _final and _all output variables if 'PMLInside' is set to false
1488 # if self.record is None:
1489 # self.record = Recorder()
1490 self.record.set_index_variables(self.kgrid, pml_size, pml_inside, is_axisymmetric)
1492 @deprecated(version="0.4.1", reason="Use TimeReversal class instead")
1493 def check_time_reversal(self) -> bool:
1494 return self.time_rev