Hi Brad,
I didn't find a solution in the meantime. I even switched off smoothing, and smoothed manually. Same result. I replaced the external source file by a simple ball, and got the same instability again. So yes, I can create an example that doesn't depend on external functions or data that has the same problem. The only thing you'll have to change is the path of the created HDF5 input file (if you want to create one at all). The code follows...
%% Input parameters
grid_length = 0.1; % length of kgrid in m
grid_width = 0.1; % width of kgrid in m
axial_spacing = 250e-6;
% spacing between neighbouring points in axial direction (dx) in m
lateral_spacing = 250e-6;
% spacing between neighbouring points in lateral direction (dy,dz) in m
time_step = 25e-9; % simulation time step in s
total_time = 35e-6; % total simulation time in s
T = 25; % temperature in °C
d_geant_grid = 100e-6; % grid spacing of the above source file
smooth_p0 = true;
% Boolean controlling whether initial pressure distribution gets smoothed automatically
smooth_density = true;
% Boolean controlling whether density distribution gets smoothed automatically
smooth_soundspeed = true;
% Boolean controlling whether sound speed distribution gets smoothed automatically
StreamToDisk = true;
% Boolean controlling whether simulation data are stored on disk during simulation
DataCast = 'gpuArray-double';
% Variable precision & simulation mode. Possible: 'single', 'double', 'gpuArray-double', 'gpuArray-single'
SaveToDisk = 'sirmiobeam-double-250-in';
% File name of the HDF5 input file
PMLInside = false;
%% Creating kWaveGrid
% In all dimensions, the grid should have as small prime factors as
% possible. That's why I chose to increase the grid to the next power of 2.
% The 20 points subtracted will get re-added in the end as PML.
Nx = 2^nextpow2(grid_length/axial_spacing)-20;
Ny = 2^nextpow2(grid_width/lateral_spacing)-20;
Nz = Ny;
kgrid = kWaveGrid(Nx,axial_spacing,Ny,lateral_spacing,Nz,lateral_spacing);
kgrid.dt = time_step;
kgrid.Nt = round(total_time/time_step);
%% Creating medium
sound_speedwater = 1.40238677e3 + 5.03798765*T - 5.80980033e-2*T.^2 ...
+ 3.34296650e-4*T.^3 - 1.47936902e-6*T.^4 + 3.14893508e-9*T.^5; %[m/s]
densitywater = 999.9 + T.*(0.06079 + T.*(-0.008313 + T.*(6.595e-5 + T.* ...
(-4.437e-7 + 1.369e-9.*T))));%[kg/m^3]
medium = struct('sound_speed', sound_speedwater,...
'density', densitywater, 'sound_speed_ref', sound_speedwater);
%% Creating source
raw_source.p0 = makeBall(500,1000,1000,200,500,500,50);
% This is similar to what my function new_read_Geant does - it creates a
% 500x1000x1000 source matrix. I simplified this to a ball source instead of
% a Geant4 simulation output for you. The behaviour in terms of instablilty
% is the same.
%% Adjusting the source to the spacing of the k-Wave grid
% This part corresponds to my adjust_grids function, in an attempt to
% become independent of the inherent grid spacing of the source file.
[x_rawlength, y_rawlength, z_rawlength] = size(raw_source.p0);
x_newlength = ceil(x_rawlength * d_geant_grid / kgrid.dx);
y_newlength = ceil(y_rawlength * d_geant_grid / kgrid.dy);
z_newlength = ceil(z_rawlength * d_geant_grid / kgrid.dz);
[y_geant, x_geant, z_geant] = meshgrid((1 : y_rawlength) * d_geant_grid,...
(1 : x_rawlength) * d_geant_grid, (1 : z_rawlength) * d_geant_grid);
[y_kgrid, x_kgrid, z_kgrid] = meshgrid((1 : y_newlength) * kgrid.dy,...
(1 : x_newlength) * kgrid.dx, (1 : z_newlength) * kgrid.dz);
raw_source.p0 = interp3(y_geant, x_geant, z_geant, raw_source.p0,...
y_kgrid, x_kgrid, z_kgrid);
%% Put the source into the kgrid
% This is similar to my makeSource_spatial3D function. It creates an empty
% container of the size of the k-Wave grid, and puts the raw source into
% this container.
source.p0 = zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz);
y_shift = floor((kgrid.Ny - y_newlength) / 2);
z_shift = floor((kgrid.Nz - z_newlength) / 2);
source.p0(1 : x_newlength, y_shift + (1 : y_newlength),...
z_shift + (1 : z_newlength)) = raw_source.p0;
%% Creating sensor
y_bragg = floor(kgrid.Ny/2)+1; % index
z_bragg = floor(kgrid.Nz/2)+1; % index
[~,x_bragg] = max(source.p0(:,y_bragg,z_bragg)); %index
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(x_bragg + round(0.025/kgrid.dx), y_bragg, z_bragg) = 1;
%% Additional input parameters
input_fullpath = strcat('/project/med4/IONOAKUSTIK/SimulationAndreas/',...
SaveToDisk,'.h5');
cvsmooth = [smooth_p0, smooth_soundspeed, smooth_density];
input_args={'PMLInside', PMLInside, 'Smooth', cvsmooth, 'DataCast',...
DataCast, 'StreamToDisk', StreamToDisk, 'SaveToDisk', input_fullpath};
kspaceFirstOrder3D(kgrid,medium,source,sensor,input_args{:});