xolotl/solver/src/monitor/Monitor3D.cpp
/**
 * This is a monitoring method that will update an hdf5 file every given time.
 */
This is a monitoring method that will update an hdf5 file every given time.
This is a monitoring method that will update an hdf5 file every given time .
PetscErrorCode
startStop3D(TS ts, PetscInt timestep, PetscReal time, Vec solution, void*)
{
	PetscErrorCode ierr;
	const double ****solutionArray, *gridPointSolution;
	IdType xs, xm, Mx, ys, ym, My, zs, zm, Mz;
	PetscFunctionBeginUser;
	auto& solverHandler = PetscSolver::getSolverHandler();
	solverHandler.getLocalCoordinates(xs, xm, Mx, ys, ym, My, zs, zm, Mz);
	double previousTime = solverHandler.getPreviousTime();
	double dt = time - previousTime;
	if (((PetscInt)((time + dt / 10.0) / hdf5Stride3D) <= hdf5Previous3D) &&
		timestep > 0)
		PetscFunctionReturn(0);
	if ((PetscInt)((time + dt / 10.0) / hdf5Stride3D) > hdf5Previous3D)
		hdf5Previous3D++;
	auto xolotlComm = util::getMPIComm();
	int procId;
	MPI_Comm_rank(xolotlComm, &procId);
	DM da;
	ierr = TSGetDM(ts, &da);
	CHKERRQ(ierr);
	ierr = DMDAVecGetArrayDOFRead(da, solution, &solutionArray);
	CHKERRQ(ierr);
	auto& network = solverHandler.getNetwork();
	const auto dof = network.getDOF();
	double concArray[dof + 1][2];
	std::vector<std::vector<int>> surfaceIndices;
	for (auto i = 0; i < My; i++) {
		std::vector<int> temp;
		for (auto j = 0; j < Mz; j++) {
			temp.push_back(solverHandler.getSurfacePosition(i, j));
		}
		surfaceIndices.push_back(temp);
	}
	io::XFile checkpointFile(
		hdf5OutputName3D, xolotlComm, io::XFile::AccessMode::OpenReadWrite);
	double currentTimeStep;
	ierr = TSGetTimeStep(ts, &currentTimeStep);
	CHKERRQ(ierr);
	auto concGroup = checkpointFile.getGroup<io::XFile::ConcentrationGroup>();
	assert(concGroup);
	auto tsGroup = concGroup->addTimestepGroup(
		timestep, time, previousTime, currentTimeStep);
	if (solverHandler.moveSurface()) {
		tsGroup->writeSurface3D(
			surfaceIndices, nInterstitial3D, previousIFlux3D);
	}
	for (auto k = 0; k < Mz; k++) {
		for (auto j = 0; j < My; j++) {
			for (auto i = 0; i < Mx; i++) {
				MPI_Barrier(xolotlComm);
				int concSize = -1;
				int concId = 0;
				bool write = false;
				if (i >= xs && i < xs + xm && j >= ys && j < ys + ym &&
					k >= zs && k < zs + zm) {
					write = true;
					gridPointSolution = solutionArray[k][j][i];
					for (auto l = 0; l < dof + 1; l++) {
						if (std::fabs(gridPointSolution[l]) > 1.0e-16) {
							concSize++;
							concArray[concSize][0] = (double)l;
							concArray[concSize][1] = gridPointSolution[l];
						}
					}
					 one last time
					concSize++;
					concId = procId;
				}
				int concProc = 0;
				MPI_Allreduce(
					&concId, &concProc, 1, MPI_INT, MPI_SUM, xolotlComm);
				MPI_Bcast(&concSize, 1, MPI_INT, concProc, xolotlComm);
				if (concSize == 0)
					continue;
				tsGroup->writeConcentrationDataset(
					concSize, concArray, write, i, j, k);
			}
		}
	}
	ierr = DMDAVecRestoreArrayDOFRead(da, solution, &solutionArray);
	CHKERRQ(ierr);
	PetscFunctionReturn(0);
}
PetscErrorCode startStop3D ( TS ts , PetscInt timestep , PetscReal time , Vec solution , void * ) { PetscErrorCode ierr ; const double * * * * solutionArray , * gridPointSolution ; IdType xs , xm , Mx , ys , ym , My , zs , zm , Mz ; PetscFunctionBeginUser ; auto & solverHandler = PetscSolver :: getSolverHandler ( ) ; solverHandler . getLocalCoordinates ( xs , xm , Mx , ys , ym , My , zs , zm , Mz ) ; double previousTime = solverHandler . getPreviousTime ( ) ; double dt = time - previousTime ; if ( ( ( PetscInt ) ( ( time + dt / 10.0 ) / hdf5Stride3D ) <= hdf5Previous3D ) && timestep > 0 ) PetscFunctionReturn ( 0 ) ; if ( ( PetscInt ) ( ( time + dt / 10.0 ) / hdf5Stride3D ) > hdf5Previous3D ) hdf5Previous3D ++ ; auto xolotlComm = util :: getMPIComm ( ) ; int procId ; MPI_Comm_rank ( xolotlComm , & procId ) ; DM da ; ierr = TSGetDM ( ts , & da ) ; CHKERRQ ( ierr ) ; ierr = DMDAVecGetArrayDOFRead ( da , solution , & solutionArray ) ; CHKERRQ ( ierr ) ; auto & network = solverHandler . getNetwork ( ) ; const auto dof = network . getDOF ( ) ; double concArray [ dof + 1 ] [ 2 ] ; std :: vector < std :: vector < int > > surfaceIndices ; for ( auto i = 0 ; i < My ; i ++ ) { std :: vector < int > temp ; for ( auto j = 0 ; j < Mz ; j ++ ) { temp . push_back ( solverHandler . getSurfacePosition ( i , j ) ) ; } surfaceIndices . push_back ( temp ) ; } io :: XFile checkpointFile ( hdf5OutputName3D , xolotlComm , io :: XFile :: AccessMode :: OpenReadWrite ) ; double currentTimeStep ; ierr = TSGetTimeStep ( ts , & currentTimeStep ) ; CHKERRQ ( ierr ) ; auto concGroup = checkpointFile . getGroup < io :: XFile :: ConcentrationGroup > ( ) ; assert ( concGroup ) ; auto tsGroup = concGroup -> addTimestepGroup ( timestep , time , previousTime , currentTimeStep ) ; if ( solverHandler . moveSurface ( ) ) { tsGroup -> writeSurface3D ( surfaceIndices , nInterstitial3D , previousIFlux3D ) ; } for ( auto k = 0 ; k < Mz ; k ++ ) { for ( auto j = 0 ; j < My ; j ++ ) { for ( auto i = 0 ; i < Mx ; i ++ ) { MPI_Barrier ( xolotlComm ) ; int concSize = -1 ; int concId = 0 ; bool write = false ; if ( i >= xs && i < xs + xm && j >= ys && j < ys + ym && k >= zs && k < zs + zm ) { write = true ; gridPointSolution = solutionArray [ k ] [ j ] [ i ] ; for ( auto l = 0 ; l < dof + 1 ; l ++ ) { if ( std :: fabs ( gridPointSolution [ l ] ) > 1.0e-16 ) { concSize ++ ; concArray [ concSize ] [ 0 ] = ( double ) l ; concArray [ concSize ] [ 1 ] = gridPointSolution [ l ] ; } } concSize ++ ; concId = procId ; } int concProc = 0 ; MPI_Allreduce ( & concId , & concProc , 1 , MPI_INT , MPI_SUM , xolotlComm ) ; MPI_Bcast ( & concSize , 1 , MPI_INT , concProc , xolotlComm ) ; if ( concSize == 0 ) continue ; tsGroup -> writeConcentrationDataset ( concSize , concArray , write , i , j , k ) ; } } } ierr = DMDAVecRestoreArrayDOFRead ( da , solution , & solutionArray ) ; CHKERRQ ( ierr ) ; PetscFunctionReturn ( 0 ) ; }
This is a monitoring method that will update an hdf5 file every given time.
This is a monitoring method that will update an hdf5 file every given time .
// Initial declarations // Get the solver handler and local coordinates // Compute the dt // Don't do anything if it is not on the stride // Update the previous time // Gets the process ID (important when it is running in parallel) // Get the da from ts // Get the solutionArray // Get the network and dof // Create an array for the concentration // Get the vector of positions of the surface // Create a temporary vector // Add the temporary vector to the vector of surface indices // Open the existing HDF5 file. // Get the current time step // Add a concentration sub group // Write the surface positions in the concentration sub group // Loop on the full grid // Wait for all the processes // Size of the concentration that will be stored // To save which proc has the information // To know which process should write // If it is the locally owned part of the grid // Get the pointer to the beginning of the solution data for // this grid point // Loop on the concentrations // Increase concSize // Fill the concArray // Increase concSize one last time // Save the procId // Get which processor will send the information // Broadcast the size // Skip the grid point if the size is 0 // All processes create the dataset and fill it // Restore the solutionArray