diff --git a/Src/GNUmakefile b/Src/GNUmakefile index 7df5fe0..5d4ac7d 100644 --- a/Src/GNUmakefile +++ b/Src/GNUmakefile @@ -10,7 +10,6 @@ USE_OMP = FALSE #EBASE = combinePlts #EBASE = isosurface #EBASE = stream -EBASE = stream2plt #EBASE = sampleStreamlines #EBASE = streamScatter #EBASE = streamSub diff --git a/Src/integrateNonEBUpdate.cpp b/Src/integrateNonEBUpdate.cpp new file mode 100644 index 0000000..6d40640 --- /dev/null +++ b/Src/integrateNonEBUpdate.cpp @@ -0,0 +1,153 @@ +/* + * ===================================================================================== + * + * Filename: integrateNonEBUpdate + * + * Description: Integral tool that takes variables from a plotfile and returns volume integral of that variable(s)-can run in parallel + * + * Version: 1.0 + * Created: 03/10/2020 01:37:30 PM + * Revision: none + * Compiler: gcc + * + * Author: Victor H. Zendejas Lopez + * Organization: CCSE @LBL + * + * ===================================================================================== + */ + + +#include +#include +#include +#include +#include + +using namespace amrex; +using namespace std; + +static void PrintUsage (char* progName) +{ + cout << "\nUsage:\n" + << progName + << "\n\tinfile = inputFileName comps=[n1 n2 ...]" + << "\n\t[-help]" + << "\n\n"; + exit(1); +} + + +static Real SumThisComp(AmrData &amrData, string whichVar , int finestLevel) +{ + + Real sum(0.0); + + // Get levels of refinement + int nLevels = finestLevel+1; + // Intersect variablenumber (must come after stoichiometry) + Vector> mf(nLevels); + for (int iLevel=0; iLevel < nLevels; iLevel++){ + Print()<< "Reading Level " << iLevel << "\n"; + int ngrow(0); + + // Make space for each component and one for the intersect flag + DistributionMapping dm(amrData.boxArray(iLevel)); + mf[iLevel].reset(new MultiFab(amrData.boxArray(iLevel), dm, 1, ngrow)); + amrData.FillVar(*mf[iLevel], iLevel, whichVar); + } + + for (int iLevel=0; iLevelboxArray(); + baf.coarsen(amrData.RefRatio()[iLevel]); + + for (MFIter mfi(*mf[iLevel]); mfi.isValid(); ++mfi) { + FArrayBox& myFab = (*mf[iLevel])[mfi]; + int idx = mfi.index(); + const Box& bx = mfi.validbox(); + + vector < pair < int,Box > > isects = baf.intersections(mf[iLevel]->boxArray()[idx]); + for (int ii = 0; ii < isects.size(); ii++){ + myFab.setVal(0,isects[ii].second,0,1); + } + + sum +=myFab.sum(0,1); + } + } + + ParallelDescriptor::ReduceRealSum(sum); + amrData.FlushGrids(amrData.StateNumber(whichVar)); + return sum; +} + + +int main(int argc, char* argv[]) +{ + Initialize(argc,argv); + { + if(argc == 1){ + // PrintUsage(argv[0]); + } + + ParallelDescriptor::StartParallel(&argc, &argv); + ParmParse pp; + + string Output_Variables; + //Count the number of terms in the Output_Varibales + int nOutput_Variables = (pp.countval("Output_Variables")); + Vector < string > Out_Var(nOutput_Variables); // Creates a vector Out_var length + // Loop over Output_Variables to extract information from each element in Out_Var + + for(int ivarOut = 0; ivarOut < nOutput_Variables; ++ivarOut){ + pp.get("Output_Variables", Out_Var[ivarOut],ivarOut); + } + string Output_Format; + pp.get("Output_Format", Output_Format); + string Plot_File; + pp.get("Plot_File", Plot_File); + string Output_File; + pp.get("Output_File", Output_File); + + // Extract Plot data + DataServices::SetBatchMode(); + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices plotDataType(Plot_File, fileType); + + if (!plotDataType.AmrDataOk()) // If nothing is stored within plotDataType then run this loop + DataServices::Dispatch(DataServices::ExitRequest, NULL); + AmrData& amrData = plotDataType.AmrDataRef(); + cout << " ...done." << endl; + + // Most likely, the variable was averaged down conservatively, lev 0 may be all that is reqd + int finestLevel = amrData.FinestLevel(); + + Vector whichVar(nOutput_Variables); + // Read in variable list - adjusting for stoichiometry + for(int v = 0; v < nOutput_Variables; v++) { + pp.get("Output_Variables", whichVar[v], v); + } + + // Check the names of the variables are present in the plotfile + Vector < int > destFills(nOutput_Variables); + for (int v = 0; v < nOutput_Variables; v++){ + destFills[v] = v; + if (amrData.StateNumber(whichVar[v])<0) { + string message="Bad variable name ("+whichVar[v]+")"; + Error(message.c_str()); + } + } + + if (ParallelDescriptor::IOProcessor()){ + cout << "TIME = " << setw(15) << amrData.Time(); + } + + for (int v = 0; v Line; + +Line::iterator FindMySeg(Line& segs, int idx) +{ + for (Line::iterator it=segs.begin(); it!=segs.end(); ++it) + { + if ( ((*it).ID_l() == idx) || ((*it).ID_r() == idx) ) + return it; + } + return segs.end(); +} + +bool operator==(const Seg& lhs, const Seg& rhs) +{ + return lhs.ID_l() == rhs.ID_l() && lhs.ID_r() == rhs.ID_r(); +} + +bool operator!=(const Seg& lhs, const Seg& rhs) +{ + return !operator==(lhs,rhs); +} + +typedef Vector Point; + +std::pair +nodesAreClose(const FArrayBox& nodes, + int lhs, + int rhs, + Real eps) +{ + Real sep = 0; + int nComp = nodes.nComp(); + Point ret(nComp,0); + for (int i=0; i(std::sqrt(sep) <= eps, ret); +} + +void +join(Line& l1, + Line& l2, + FArrayBox& nodes, + Point& pt, + list& newPts, + list& deadPts) +{ + int newPtID = nodes.box().numPts() + newPts.size(); + newPts.push_back(pt); + + int rhs = l1.back().ID_l(); + deadPts.push_back(l1.back().ID_r()); + BL_ASSERT(l1.size()>0); + l1.pop_back(); + l1.push_back(Seg(rhs,newPtID)); + + int lhs = l2.front().ID_r(); + deadPts.push_back(l2.front().ID_l()); + BL_ASSERT(l2.size()>0); + l2.pop_front(); + l2.push_front(Seg(newPtID,lhs)); + + l1.splice(l1.end(),l2); +} + +// Provides an ordering to uniquely sort segments +struct SegLT +{ + bool operator()(const Seg& lhs, + const Seg& rhs) + { + int smL = std::min(lhs.ID_l(),lhs.ID_r()); + int smR = std::min(rhs.ID_l(),rhs.ID_r()); + if (smL < smR) + { + return true; + } + else if (smL == smR) + { + int lgL = std::max(lhs.ID_l(),lhs.ID_r()); + int lgR = std::max(rhs.ID_l(),rhs.ID_r()); + return lgL < lgR; + } + return false; + } +}; + +using std::list; + +list +MakeCLines(const FArrayBox& nodes, + Vector& faceData, + int nodesPerElt) +{ + // We want to clean things up a little...the contour is now in a form + // of a huge list of itty-bitty segments. We want to connect up the segments + // to form a minimal number of polylines. + + // First build up a Line + int nElts = faceData.size() / nodesPerElt; + BL_ASSERT(nElts * nodesPerElt == faceData.size()); + Line segList; + for (int i=0; i cLines; + cLines.push_back(Line()); + + while (segList.begin() != segList.end()) + { + Line::iterator segIt = FindMySeg(segList,idx); + if (segIt != segList.end()) + { + int idx_l = (*segIt).ID_l(); + int idx_r = (*segIt).ID_r(); + + if ( idx_l == idx ) + { + idx = idx_r; + cLines.back().push_back(*segIt); + } + else + { + idx = idx_l; + Seg nseg(*segIt); nseg.flip(); + cLines.back().push_back(nseg); + } + + segList.erase(segIt); + } + else + { + cLines.push_back(Line()); + idx = segList.front().ID_r(); + segList.pop_front(); + } + } + + // Connect up the line fragments as much as possible + bool changed; + do + { + changed = false; + for (std::list::iterator it = cLines.begin(); it!=cLines.end(); ++it) + { + if (!(*it).empty()) + { + const int idx_l = (*it).front().ID_l(); + const int idx_r = (*it).back().ID_r(); + for (std::list::iterator it1 = cLines.begin(); it1!=cLines.end(); ++it1) + { + if (!(*it1).empty() && (*it).front()!=(*it1).front()) + { + if (idx_r == (*it1).front().ID_l()) + { + (*it).splice((*it).end(),*it1); + changed = true; + } + else if (idx_r == (*it1).back().ID_r()) + { + (*it1).reverse(); + for (Line::iterator it2=(*it1).begin(); it2!=(*it1).end(); ++it2) + (*it2).flip(); + (*it).splice((*it).end(),*it1); + changed = true; + } + else if (idx_l == (*it1).front().ID_l()) + { + (*it1).reverse(); + for (Line::iterator it2=(*it1).begin(); it2!=(*it1).end(); ++it2) + (*it2).flip(); + (*it).splice((*it).begin(),*it1); + changed = true; + } + } + } + } + } + } while(changed); + + // Clear out empty placeholders for lines we connected up to others. + for (std::list::iterator it = cLines.begin(); it!=cLines.end();) + { + if (it->empty()) + cLines.erase(it++); + else + it++; + } + + // Dump out the final number of contours + cerr << " number of contour lines: " << cLines.size() << endl; + + return cLines; +} + int main (int argc, char* argv[]) @@ -1163,7 +1388,7 @@ main (int argc, { const auto& f = state.array(mfi); const Box& gbox = state[mfi].box(); - const int ratio = pf.refRatio(lev-1); + const int ratio = (lev == 0 ? 1 : pf.refRatio(lev-1)); Real loc[BL_SPACEDIM]; if (lev!=0) { @@ -1171,18 +1396,18 @@ main (int argc, AMREX_PARALLEL_FOR_3D ( gbox, i, j, k, { - f(i,j,k,0) = ((i < 0 ? i*ri-1 : i*ri) + 0.5)*dxf[0]*ratio + plo[0]; - f(i,j,k,1) = ((j < 0 ? j*ri-1 : j*ri) + 0.5)*dxf[1]*ratio + plo[1]; - f(i,j,k,2) = ((k < 0 ? k*ri-1 : k*ri) + 0.5)*dxf[2]*ratio + plo[2]; + AMREX_D_TERM(f(i,j,k,0) = (static_cast(i < 0 ? i*ri-1 : i*ri) + 0.5)*dxf[0]*ratio + plo[0];, + f(i,j,k,1) = (static_cast(j < 0 ? j*ri-1 : j*ri) + 0.5)*dxf[1]*ratio + plo[1];, + f(i,j,k,2) = (static_cast(k < 0 ? k*ri-1 : k*ri) + 0.5)*dxf[2]*ratio + plo[2];); }); } const Box& vbox = gridArray[lev][mfi.index()]; AMREX_PARALLEL_FOR_3D ( vbox, i, j, k, { - f(i,j,k,0) = (i + 0.5)*dxf[0] + plo[0]; - f(i,j,k,1) = (j + 0.5)*dxf[1] + plo[1]; - f(i,j,k,2) = (k + 0.5)*dxf[2] + plo[2]; + AMREX_D_TERM(f(i,j,k,0) = (i + 0.5)*dxf[0] + plo[0];, + f(i,j,k,1) = (j + 0.5)*dxf[1] + plo[1];, + f(i,j,k,2) = (k + 0.5)*dxf[2] + plo[2];); }); } @@ -1716,6 +1941,57 @@ main (int argc, const FABdata& tmpData = *tmpDataP; + /* + In 2D we can connect up the line segments in order to get a set of + lines. We do this and hack in an output function that creates a + file that ParaView can read. May need to split out one zone per + file...not sure. Also, ParaView doesn't seem to like variables + with names like Y(XX), so we translate them to XX. + */ +#if AMREX_SPACEDIM==2 + list contours = MakeCLines(tmpData.getFab(),eltRaw,nodesPerElt); + + std::string fileName = "MARC.dat"; + std::ofstream ofm(fileName.c_str()); + + // Build variable name string + std::string mvars; + for (int j=0; jID_l(); + const Real* vec = sortedNodes[nodeIdx]->m_vec; + for (int n=0; nm_vec; + for (int n=0; n +#include +#include + +#include +#include +#include +#include +#include + +using namespace amrex; +using namespace std; + + +int main (int argc, char* argv[]) +{ + + Initialize(argc,argv); + { + + // ParmParse is a class that has stored functions + ParmParse pp; + // declares Output_Variables as an string-comes from the input file + string Output_Variables; + //Count the number of terms in the Output_Varibales + int nOutput_Variables = (pp.countval("Output_Variables")); + Vector < string > Out_Var(nOutput_Variables); // Creates a vector Out_var length nOU..es + // Loop over Output_Variables to extract information from each element in Out_Var + + for(int ivarOut = 0; ivarOut < nOutput_Variables; ++ivarOut){ + pp.get("Output_Variables", Out_Var[ivarOut],ivarOut); + } + string Output_Format; + pp.get("Output_Format", Output_Format); + string Plot_File; + pp.get("Plot_File", Plot_File); + + string Output_File; + pp.get("Output_File", Output_File); + // Extract Plot data + // 1. DataService-Class-setBatchMode is function in that class + // 2. FileType is a class, plotType is a variable that takes in arguments Amrvis::NEWPLT + //3. DataServices creates a variable that takesin Plot_File and plotType as an arguement + //4. amrData is a variable that contains the directory of where the data from + //the plot files are located. AmrDataRef is a fucntion acting on plotDataType. + + DataServices::SetBatchMode(); + Amrvis::FileType plotType(Amrvis::NEWPLT); + DataServices plotDataType(Plot_File, plotType); + + if (!plotDataType.AmrDataOk()) // If nothing is stored within plotDataType then run this loop + DataServices::Dispatch(DataServices::ExitRequest, NULL); + + AmrData& amrData = plotDataType.AmrDataRef(); + cout << " ...done." << endl; + + Vector whichVar(nOutput_Variables); + // Read in variable list - adjusting for stoichiometry + + for(int v = 0; v < nOutput_Variables; v++) { + pp.get("Output_Variables", whichVar[v], v); + } + // Check the names of the variables are present in the plotfile + Vector < int > destFills(nOutput_Variables); + + for (int v = 0; v < nOutput_Variables; v++){ + destFills[v] = v; + + if (amrData.StateNumber(whichVar[v])<0) { + string message="Bad variable name ("+whichVar[v]+")"; + Error(message.c_str()); + } + } + // Get domain size + Vector probLo=amrData.ProbLo(); + Vector probHi=amrData.ProbHi(); + // Get probDomain + Vector probDomain = amrData.ProbDomain(); + // Get levels of refinement + int finestLevel = amrData.FinestLevel(); + int nLevels = finestLevel+1; + // Intersect variable number (must come after stoichiometry) + int isVar=nOutput_Variables; + + // Load the data + // + // 1. boxArray is a class for storing collection of boxes on a single AMR level + // 2. DistributionMapping is a class that describes which process owns the data + // on the domain specified by the boxes in boxArray + // + // Vector> mf(nLevels); + Vector> mf(nLevels); + for (int iLevel=0; iLevel < nLevels; iLevel++){ + Print()<< "Reading Level " << iLevel << "\n"; + int ngrow(0); + + // Make space for each component and one for the intersect flag + DistributionMapping dm(amrData.boxArray(iLevel)); + mf[iLevel].reset(new MultiFab(amrData.boxArray(iLevel), dm, nOutput_Variables+1, ngrow)); + // Put the value 1 in the intersect flag + mf[iLevel]->setVal(1.,isVar,1); + // Load the data (make a copy of whichVar to make sure it's the same size as destFills) + + Vector loadWhichVar(nOutput_Variables); + for (int v = 0; v < nOutput_Variables; v++){ + loadWhichVar[v] = whichVar[v]; + } + amrData.FillVar(*mf[iLevel], iLevel, loadWhichVar, destFills); + } + + // Set the intersect flag to zero if necessary + for (int iLevel=0; iLevelboxArray(); + baf.coarsen(amrData.RefRatio()[iLevel]); + Print()<< "Tagging Level " << iLevel << "\n"; + + for (MFIter mfi(*mf[iLevel]); mfi.isValid(); ++mfi) { + FArrayBox& myFab = (*mf[iLevel])[mfi]; + int idx = mfi.index(); + vector < pair > isects = baf.intersections(mf[iLevel]->boxArray()[idx]); + for (int ii = 0; ii < isects.size(); ii++){ + myFab.setVal(0,isects[ii].second,isVar,1); + } + } + } +// for (int iLevel=0; iLevel < nLevels; iLevel++){ +// VisMF::Write(*mf[iLevel],"MFvis_Level"+to_string(iLevel)); +// } + + // Create an outputfile to write data to + ofstream ofs; + ofs.open(Output_File); + + // Load Header file for output file + Vector loadWhichVar(nOutput_Variables); + for (int v = 0; v < nOutput_Variables; v++){ + loadWhichVar[v] = whichVar[v]; + ofs << " " << loadWhichVar [v]; + } + ofs << "\n"; + + for (int iLevel=0; iLevel < nLevels; iLevel++){ + + for (MFIter mfi(*mf[iLevel]); mfi.isValid(); ++mfi) { + const Box& box = mfi.validbox(); + const auto& Tagfab = mf[iLevel]->array(mfi,isVar); + const auto& DataFab = mf[iLevel]->array(mfi); + + // Write data to output file + + AMREX_HOST_DEVICE_PARALLEL_FOR_3D(box,i,j,k, + { + + if( Tagfab(i,j,k) > 0.5){ + for(int n = 0; n < nOutput_Variables; ++n){ + ofs << " " << DataFab(i,j,k,n); + } + ofs << "\n"; + } + }); + } + } + ofs.close(); + } + Finalize(); + return 0; +} +