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{:});
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%