xolotl/solver/src/monitor/Monitor3D.cpp
/**
 * This is a monitoring method that will compute the helium retention.
 */
This is a monitoring method that will compute the helium retention.
This is a monitoring method that will compute the helium retention .
PetscErrorCode
computeHeliumRetention3D(TS ts, PetscInt, PetscReal time, Vec solution, void*)
{
	PetscErrorCode ierr;
	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);
	auto fluxHandler = solverHandler.getFluxHandler();
	DM da;
	ierr = TSGetDM(ts, &da);
	CHKERRQ(ierr);
	auto grid = solverHandler.getXGrid();
	using NetworkType = core::network::IPSIReactionNetwork;
	auto& network = dynamic_cast<NetworkType&>(solverHandler.getNetwork());
	const auto dof = network.getDOF();
	double hy = solverHandler.getStepSizeY();
	double hz = solverHandler.getStepSizeZ();
	double ****solutionArray, *gridPointSolution;
	ierr = DMDAVecGetArrayDOFRead(da, solution, &solutionArray);
	CHKERRQ(ierr);
	auto numSpecies = network.getSpeciesListSize();
	auto myConcData = std::vector<double>(numSpecies, 0.0);
	for (auto zk = zs; zk < zs + zm; zk++) {
		for (auto yj = ys; yj < ys + ym; yj++) {
			auto surfacePos = solverHandler.getSurfacePosition(yj, zk);
			for (auto xi = xs; xi < xs + xm; xi++) {
				if (xi < surfacePos + solverHandler.getLeftOffset() ||
					xi >= Mx - solverHandler.getRightOffset())
					continue;
				gridPointSolution = solutionArray[zk][yj][xi];
				double hx = grid[xi + 1] - grid[xi];
				using HostUnmanaged = Kokkos::View<double*, Kokkos::HostSpace,
					Kokkos::MemoryUnmanaged>;
				auto hConcs = HostUnmanaged(gridPointSolution, dof);
				auto dConcs = Kokkos::View<double*>("Concentrations", dof);
				deep_copy(dConcs, hConcs);
				for (auto id = core::network::SpeciesId(numSpecies); id; ++id) {
					myConcData[id()] +=
						network.getTotalAtomConcentration(dConcs, id, 1) * hx *
						hy * hz;
				}
			}
		}
	}
	auto xolotlComm = util::getMPIComm();
	int procId;
	MPI_Comm_rank(xolotlComm, &procId);
	auto totalConcData = std::vector<double>(numSpecies, 0.0);
	MPI_Reduce(myConcData.data(), totalConcData.data(), myConcData.size(),
		MPI_DOUBLE, MPI_SUM, 0, xolotlComm);
	if (procId == 0) {
		PetscInt Mx, My, Mz;
		ierr = DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, &Mz, PETSC_IGNORE,
			PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE,
			PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);
		CHKERRQ(ierr);
		double surface = (double)(My * Mz) * hy * hz;
		for (auto i = 0; i < numSpecies; ++i) {
			totalConcData[i] /= surface;
		}
		double fluence = fluxHandler->getFluence();
		std::cout << "\nTime: " << time << std::endl;
		for (auto id = core::network::SpeciesId(numSpecies); id; ++id) {
			std::cout << network.getSpeciesName(id)
					  << " content = " << totalConcData[id()] << '\n';
		}
		std::cout << "Fluence = " << fluence << "\n" << std::endl;
		std::ofstream outputFile;
		outputFile.open("retentionOut.txt", std::ios::app);
		outputFile << fluence << " ";
		for (auto i = 0; i < numSpecies; ++i) {
			outputFile << totalConcData[i] << " ";
		}
		outputFile << std::endl;
		outputFile.close();
	}
	ierr = DMDAVecRestoreArrayDOFRead(da, solution, &solutionArray);
	CHKERRQ(ierr);
	PetscFunctionReturn(0);
}
PetscErrorCode computeHeliumRetention3D ( TS ts , PetscInt , PetscReal time , Vec solution , void * ) { PetscErrorCode ierr ; 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 ) ; auto fluxHandler = solverHandler . getFluxHandler ( ) ; DM da ; ierr = TSGetDM ( ts , & da ) ; CHKERRQ ( ierr ) ; auto grid = solverHandler . getXGrid ( ) ; using NetworkType = core :: network :: IPSIReactionNetwork ; auto & network = dynamic_cast < NetworkType & > ( solverHandler . getNetwork ( ) ) ; const auto dof = network . getDOF ( ) ; double hy = solverHandler . getStepSizeY ( ) ; double hz = solverHandler . getStepSizeZ ( ) ; double * * * * solutionArray , * gridPointSolution ; ierr = DMDAVecGetArrayDOFRead ( da , solution , & solutionArray ) ; CHKERRQ ( ierr ) ; auto numSpecies = network . getSpeciesListSize ( ) ; auto myConcData = std :: vector < double > ( numSpecies , 0.0 ) ; for ( auto zk = zs ; zk < zs + zm ; zk ++ ) { for ( auto yj = ys ; yj < ys + ym ; yj ++ ) { auto surfacePos = solverHandler . getSurfacePosition ( yj , zk ) ; for ( auto xi = xs ; xi < xs + xm ; xi ++ ) { if ( xi < surfacePos + solverHandler . getLeftOffset ( ) || xi >= Mx - solverHandler . getRightOffset ( ) ) continue ; gridPointSolution = solutionArray [ zk ] [ yj ] [ xi ] ; double hx = grid [ xi + 1 ] - grid [ xi ] ; using HostUnmanaged = Kokkos :: View < double * , Kokkos :: HostSpace , Kokkos :: MemoryUnmanaged > ; auto hConcs = HostUnmanaged ( gridPointSolution , dof ) ; auto dConcs = Kokkos :: View < double * > ( " " , dof ) ; deep_copy ( dConcs , hConcs ) ; for ( auto id = core :: network :: SpeciesId ( numSpecies ) ; id ; ++ id ) { myConcData [ id ( ) ] += network . getTotalAtomConcentration ( dConcs , id , 1 ) * hx * hy * hz ; } } } } auto xolotlComm = util :: getMPIComm ( ) ; int procId ; MPI_Comm_rank ( xolotlComm , & procId ) ; auto totalConcData = std :: vector < double > ( numSpecies , 0.0 ) ; MPI_Reduce ( myConcData . data ( ) , totalConcData . data ( ) , myConcData . size ( ) , MPI_DOUBLE , MPI_SUM , 0 , xolotlComm ) ; if ( procId == 0 ) { PetscInt Mx , My , Mz ; ierr = DMDAGetInfo ( da , PETSC_IGNORE , & Mx , & My , & Mz , PETSC_IGNORE , PETSC_IGNORE , PETSC_IGNORE , PETSC_IGNORE , PETSC_IGNORE , PETSC_IGNORE , PETSC_IGNORE , PETSC_IGNORE , PETSC_IGNORE ) ; CHKERRQ ( ierr ) ; double surface = ( double ) ( My * Mz ) * hy * hz ; for ( auto i = 0 ; i < numSpecies ; ++ i ) { totalConcData [ i ] /= surface ; } double fluence = fluxHandler -> getFluence ( ) ; std :: cout << " \n " << time << std :: endl ; for ( auto id = core :: network :: SpeciesId ( numSpecies ) ; id ; ++ id ) { std :: cout << network . getSpeciesName ( id ) << " " << totalConcData [ id ( ) ] << ' \n ' ; } std :: cout << " " << fluence << " \n " << std :: endl ; std :: ofstream outputFile ; outputFile . open ( " " , std :: ios :: app ) ; outputFile << fluence << " " ; for ( auto i = 0 ; i < numSpecies ; ++ i ) { outputFile << totalConcData [ i ] << " " ; } outputFile << std :: endl ; outputFile . close ( ) ; } ierr = DMDAVecRestoreArrayDOFRead ( da , solution , & solutionArray ) ; CHKERRQ ( ierr ) ; PetscFunctionReturn ( 0 ) ; }
This is a monitoring method that will compute the helium retention.
This is a monitoring method that will compute the helium retention .
// Initial declarations // Get the solver handler and local coordinates // Get the flux handler that will be used to compute fluxes. // Get the da from ts // Get the physical grid in the x direction // Get the network // Setup step size variables // Get the array of concentration // Store the concentration over the grid // Loop on the grid // Get the surface position // Boundary conditions // Get the pointer to the beginning of the solution data for // this grid point // Get the total concentrations at this grid point // Get the current process ID // Sum all the concentrations through MPI reduce // Master process // Get the total size of the grid rescale the concentrations // Compute the total surface irradiated by the helium flux // Rescale the concentration // Get the fluence // Print the result // Uncomment to write the retention and the fluence in a file // Restore the solutionArray