Commit 6473dfb0 authored by winklehner_d's avatar winklehner_d

Tried to merge master and local changes to ParallelCyclotronTracker and...

Tried to merge master and local changes to ParallelCyclotronTracker and MGPoissonSolver...hope it worked :) -DW
parent 40dc4210
......@@ -2689,6 +2689,8 @@ void ParallelCyclotronTracker::Tracker_Generic() {
// Total number of particles determines single particle mode (SPM, 1 particle) or
// Static Equilibrium Orbit Mode (SEO, 2 particles)
// *gmsg << "*** Running Tracker_Generic ***" << endl;
//if (timeIntegrator_m == 1) BorisPusher pusher; // Create a BorisPusher only for LF-2 method
BorisPusher pusher; // TEMP: OPAL will not compile without a pusher. -DW
......@@ -3589,35 +3591,37 @@ void ParallelCyclotronTracker::Tracker_Generic() {
// Here is global frame, don't do itsBunch->boundp();
// Check separately for phase space (ps) and statistics (stat) data dump frequency
if((((step_m + 1) % Options::psDumpFreq == 0) && initialTotalNum_m != 2)
|| (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) {
if (itsBunch->getTotalNum()>0) { // Only dump last step if we have particles left.
// Check separately for phase space (ps) and statistics (stat) data dump frequency
if((((step_m + 1) % Options::psDumpFreq == 0) && initialTotalNum_m != 2)
|| (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) {
// Write phase space data to h5 file
bunchDumpPhaseSpaceData();
}
// Write phase space data to h5 file
bunchDumpPhaseSpaceData();
}
if((((step_m + 1) % Options::statDumpFreq == 0) && initialTotalNum_m != 2)
|| (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) {
if((((step_m + 1) % Options::statDumpFreq == 0) && initialTotalNum_m != 2)
|| (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) {
// Write statistics data to stat file
bunchDumpStatData();
// Write statistics data to stat file
bunchDumpStatData();
}
}
if(!(step_m + 1 % 1000))
*gmsg << "Step " << step_m + 1 << endl;
} // end for: the integration is DONE after maxSteps_m steps!
} // end for: the integration is DONE after maxSteps_m steps or if all particles are lost!
// Some post-integration stuff
*gmsg << endl;
*gmsg << "* ---------------------------- DONE TRACKING PARTICLES -------------------------------- * " << endl;
if (!(itsBunch->getTotalNum()>0))
throw OpalException("ParallelCyclotronTracker::Tracker_MTS()", "No particles anymore, give up now!");
// TODO: This is totally wrecking saving paricles that were lost on collimators, etc.
// if no particles remain! -DW
//if (!(itsBunch->getTotalNum()>0))
// throw OpalException("ParallelCyclotronTracker::Tracker_Generic()", "No particles anymore, give up now!");
// Calculate tunes after tracking.
for(size_t ii = 0; ii < (itsBunch->getLocalNum()); ii++) {
if(itsBunch->ID[ii] == 0) {
double FinalMomentum2 = pow(itsBunch->P[ii](0), 2.0) + pow(itsBunch->P[ii](1), 2.0) + pow(itsBunch->P[ii](2), 2.0);
......@@ -3651,17 +3655,26 @@ void ParallelCyclotronTracker::Tracker_Generic() {
if(initialTotalNum_m == 1) closeFiles();
// Print out the Bunch information at end of the run. Because the bunch information
// displays in units of m we have to change back and forth one more time.
// Furthermore it is my opinion that the same units should be used throughout OPAL. -DW
itsBunch->R *= Vector_t(0.001); // mm --> m
*gmsg << endl << "* *********************** Bunch information in global frame: ***********************";
itsBunch->calcBeamParameters_cycl();
if (itsBunch->getTotalNum()>0){
// Print out the Bunch information at end of the run. Because the bunch information
// displays in units of m we have to change back and forth one more time.
// Furthermore it is my opinion that the same units should be used throughout OPAL. -DW
itsBunch->R *= Vector_t(0.001); // mm --> m
*gmsg << endl << "* *********************** Bunch information in global frame: ***********************";
*gmsg << *itsBunch << endl;
itsBunch->calcBeamParameters_cycl();
itsBunch->R *= Vector_t(1000.0); // m --> mm
*gmsg << *itsBunch << endl;
itsBunch->R *= Vector_t(1000.0); // m --> mm
} else {
*gmsg << endl << "* No Particles left in bunch!" << endl;
*gmsg << "* **********************************************************************************" << endl;
}
}
bool ParallelCyclotronTracker::getFieldsAtPoint(const double &t, const size_t &Pindex, Vector_t &Efield, Vector_t &Bfield) {
......
......@@ -265,14 +265,17 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
NDIndex<3> localId = layout_m->getLocalNDIndex();
IpplTimings::startTimer(FunctionTimer1_m);
// Compute boundary intersections (only do on the first step)
if(!itsBunch_m->getLocalTrackStep())
bp->Compute(hr, localId);
IpplTimings::stopTimer(FunctionTimer1_m);
// Define the Map
INFOMSG(level2 << "* Computing Map..." << endl);
IpplTimings::startTimer(FunctionTimer2_m);
computeMap(localId);
IpplTimings::stopTimer(FunctionTimer2_m);
INFOMSG(level2 << "* Done." << endl);
// Allocate the RHS with the new Epetra Map
if (Teuchos::is_null(RHS))
......@@ -280,30 +283,50 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
RHS->PutScalar(0.0);
// get charge densities from IPPL field and store in Epetra vector (RHS)
Ippl::Comm->barrier();
std::cout << "* Node:" << Ippl::myNode() << ", Filling RHS..." << std::endl;
Ippl::Comm->barrier();
IpplTimings::startTimer(FunctionTimer3_m);
int id = 0;
float scaleFactor = itsBunch_m->getdT();
std::cout << "* Node:" << Ippl::myNode() << ", Rho for final element: " << rho[localId[0].last()][localId[1].last()][localId[2].last()].get() << std::endl;
Ippl::Comm->barrier();
std::cout << "* Node:" << Ippl::myNode() << ", Local nx*ny*nz = " << localId[2].last() * localId[0].last() * localId[1].last() << std::endl;
std::cout << "* Node:" << Ippl::myNode() << ", Number of reserved local elements in RHS: " << RHS->MyLength() << std::endl;
std::cout << "* Node:" << Ippl::myNode() << ", Number of reserved global elements in RHS: " << RHS->GlobalLength() << std::endl;
Ippl::Comm->barrier();
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
if (bp->isInside(idx, idy, idz))
RHS->Values()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
if (bp->isInside(idx, idy, idz))
RHS->Values()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
}
}
}
std::cout << "* Node:" << Ippl::myNode() << ", Number of Local Inside Points " << id << std::endl;
Ippl::Comm->barrier();
IpplTimings::stopTimer(FunctionTimer3_m);
std::cout << "* Node:" << Ippl::myNode() << ", Done." << std::endl;
Ippl::Comm->barrier();
// build discretization matrix
INFOMSG(level2 << "* Building Discretization Matrix..." << endl);
IpplTimings::startTimer(FunctionTimer4_m);
if(Teuchos::is_null(A))
A = rcp(new Epetra_CrsMatrix(Copy, *Map, 7, true));
ComputeStencil(hr, RHS);
IpplTimings::stopTimer(FunctionTimer4_m);
INFOMSG(level2 << "* Done." << endl);
#ifdef DBG_STENCIL
EpetraExt::RowMatrixToMatlabFile("DiscrStencil.dat", *A);
#endif
INFOMSG(level2 << "* Computing Preconditioner..." << endl);
IpplTimings::startTimer(FunctionTimer5_m);
if(!MLPrec) {
MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
......@@ -312,11 +335,12 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
} else if (precmode_m == REUSE_PREC){
}
IpplTimings::stopTimer(FunctionTimer5_m);
INFOMSG(level2 << "* Done." << endl);
// setup preconditioned iterative solver
// use old LHS solution as initial guess
INFOMSG(level2 << "* Final Setup of Solver..." << endl);
IpplTimings::startTimer(FunctionTimer6_m);
problem_ptr->setOperator(A);
problem_ptr->setLHS(LHS);
problem_ptr->setRHS(RHS);
......@@ -330,15 +354,19 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
}
}
IpplTimings::stopTimer(FunctionTimer6_m);
INFOMSG(level2 << "* Done." << endl);
std::ofstream timings;
char filename[50];
sprintf(filename, "timing_MX%d_MY%d_MZ%d_nProc%d_recB%d_numB%d_nLHS%d", orig_nr_m[0], orig_nr_m[1], orig_nr_m[2], Comm.NumProc(), recycleBlocks_m, numBlocks_m, nLHS_m);
if(Comm.MyPID() == 0) timings.open(filename, std::ios::app);
double time = MPI_Wtime();
INFOMSG(level2 << "* Solving for Space Charge..." << endl);
IpplTimings::startTimer(FunctionTimer7_m);
solver_ptr->solve();
IpplTimings::stopTimer(FunctionTimer7_m);
INFOMSG(level2 << "* Done." << endl);
time = MPI_Wtime() - time;
double minTime, maxTime, avgTime;
......@@ -394,7 +422,7 @@ void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
if (bp->isInside(idx, idy, idz))
if (bp->isInside(idx, idy, idz))
numMyGridPoints++;
}
}
......@@ -412,7 +440,7 @@ void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
if (bp->isInside(idx, idy, idz)) {
if (bp->isInside(idx, idy, idz)) {
v[0] = (double)idx;
v[stride] = (double)idy;
v[stride2] = (double)idz;
......@@ -446,13 +474,14 @@ void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
}
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
int NumMyElements = 0;
std::vector<int> MyGlobalElements;
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
if (bp->isInside(idx, idy, idz)) {
if (bp->isInside(idx, idy, idz)) {
MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
NumMyElements++;
}
......
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