Hi Bradley,

I encountered a situation similar to the one Toni described above. In my simulation, I added a collimator into the medium. With `kspaceFirstOrder3DC`

or `kspaceFirstOrder3DG`

, I got all NaN in the sensor grid as long as I added the collimator. While with common `kspaceFirstOrder3D`

, whether casting data to `single`

or `gpuArray-single`

, the results were normal. Meanwhile, without adding the collimator, `kspaceFirstOrder3DC`

and `kspaceFirstOrder3DG`

work fine with pure water medium.

I'm not sure how this would happen as I have used the function `checkStability`

to make sure my dt (3.3333e-08) is smaller than the dt_stability_limit (5.0456e-08).

Thank you very much for taking a look into this!

Best,

Artery

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% =========================================================================

% DEFINE LITERALS

% =========================================================================

% medium parameters

c0 = 1500; % sound speed [m/s]

rho0 = 1000; % density [kg/m^3]

alpha0 = 3.0227e-05; % water power law absorption coefficient [dB/(MHz^y cm)]

% collimator parameters

c_collimator = 2410;

rho_collimator = 1160;

alpha_collimator= 9.54;

% source parameters

source_f0 = 5e5; % source frequency [Hz]

source_roc = 59.3e-3; % bowl radius of curvature [m]

source_diameter = 25.4e-3; % bowl aperture diameter [m]

source_amp = 1.4*2e4; % source pressure [Pa]

focus_depth = 38.1e-3; % [m]

% grid parameters

axial_size = 70e-3; % total grid size in the axial dimension [m]

lateral_size = 70e-3; % total grid size in the lateral dimension [m]

% computational parameters

ppw = 12; % number of points per wavelength

t_end = 200e-6; % total compute time [s]

record_periods = 1; % number of periods to record

cfl = 0.2; % CFL number

source_x_offset = 20; % grid points to offset the source

bli_tolerance = 0.01; % tolerance for truncation of the off-grid source points

upsampling_rate = 10; % density of integration points relative to grid

% =========================================================================

% DEFINE KGRID

% =========================================================================

% calculate the grid spacing based on the PPW and F0

dx = c0 / (ppw * source_f0); % [m]

% compute the size of the grid

Nx = roundEven(axial_size / dx) + source_x_offset;

Ny = roundEven(lateral_size / dx);

Nz = Ny;

% create the computational grid

kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);

% compute points per temporal period

PPP = round(ppw / cfl);

% compute corresponding time spacing

dt = 1 / (PPP * source_f0);

% create the time array using an integer number of points per period

Nt = round(t_end / dt);

kgrid.setTime(Nt, dt);

% =========================================================================

% DEFINE KWAVEARRAY FOR TRANDUCER

% =========================================================================

% set bowl position and orientation

bowl_pos = [kgrid.x_vec(1) + source_x_offset * kgrid.dx, 0, 0];

focus_pos = [bowl_pos(1) + focus_depth, 0, 0];

% create empty kWaveArray

karray = kWaveArray('BLITolerance', bli_tolerance, 'UpsamplingRate', upsampling_rate, 'SinglePrecision', true);

% add bowl shaped element

karray.addBowlElement(bowl_pos, source_roc, source_diameter, focus_pos);

% =========================================================================

% DEFINE SOURCE SIGNAL

% =========================================================================

% create time varying source

source_sig = createCWSignals(kgrid.t_array, source_f0, source_amp, 0);

% assign binary mask

source.p_mask = karray.getArrayBinaryMask(kgrid);

% assign source signals

source.p = karray.getDistributedSourceSignal(kgrid, source_sig);

% =========================================================================

% DEFINE BINARY COLLIMATOR MASK

% =========================================================================

[OUTPUTgrid] = VOXELISE(166,166,152,'Collimator_14mm.STL','xyz');

OUTPUTgrid = imrotate3(double(OUTPUTgrid), 90, [1 0 0]);

OUTPUTgrid = logical(OUTPUTgrid);

collimator_mask = zeros(Nx, Ny, Nz, 'logical');

collimator_mask(ceil(Nx/2)-ceil(size(OUTPUTgrid,1)/2)-59-14:ceil(Nx/2)+ceil(size(OUTPUTgrid,1)/2)-1-59-14, ...

ceil(Ny/2)-ceil(size(OUTPUTgrid,2)/2)+1:ceil(Ny/2)+ceil(size(OUTPUTgrid,2)/2)-1+1, ...

ceil(Nz/2)-ceil(size(OUTPUTgrid,3)/2)+2:ceil(Nz/2)+ceil(size(OUTPUTgrid,3)/2)-1+2) = OUTPUTgrid;

% =========================================================================

% DEFINE MEDIUM

% =========================================================================

% assign medium properties

medium.sound_speed = c0.*ones(Nx, Ny, Nz);

medium.sound_speed(collimator_mask) = c_collimator;

medium.density = rho0.*ones(Nx, Ny, Nz);

medium.density(collimator_mask) = rho_collimator;

medium.alpha_coeff = alpha0.*ones(Nx, Ny, Nz);

medium.alpha_coeff(collimator_mask) = alpha_collimator;

medium.alpha_power = 1.001;

medium.alpha_mode = 'no_dispersion';

% =========================================================================

% DEFINE SENSOR

% =========================================================================

% set sensor mask to record central plane, not including the source point

sensor.mask = zeros(Nx, Ny, Nz);

sensor.mask(source_x_offset + 2:end, :, Nz/2 + 1) = 1; % mid-plane/central plane

% record the pressure

sensor.record = {'p', 'p_rms'};

% record only the final few periods when the field is in steady state

sensor.record_start_index = kgrid.Nt - record_periods * PPP + 1;

% =========================================================================

% SIMULATION

% =========================================================================

% set input options

input_args = {...

'PMLSize', 'auto', ...

'PMLInside', false, ...

'PlotPML', false, ...

'DisplayMask', 'off'};

% run code

sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%