Commit 9e7df55f authored by Gaurav Saxena's avatar Gaurav Saxena
Browse files

Parallelized Code -- output is 'slightly' different than the serial version of...

Parallelized Code -- output is 'slightly' different than the serial version of the code. Remember decomposition is 1-D in Z direction. To check try with : 1 process to see if output is the same...then with 2 processes then with 4...
parent 2bb697ff
......@@ -115,7 +115,7 @@ Microenvironment::Microenvironment()
thomas_setup_done = false;
diffusion_solver_setup_done = false;
diffusion_decay_solver = empty_diffusion_solver;
//diffusion_decay_solver = empty_diffusion_solver; //-->Gaurav Saxena commented this out (I think it has no effect really)
diffusion_decay_solver = diffusion_decay_solver__constant_coefficients_LOD_3D;
//This line commented by Gaurav Saxena
......@@ -661,10 +661,10 @@ std::vector<double>& Microenvironment::density_vector( int i, int j )
std::vector<double>& Microenvironment::density_vector( int n )
{ return (*p_density_vectors)[ n ]; }
void Microenvironment::simulate_diffusion_decay( double dt )
void Microenvironment::simulate_diffusion_decay( double dt, int mpi_Size, int mpi_Rank, int *mpi_Coords, int *mpi_Dims, MPI_Comm mpi_Cart_comm )
{
if( diffusion_decay_solver )
{ diffusion_decay_solver( *this, dt ); }
{ diffusion_decay_solver( *this, dt, mpi_Size, mpi_Rank, mpi_Coords, mpi_Dims, mpi_Cart_comm ); }
else
{
std::cout << "Warning: diffusion-reaction-source/sink solver not set for Microenvironment object at " << this << ". Nothing happened!" << std::endl;
......@@ -677,7 +677,9 @@ void Microenvironment::simulate_diffusion_decay( double dt )
void Microenvironment::auto_choose_diffusion_decay_solver( void )
{
// set the safest choice
diffusion_decay_solver = diffusion_decay_solver__constant_coefficients_explicit;
//diffusion_decay_solver = diffusion_decay_solver__constant_coefficients_explicit; //-->Gaurav Saxena commented this out as this was causing compilation problems
//because diffusion_decay_solver is a function pointer to a parallel version of 3D solver with constant coefficients ---> later on when we implement 'explicit' solver
//we can uncomment this
std::cout << "Warning: auto-selection of diffusion-decay-source/sink solver not fully implemented!" << std::endl;
......
......@@ -153,7 +153,7 @@ class Microenvironment
Microenvironment();
Microenvironment(std::string name);
void (*diffusion_decay_solver)( Microenvironment&, double);
void (*diffusion_decay_solver)( Microenvironment&, double, int, int, int *, int *, MPI_Comm); //-->Gaurav Saxena changed prototype
void (*bulk_supply_rate_function)( Microenvironment* pMicroenvironment, int voxel_index, std::vector<double>* write_destination );
void (*bulk_supply_target_densities_function)( Microenvironment* pMicroenvironment, int voxel_index, std::vector<double>* write_destination );
void (*bulk_uptake_rate_function)( Microenvironment* pMicroenvironment, int voxel_index, std::vector<double>* write_destination );
......@@ -224,7 +224,7 @@ class Microenvironment
std::vector<double>& density_vector( int n );
/*! advance the diffusion-decay solver by dt time */
void simulate_diffusion_decay( double dt );
void simulate_diffusion_decay( double dt, int mpi_Size, int mpi_Rank, int *mpi_Coords, int *mpi_Dims, MPI_Comm mpi_Cart_comm );//-->Gaurav Saxena changed this prototype
/*! advance the source/sink solver by dt time */
void simulate_bulk_sources_and_sinks( double dt );
......@@ -248,7 +248,7 @@ class Microenvironment
friend void diffusion_decay_solver__constant_coefficients_explicit( Microenvironment& S, double dt );
friend void diffusion_decay_solver__constant_coefficients_explicit_uniform_mesh( Microenvironment& S, double dt );
friend void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& S, double dt );
friend void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& S, double dt, int mpi_Size, int mpi_Rank, int *mpi_Coords, int *mpi_Dims, MPI_Comm mpi_Cart_comm ); //-->Gaurav Saxena changed this prototype
friend void diffusion_decay_solver__constant_coefficients_LOD_2D( Microenvironment& S, double dt );
friend void diffusion_decay_explicit_uniform_rates( Microenvironment& M, double dt );
......@@ -270,7 +270,7 @@ extern void diffusion_decay_solver__variable_coefficients_explicit( Microenviron
extern void diffusion_decay_solver__variable_coefficients_explicit_uniform_mesh( Microenvironment& S, double dt );
extern void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& S, double dt );
extern void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& S, double dt, int mpi_Size, int mpi_Rank, int *mpi_Coords, int *mpi_Dims, int *mpi_Cart_comm ); //-->Gaurav Saxena changed this prototype
extern void diffusion_decay_solver__constant_coefficients_LOD_2D( Microenvironment& S, double dt );
extern void diffusion_decay_solver__variable_coefficients_LOD_3D( Microenvironment& S, double dt );
......
......@@ -960,11 +960,11 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
std::vector<double> temp_cz(M.thomas_cz[0].size());
MPI_Irecv(&temp_cz[0], temp_cz[0].size(), MPI_DOUBLE, ser_ctr, 1111, mpi_Cart_comm, &recv_req);
MPI_Irecv(&temp_cz[0], temp_cz.size(), MPI_DOUBLE, ser_ctr, 1111, mpi_Cart_comm, &recv_req);
MPI_Wait(&recv_req, MPI_STATUS_IGNORE);
axpy( &M.thomas_denomz[0] , M.thomas_constant1 , &temp_cz[0] ); //CHECK IF &temp_cz[0] is OK, axpy() in BioFVM_vector.cpp
axpy( &M.thomas_denomz[0] , M.thomas_constant1 , temp_cz ); //CHECK IF &temp_cz[0] is OK, axpy() in BioFVM_vector.cpp
M.thomas_cz[0] /= M.thomas_denomz[0]; // the value at size-1 is not actually used
}
......@@ -1089,36 +1089,36 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
/* Then M.apply_dirichlet_conditions(); */
/*--------------------------------------------------------------------------------*/
#pragma omp parallel for
for( int j=0; j < M.mesh.y_coordinates.size() ; j++ )
{
for( int i=0; i < M.mesh.x_coordinates.size() ; i++ )
{
// remaining part of forward elimination, using pre-computed quantities
int n = M.voxel_index(i,j,0);
(*M.p_density_vectors)[n] /= M.thomas_denomz[0];
// should be an empty loop if mesh.z_coordinates.size() < 2
for( int k=1; k < M.mesh.z_coordinates.size() ; k++ )
{
n = M.voxel_index(i,j,k);
axpy( &(*M.p_density_vectors)[n] , M.thomas_constant1 , (*M.p_density_vectors)[n-M.thomas_k_jump] );
(*M.p_density_vectors)[n] /= M.thomas_denomz[k];
}
// for parallelization need to break forward elimination and back substitution into
// should be an empty loop if mesh.z_coordinates.size() < 2
for( int k = M.mesh.z_coordinates.size()-2 ; k >= 0 ; k-- )
{
n = M.voxel_index(i,j,k);
naxpy( &(*M.p_density_vectors)[n] , M.thomas_cz[k] , (*M.p_density_vectors)[n+M.thomas_k_jump] );
// n -= i_jump;
}
}
}
// #pragma omp parallel for
// for( int j=0; j < M.mesh.y_coordinates.size() ; j++ )
// {
//
// for( int i=0; i < M.mesh.x_coordinates.size() ; i++ )
// {
//
// // remaining part of forward elimination, using pre-computed quantities
// int n = M.voxel_index(i,j,0);
// (*M.p_density_vectors)[n] /= M.thomas_denomz[0];
//
// // should be an empty loop if mesh.z_coordinates.size() < 2
// for( int k=1; k < M.mesh.z_coordinates.size() ; k++ )
// {
// n = M.voxel_index(i,j,k);
// axpy( &(*M.p_density_vectors)[n] , M.thomas_constant1 , (*M.p_density_vectors)[n-M.thomas_k_jump] );
// (*M.p_density_vectors)[n] /= M.thomas_denomz[k];
// }
//
// // for parallelization need to break forward elimination and back substitution into
// // should be an empty loop if mesh.z_coordinates.size() < 2
//
// for( int k = M.mesh.z_coordinates.size()-2 ; k >= 0 ; k-- )
// {
// n = M.voxel_index(i,j,k);
// naxpy( &(*M.p_density_vectors)[n] , M.thomas_cz[k] , (*M.p_density_vectors)[n+M.thomas_k_jump] );
// // n -= i_jump;
// }
// }
// }
/*----------------------------TILL HERE----------------------------------------------*/
......@@ -1127,8 +1127,8 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
/*-----------------------------------------------------------------------------------*/
int n1;
int x_size = M.Mesh.x_coordinates.size();
int y_size = M.Mesh.y_coordinates.size();
int x_size = M.mesh.x_coordinates.size();
int y_size = M.mesh.y_coordinates.size();
int p_size = (*M.p_density_vectors)[0].size(); //All p_density_vectors elements have same size, use anyone
int snd_data_size = x_size * y_size * p_size;
......@@ -1155,9 +1155,9 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
if(rank == 0)
{
#pragma omp parallel for
for(j=0; j<= M.mesh.y_coordinates()-1; j++)
for(int j=0; j<= M.mesh.y_coordinates.size()-1; j++)
{
for(i=0; i<= M.mesh.x_coordinates.size()-1; i++)
for(int i=0; i<= M.mesh.x_coordinates.size()-1; i++)
{
int n = M.voxel_index(i,j,0);
(*M.p_density_vectors)[n] /= M.thomas_denomz[0];
......@@ -1165,7 +1165,7 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
for( int k=1; k < M.mesh.z_coordinates.size() ; k++ )
{
n = M.voxel_index(i,j,k);
n1 = M.voxel_indexi,j,k-1);
n1 = M.voxel_index(i,j,k-1);
axpy( &(*M.p_density_vectors)[n] , M.thomas_constant1 , (*M.p_density_vectors)[n1] );
(*M.p_density_vectors)[n] /= M.thomas_denomz[k];
}
......@@ -1175,9 +1175,9 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
else
{
#pragma omp parallel for
for(j=0; j<= M.mesh.y_coordinates()-1; j++)
for(int j=0; j<= M.mesh.y_coordinates.size()-1; j++)
{
for(i=0; i<= M.mesh.x_coordinates.size()-1; i++)
for(int i=0; i<= M.mesh.x_coordinates.size()-1; i++)
{
int n = M.voxel_index(i,j,0); //Need to consider case separately for k=0, as k-1 would be -1 !
axpy( &(*M.p_density_vectors)[n] , M.thomas_constant1 , block3d[j][i]);
......@@ -1186,39 +1186,38 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
for( int k=1; k < M.mesh.z_coordinates.size() ; k++ )
{
n = M.voxel_index(i,j,k);
n1 = M.voxel_indexi,j,k-1);
n1 = M.voxel_index(i,j,k-1);
axpy( &(*M.p_density_vectors)[n] , M.thomas_constant1 , (*M.p_density_vectors)[n1] );
(*M.p_density_vectors)[n] /= M.thomas_denomz[k];
}
}
}
}
if(rank < (size-1))
{
int z_end = M.mesh.z_coordinates.size()-1;
int ctr=0;
for(j=0; j<= M.mesh.y_coordinates()-1; j++)
for(int j=0; j<= M.mesh.y_coordinates.size()-1; j++)
{
for(i=0; i<= M.mesh.x_coordinates.size()-1; i++)
for(int i=0; i<= M.mesh.x_coordinates.size()-1; i++)
{
int n = M.voxel_index(i,j,z_end);
for(int ele=0; ele <= (*M.p_density_vectors)[n].size()-1; ele++)
snd_data[ctr++] = (*M.p_density_vectors)[n][ele];
}
}
MPI_Isend(snd_data, snd_data_size, MPI_DOUBLE, ser_ctr+1, 1111, mpi_Cart_comm, &send_req);
MPI_Isend(&snd_data[0], snd_data_size, MPI_DOUBLE, ser_ctr+1, 1111, mpi_Cart_comm, &send_req);
MPI_Wait(&send_req, MPI_STATUS_IGNORE);
}
}
if(rank == (ser_ctr+1) && rank <= (size-1))
{
//Receive the data here and try to put in same format as vector of vectors
MPI_Irecv(rcv_data, rcv_data_size, MPI_DOUBLE, ser_ctr, 1111, mpi_Cart_comm, &recv_req);
MPI_Irecv(&rcv_data[0], rcv_data_size, MPI_DOUBLE, ser_ctr, 1111, mpi_Cart_comm, &recv_req);
MPI_Wait(&recv_req,MPI_STATUS_IGNORE);
int ctr = 0;
for(int m=0; m<y_size; j++)
for(int m=0; m<y_size; m++)
{
for(int n=0; n<x_size; n++)
{
......@@ -1237,22 +1236,22 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
/* CODE FOR BACK SUBSITUTION */
/*-----------------------------------------------------------------------------------*/
int n2;
for(int ser_ctr=size-1; ser_ctr >= 0; ser_ctr++)
for(int ser_ctr=size-1; ser_ctr >= 0; ser_ctr--)
{
if(rank == ser_ctr)
{
if(rank == (size-1))
{
#pragma omp parallel for
for(j=0; j<= M.mesh.y_coordinates()-1; j++)
for(int j=0; j<= M.mesh.y_coordinates.size()-1; j++)
{
for(i=0; i<= M.mesh.x_coordinates.size()-1; i++)
for(int i=0; i<= M.mesh.x_coordinates.size()-1; i++)
{
for( int k=M.mesh.z_coordinates.size()-2; k>=0 ; k-- )
{
n = M.voxel_index(i,j,k);
n2 = M.voxel_index(i,j,k+1);
int n = M.voxel_index(i,j,k);
int n2 = M.voxel_index(i,j,k+1);
naxpy( &(*M.p_density_vectors)[n] , M.thomas_cz[k] , (*M.p_density_vectors)[n2] );
}
}
......@@ -1261,48 +1260,47 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
else
{
#pragma omp parallel for
for(j=0; j<= M.mesh.y_coordinates()-1; j++)
for(int j=0; j<= M.mesh.y_coordinates.size()-1; j++)
{
for(i=0; i<= M.mesh.x_coordinates.size()-1; i++)
for(int i=0; i<= M.mesh.x_coordinates.size()-1; i++)
{
int n = M.voxel_index(i,j,M.mesh.z_coordinates.size()-1);
naxpy( &(*M.p_density_vectors)[n] , M.thomas_cz[k] , block3d[j][i]);
naxpy( &(*M.p_density_vectors)[n] , M.thomas_cz[M.mesh.z_coordinates.size()-1] , block3d[j][i]);
for( int k=M.mesh.z_coordinates.size()-2; k >=0; k-- )
{
n = M.voxel_index(i,j,k);
n2 = M.voxel_index(i,j,k+1);
int n = M.voxel_index(i,j,k);
int n2 = M.voxel_index(i,j,k+1);
naxpy( &(*M.p_density_vectors)[n] , M.thomas_cz[k] , (*M.p_density_vectors)[n2] );
}
}
}
}
if(rank > 0)
{
int z_start = 0;
int ctr=0;
for(j=0; j<= M.mesh.y_coordinates()-1; j++)
for(int j=0; j<= M.mesh.y_coordinates.size()-1; j++)
{
for(i=0; i<= M.mesh.x_coordinates.size()-1; i++)
for(int i=0; i<= M.mesh.x_coordinates.size()-1; i++)
{
int n = M.voxel_index(i,j,z_start);
for(int ele=0; ele <= (*M.p_density_vectors)[n].size()-1; ele++)
snd_data[ctr++] = (*M.p_density_vectors)[n][ele];
}
}
MPI_Isend(snd_data, snd_data_size, MPI_DOUBLE, ser_ctr-1, 1111, mpi_Cart_comm, &send_req);
MPI_Isend(&snd_data[0], snd_data_size, MPI_DOUBLE, ser_ctr-1, 1111, mpi_Cart_comm, &send_req);
MPI_Wait(&send_req, MPI_STATUS_IGNORE);
}
}
if(rank == (ser_ctr-1) && rank >= 0)
{
//Receive the data here and try to put in same format as vector of vectors
MPI_Irecv(rcv_data, rcv_data_size, MPI_DOUBLE, ser_ctr, 1111, mpi_Cart_comm, &recv_req);
MPI_Irecv(&rcv_data[0], rcv_data_size, MPI_DOUBLE, ser_ctr, 1111, mpi_Cart_comm, &recv_req);
MPI_Wait(&recv_req,MPI_STATUS_IGNORE);
int ctr = 0;
for(int m=0; m<y_size; j++)
for(int m=0; m<y_size; m++)
{
for(int n=0; n<x_size; n++)
{
......@@ -1316,8 +1314,6 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
MPI_Barrier(mpi_Cart_comm);
}
M.apply_dirichlet_conditions();
// reset gradient vectors
......@@ -1326,6 +1322,13 @@ void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M,
return;
}
//-->Had to create a function like this else compiler complains of undefined reference to this function due to call in initialize_microenvironment() in microenvironment.cpp
//-->Remember I have commented out LOD_2D and LOD_3D in solvers.cpp
void diffusion_decay_solver__constant_coefficients_LOD_2D( Microenvironment& M, double dt, int mpi_Size, int mpi_Rank, int *mpi_Coords, int *mpi_Dims, MPI_Comm mpi_Cart_comm )
{
}
};
......
This diff is collapsed.
......@@ -55,9 +55,9 @@ namespace BioFVM{
// /*! diffusion-decay solvers for the equation du/dt = D*Laplacian(u) - lambda*u - U(x)*u + M(X)*(uT-u) */
// /*! diffusion-decay solver: 3D LOD implicit (stable method). D and r uniform */
void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M, double dt ); // done
void diffusion_decay_solver__constant_coefficients_LOD_3D( Microenvironment& M, double dt, int mpi_Size, int mpi_Rank, int *mpi_Coords, int *mpi_Dims, MPI_Comm mpi_Cart_comm );//-->Gaurav Saxena changed this prototype // done
// /*! diffusion-decay solver: 2D LOD implicit (stable method). D and r uniform */
void diffusion_decay_solver__constant_coefficients_LOD_2D( Microenvironment& M, double dt ); // done
void diffusion_decay_solver__constant_coefficients_LOD_2D( Microenvironment& M, double dt, int mpi_Size, int mpi_Rank, int *mpi_Coords, int *mpi_Dims, MPI_Comm mpi_Cart_comm );//-->Gaurav Saxena changed this function prototype to avoid compilation problems (not yet implemented parallel version) // done
/*! This solves for constant diffusion coefficients on a general mesh using the
explicit stepping for the diffusion operator, and implicit stepping for all
......@@ -70,4 +70,4 @@ void diffusion_decay_solver__constant_coefficients_explicit( Microenvironment& M
void diffusion_decay_solver__constant_coefficients_explicit_uniform_mesh( Microenvironment& M, double dt );
};
#endif
\ No newline at end of file
#endif
No preview for this file type
......@@ -446,13 +446,14 @@ int main( int argc, char* argv[] )
/*------------------------------------------------------------*/
/* MPI Comm Size and Rank */
/* Right now it should be 1-D Z decomposition */
/*------------------------------------------------------------*/
mpi_Dims[0] = 0;
mpi_Dims[1] = 0;
mpi_Dims[2] = 0;
MPI_Dims_create(mpi_Size,3,mpi_Dims);
mpi_Dims[0] = 1;
mpi_Dims[1] = 1;
mpi_Dims[2] = 4;
//MPI_Dims_create(mpi_Size,3,mpi_Dims);
mpi_Is_periodic[0]=0;
mpi_Is_periodic[1]=0;
......@@ -498,7 +499,7 @@ int main( int argc, char* argv[] )
/* Resize local Microenvironment on each process */
/*------------------------------------------------------------*/
microenvironment.resize_space_uniform( minX,maxX,minY,maxY,minZ,maxZ, mesh_resolution,mpi_Dims, mpi_Coords);
microenvironment.resize_space_uniform( minX,maxX,minY,maxY,minZ,maxZ, mesh_resolution, mpi_Dims, mpi_Coords);
microenvironment.display_information( std::cout );
......@@ -564,7 +565,7 @@ int main( int argc, char* argv[] )
(*temp_point_sink->uptake_rates)[0]=0.8;
temp_point_sink->set_internal_uptake_constants(dt);
}
*/
double t = 0.0;
double t_max=5;
......@@ -572,12 +573,12 @@ int main( int argc, char* argv[] )
while( t < t_max )
{
microenvironment.simulate_cell_sources_and_sinks( dt );
microenvironment.simulate_diffusion_decay( dt );
microenvironment.simulate_diffusion_decay( dt, mpi_Size, mpi_Rank, mpi_Coords, mpi_Dims, mpi_Cart_comm );
t += dt;
}
microenvironment.write_to_matlab( "final.mat" );
microenvironment.write_to_matlab( "final.mat", mpi_Rank, mpi_Size, mpi_Cart_comm ); //Remember to use Parallel Version !
std::cout<<"done!"<<std::endl;
*/
/*-------------------------------------------------------------------------*/
/* MPI Finalize */
......
File added
This diff is collapsed.
#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=8
#SBATCH --ntasks-per-node=4
#SBATCH --cpus-per-task=1
#SBATCH -t 00:05:00
#SBATCH -o output-%j
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment