k-Wave

1. samantp
   

   Member
   Joined: Jun '16
   Posts: 17

   Hello,
   I want to run 3D simulations using extremely high grid sizes (2^11 points in the z direction representing a 10 micron physical distance). I have a locally compiled C++ script that I am running on a supercomputer to this end and I wanted to model the 3D propagation of a time varying Gaussian source and its attenuation through space as time progresses.

   Now, in order to avoid overloading MATLAB's memory in the generation of an input .h5 file, I have been trying to keep the number of gridpoints in the x and y directions to a minimum. I was wondering, does this have any effect on the output signal in the z direction? I essentially want to generate a signal at z=0 and measure it as it propagates all the way to z=10um, so in principle I don't really need the x and y spaces. However using the 1D or 2D model corresponds to a plane or line source in 3D, and so I am forced to use 3D simulations and have large computation times. I have been trying to keep the number of gridpoints in x and y to a minimum, but have noticed that the wave from my point source keeps bouncing over to the other side of the simulation space (despite my initial source being a few grid points from the PML).

   My question is, does it matter what the number of grid points in x and y are if I only care about 3D propagation in the z direction? Is there any benefit in setting Nx or Ny to be anything higher than the lowest possible power of 2 that is still higher than PMLsize+1 or does it not matter?

   Thank you!!

   Posted 6 years ago #

2. samantp
   

   Member
   Joined: Jun '16
   Posts: 17

   An additional question,
   I am running a simulation with the following parameters, essentially a gaussian time varying source resembling a picosecond laser pulse for photoacoustic simulations, my parameters are below

   % create the computational grid
   Nx = 22;            % number of grid points in the x direction
   Ny = 22;            % number of grid points in the y direction
   Nz = 70;            % number of grid points in the z direction
   dx = 1e-14;        % grid point spacing in the x direction [m]
   dy = 1e-14;        % grid point spacing in the y direction [m]
   dz = 1e-14;        % grid point spacing in the z direction [m]
   kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
   % define the properties of the propagation medium
   medium.sound_speed = 1480*ones(Nx, Ny, Nz);	% [m/s]
   %medium.sound_speed(1:Nx/2, :, :) = 1800;    % [m/s]
   medium.density = 10000*ones(Nx, Ny, Nz);     % [kg/m^3]
   %medium.density(:, Ny/4:end, :) = 1200;      % [kg/m^3]
   medium.alpha_coeff = 0.2;  % [dB/(MHz^y cm)]
   medium.alpha_power = 2;
   %dt=0.01e-12;
   [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
   PhaseShift=7e-12; %[s] Phase shift in time of pressure pulse
   % create the time array
   
   % define two point sources at (Nx/2, Ny/2, Nz/2) and Nx/2, Ny/2kgrid.t_array), Nz/2+5
   %source_radius = 5;  % [grid points]
   source.p_mask = zeros(Nx, Ny, Nz);
   %source.p_mask(round(Nx/2), round(Ny/2), 50) = 1;
   source.p_mask(1, 1, 10) = 1;
   
   %source.p_mask(Nx/2, Ny/2, Nz/2+5) = 1;
   % define a time varying source
   %source_freq = 2e11;  % [Hz]
   source_mag = 2500;     % [Pa] Max Amplitude
   Pulsewidth=7e-12; %Pulse Width [s]
   C=source_mag*Pulsewidth^2/(exp(-0.5)); %Constant adjusted so that peaks of pressure are equal to specified Source magnitude
   
   source.p = (-C/Pulsewidth^3)*(kgrid.t_array-PhaseShift).*exp(-((kgrid.t_array-PhaseShift).^2)/(2*Pulsewidth^2)); %Derivation result as listed in .pdf
   % filter the source to remove high frequencies not supported by the grid
   source.p = filterTimeSeries(kgrid, medium, source.p);

   Which is adapted code from one of the examples (the time varying monopole). However I keep getting the warning "WARNING: time step may be too large for a stable simulation", even though the timestep is decided by MakeTime. What is the general procedure for determining an appropriate timestep if not MakeTime? I have also tried manually making the timestep as small as
   

   dt=0.001e-12;
   kgrid.t_array=[0:dt:10000*dt];

   but I still get the same warning and the simulation blows up. Any idea why this might be happening?

   Thanks again!

   Posted 6 years ago #

3. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Samantp,

   Your approach sounds sensible, and I can't immediately think of any problems, provided your grid spacing is equal (dx = dy = dz). One possible thing to watch out for - the PML doesn't work very well below 3 points per wavelength (see here). So you could try reducing the bandwidth of your source to at least 3 (possibly 4) points per wavelength. Another thing to try is increasing the PML size to 20 (leaving the PML alpha at its default value).

   The makeTime method doesn't consider the stability of the simulation, it just sets the time step based on the CFL. In an absorbing medium, the simulation is no longer unconditionally stable, and you can find the largest stable time step using the function checkStability. This is the function used by k-Wave to print the warning to the screen.

   Hope that helps,

   Brad.

   Posted 6 years ago #

4. samantp
   

   Member
   Joined: Jun '16
   Posts: 17

   Thank you Dr. Treeby! One last question, is it possible to have assymetrical time-steps in simulations? For example my simulation consists of a very short laser pulse (7ps) followed by the propagation of the acoustic wave, which occurs on much larger timescales. Is it possible then to run (for example) the first 20ps of the simulation using a very short timestep and the next 100ns at a different timestep?

   Thanks

   Posted 6 years ago #

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   In principle, yes, in practice, no. The assumption of a fixed time step is more or less built into the fabric of k-Wave at the moment.

   Brad.

   Posted 6 years ago #

6. samantp
   

   Member
   Joined: Jun '16
   Posts: 17

   Thank you Dr Treeby!

   Posted 6 years ago #

RSS feed for this topic

Reply

You must log in to post.