k-Wave

1. fantasyhaha
   

   Member
   Joined: Sep '16
   Posts: 9

   Hello developers,
   We meet a problem in simulation in a relatively large space.
   Parameters are as follows:
   PML_size = 10; % size of the PML in grid points
   Nx = 420; % number of grid points in the x direction
   Ny = Nx; % number of grid points in the y direction
   Nz = Nx; % number of grid points in the z direction
   dx = 0.2e-3; % grid point spacing in the x direction [m]
   dy = dx; % grid point spacing in the y direction [m]
   dz = dx; % grid point spacing in the z direction [m]
   kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
   % define the properties of the propagation medium
   medium.sound_speed = 1500; % [m/s]
   sz=128;
   p0_init=zeros(Nx,Ny,Nz);
   p0_init((Nx/2-sz/2+1):(Nx/2+sz/2+1),(Nx/2-sz/2+1):(Nx/2+sz/2+1),(Nx/2-sz/2+1):(Nx/2+sz/2+1))=10;
   % smooth the initial pressure distribution and restore the magnitude
   p0 = smooth(p0_init, true);

   % assign to the source structure
   source.p0 = p0;

   % define a Cartesian spherical sensor
   sensor_radius = 40e-3; % [m]
   center_pos = [0, 0, 0]; % [m]
   num_sensor_points = 1024;
   sensor_mask = makeCartSphere(sensor_radius, num_sensor_points, center_pos, true);

   % assign to the sensor structure
   sensor.mask = sensor_mask;

   % create the time array
   cfl=0.15;
   % kgrid.makeTime(medium.sound_speed);
   kgrid.makeTime(medium.sound_speed);
   % kgrid.t_array(3470:end)=[];
   % set the input arguements
   input_args = {'PMLSize', PML_size, 'PMLInside', false, 'PlotPML', true, ...
   'Smooth', false, 'DataCast', 'single', 'CartInterp', 'nearest'};

   % run the simulation
   sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

   We meet some problems during the simulation and the signal then the 'sensor_data' becomes continuous oscillation with very high amplitude.(I hope I can share you a figure to show the results) This happened before the wave reach our transducer or PML. We have already tested different 'cfl' values but can not solve this problem.

   Posted 12 months ago #

2. fantasyhaha
   

   Member
   Joined: Sep '16
   Posts: 9

   This problem seems to happen on the server with a Linux system(Ubuntu 22.04). When I test it on a high-performance PC with Windows 10 with the same Matlab version, it's stable.

   Posted 12 months ago #

3. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Hi fantasyhaha,

   That sounds strange. You run the same simulation script on the same version of matlab and it's different if you're running Linux or Windows? Please can you send more details about the differences in the output.

   Best wishes
   Ben

   Posted 11 months ago #

4. fantasyhaha
   

   Member
   Joined: Sep '16
   Posts: 9

   Thank you sir.
   I have checked the code and I can show you the test code as follow:
   clearvars;
   close all
   % =========================================================================
   % SIMULATION
   % =========================================================================

   % create the computational grid
   PML_size = 20; % size of the PML in grid points
   Nx = 446; % number of grid points in the x direction
   Ny = Nx; % number of grid points in the y direction
   Nz = Nx; % number of grid points in the z direction
   dx = 0.2e-3; % grid point spacing in the x direction [m]
   dy = dx; % grid point spacing in the y direction [m]
   dz = dx; % grid point spacing in the z direction [m]
   kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);

   % define the properties of the propagation medium
   medium.sound_speed = 1500; % [m/s]

   % create initial pressure distribution using makeBall
   % ball_magnitude = 10; % [Pa]
   ball_x_pos = Nx/2; % [grid points]
   ball_y_pos = Ny/2; % [grid points]
   ball_z_pos = Nz/2; % [grid points]
   % ball_radius = 30; % [grid points]
   % p0_binary = ball_magnitude * makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);
   sz=80;
   p0_init=zeros(Nx,Ny,Nz);
   p0_init((Nx/2-sz/2+1):(Nx/2+sz/2+1),(Nx/2-sz/2+1):(Nx/2+sz/2+1),(Nx/2-sz/2+1):(Nx/2+sz/2+1))=10;
   load('brain_up_2.mat')
   brain=imresize3(V,[Nx Ny Nz]);
   brain_BW=imbinarize(brain);
   p0_init=p0_init.*brain_BW;
   % smooth the initial pressure distribution and restore the magnitude
   p0 = smooth(p0_init, true);

   % assign to the source structure
   source.p0 = p0;

   % define a Cartesian spherical sensor
   sensor_radius = 40e-3; % [m]
   center_pos = [0.5, 0.5, 0.5]*Nx*dx; % [m]
   num_sensor_points = 1024;
   sensor_mask = makeCartSphere(sensor_radius, num_sensor_points, center_pos, true);

   % assign to the sensor structure
   sensor.mask = sensor_mask;

   % create the time array
   cfl=0.3;
   % kgrid.makeTime(medium.sound_speed);
   kgrid.makeTime(medium.sound_speed,cfl);
   time_limit=(Nx+sz)*dx/medium.sound_speed;
   kgrid.t_array(kgrid.t_array>time_limit)=[];
   % set the input arguements
   input_args = {'PMLSize', PML_size, 'PMLInside', false, 'PlotPML', true, ...
   'Smooth', false, 'DataCast', 'single', 'CartInterp', 'nearest'};
   % % plot the initial pressure and sensor surface in voxel form
   % voxelPlot(double(p0_init | cart2grid(kgrid, sensor_mask)));
   % view([60, 20]);
   dt_stability_limit = checkStability(kgrid, medium);
   % run the simulation
   sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
   %%

   Posted 11 months ago #

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,588

   Thanks fantasyhaha. Can you send more details about the differences in the output between the linux and windows simulations?

   Posted 10 months ago #

RSS feed for this topic

Reply

You must log in to post.