diff --git a/CMakeLists.txt b/CMakeLists.txt index 91e9c2ef..efe167d6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -96,7 +96,7 @@ include("cmake/cpm.cmake") CPMAddPackage( NAME ViennaCore - VERSION 1.2.0 + VERSION 1.2.1 GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore" OPTIONS "VIENNACORE_FORMAT_EXCLUDE app/" EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}) @@ -108,13 +108,13 @@ CPMAddPackage( CPMFindPackage( NAME ViennaRay - VERSION 3.1.4 + VERSION 3.2.0 GIT_REPOSITORY "https://github.com/ViennaTools/ViennaRay" EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}) CPMFindPackage( NAME ViennaLS - VERSION 4.2.1 + VERSION 4.2.2 GIT_REPOSITORY "https://github.com/ViennaTools/ViennaLS" EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}) diff --git a/include/viennaps/psProcess.hpp b/include/viennaps/psProcess.hpp index 619089db..d4367a0a 100644 --- a/include/viennaps/psProcess.hpp +++ b/include/viennaps/psProcess.hpp @@ -1,19 +1,12 @@ #pragma once +#include "psProcessBase.hpp" #include "psProcessModel.hpp" -#include "psProcessParams.hpp" #include "psTranslationField.hpp" -#include "psUnits.hpp" #include "psUtils.hpp" #include -#include -#include -#include -#include -#include -#include #include namespace viennaps { @@ -23,758 +16,42 @@ using namespace viennacore; /// This class server as the main process tool, applying a user- or pre-defined /// process model to a domain. Depending on the user inputs surface advection, a /// single callback function or a geometric advection is applied. -template class Process { - using translatorType = std::unordered_map; - using psDomainType = SmartPointer>; +template +class Process final : public ProcessBase { + // Typedefs + using TranslatorType = std::unordered_map; + using DomainType = SmartPointer>; public: - Process() {} - Process(psDomainType passedDomain) : domain(passedDomain) {} + Process() = default; + Process(DomainType domain) : ProcessBase(domain) {} // Constructor for a process with a pre-configured process model. - Process(psDomainType passedDomain, - SmartPointer> passedProcessModel, - const NumericType passedDuration = 0.) - : domain(passedDomain), model(passedProcessModel), - processDuration(passedDuration) {} + Process(DomainType domain, + SmartPointer> processModel, + const NumericType duration = 0.) + : ProcessBase(domain, processModel, duration), + processModel_(processModel) {} // Set the process model. This can be either a pre-configured process model or // a custom process model. A custom process model must interface the // psProcessModel class. - void setProcessModel( - SmartPointer> passedProcessModel) { - model = passedProcessModel; + void + setProcessModel(SmartPointer> processModel) { + processModel_ = processModel; + this->model_ = processModel; } - // Set the process domain. - void setDomain(psDomainType passedDomain) { domain = passedDomain; } - - /* ----- Process parameters ----- */ - - // Set the duration of the process. - void setProcessDuration(NumericType passedDuration) { - processDuration = passedDuration; - } - - // Returns the duration of the recently run process. This duration can - // sometimes slightly vary from the set process duration, due to the maximum - // time step according to the CFL condition. - NumericType getProcessDuration() const { return processTime; } - - /* ----- Ray tracing parameters ----- */ - - // Set the number of iterations to initialize the coverages. - void setMaxCoverageInitIterations(unsigned maxIt) { maxIterations = maxIt; } - - // Set the threshold for the coverage delta metric to reach convergence. - void setCoverageDeltaThreshold(NumericType treshold) { - coverageDeltaThreshold = treshold; - } - - // Set the source direction, where the rays should be traced from. - void setSourceDirection(viennaray::TraceDirection passedDirection) { - rayTracingParams.sourceDirection = passedDirection; - } - - // Specify the number of rays to be traced for each particle throughout the - // process. The total count of rays is the product of this number and the - // number of points in the process geometry. - void setNumberOfRaysPerPoint(unsigned numRays) { - rayTracingParams.raysPerPoint = numRays; - } - - // Enable flux smoothing. The flux at each surface point, calculated - // by the ray tracer, is averaged over the surface point neighbors. - void enableFluxSmoothing() { rayTracingParams.smoothingNeighbors = 1; } - - // Disable flux smoothing. - void disableFluxSmoothing() { rayTracingParams.smoothingNeighbors = 0; } - void setRayTracingDiskRadius(NumericType radius) { - rayTracingParams.diskRadius = radius; - if (rayTracingParams.diskRadius < 0.) { + rayTracingParams_.diskRadius = radius; + if (rayTracingParams_.diskRadius < 0.) { Logger::getInstance() .addWarning("Disk radius must be positive. Using default value.") .print(); - rayTracingParams.diskRadius = 0.; + rayTracingParams_.diskRadius = 0.; } } - void enableFluxBoundaries() { rayTracingParams.ignoreFluxBoundaries = false; } - - // Ignore boundary conditions during the flux calculation. - void disableFluxBoundaries() { rayTracingParams.ignoreFluxBoundaries = true; } - - // Enable the use of random seeds for ray tracing. This is useful to - // prevent the formation of artifacts in the flux calculation. - void enableRandomSeeds() { rayTracingParams.useRandomSeeds = true; } - - // Disable the use of random seeds for ray tracing. - void disableRandomSeeds() { rayTracingParams.useRandomSeeds = false; } - - void setRayTracingParameters( - const RayTracingParameters &passedRayTracingParams) { - rayTracingParams = passedRayTracingParams; - } - - auto &getRayTracingParameters() { return rayTracingParams; } - - /* ----- Advection parameters ----- */ - - // Set the integration scheme for solving the level-set equation. - // Possible integration schemes are specified in - // viennals::IntegrationSchemeEnum. - void setIntegrationScheme(IntegrationScheme passedIntegrationScheme) { - advectionParams.integrationScheme = passedIntegrationScheme; - } - - // Enable the output of the advection velocities on the level-set mesh. - void enableAdvectionVelocityOutput() { - advectionParams.velocityOutput = true; - } - - // Disable the output of the advection velocities on the level-set mesh. - void disableAdvectionVelocityOutput() { - advectionParams.velocityOutput = false; - } - - // Set the CFL (Courant-Friedrichs-Levy) condition to use during surface - // advection in the level-set. The CFL condition defines the maximum distance - // a surface is allowed to move in a single advection step. It MUST be below - // 0.5 to guarantee numerical stability. Defaults to 0.4999. - void setTimeStepRatio(NumericType cfl) { - advectionParams.timeStepRatio = cfl; - } - - void setAdvectionParameters( - const AdvectionParameters &passedAdvectionParams) { - advectionParams = passedAdvectionParams; - } - - auto &getAdvectionParameters() { return advectionParams; } - - /* ----- Process execution ----- */ - - // A single flux calculation is performed on the domain surface. The result is - // stored as point data on the nodes of the mesh. - SmartPointer> calculateFlux() { - model->initialize(domain, 0.); - const auto name = model->getProcessName().value_or("default"); - const auto logLevel = Logger::getLogLevel(); - - // Generate disk mesh from domain - auto mesh = SmartPointer>::New(); - viennals::ToDiskMesh meshConverter(mesh); - for (auto dom : domain->getLevelSets()) { - meshConverter.insertNextLevelSet(dom); - } - meshConverter.apply(); - - viennaray::Trace rayTracer; - initializeRayTracer(rayTracer); - - auto points = mesh->getNodes(); - auto normals = *mesh->getCellData().getVectorData("Normals"); - auto materialIds = *mesh->getCellData().getScalarData("MaterialIds"); - - if (rayTracingParams.diskRadius == 0.) { - rayTracer.setGeometry(points, normals, domain->getGrid().getGridDelta()); - } else { - rayTracer.setGeometry(points, normals, domain->getGrid().getGridDelta(), - rayTracingParams.diskRadius); - } - rayTracer.setMaterialIds(materialIds); - - const bool useProcessParams = - model->getSurfaceModel()->getProcessParameters() != nullptr; - bool useCoverages = false; - - // Initialize coverages - model->getSurfaceModel()->initializeSurfaceData(points.size()); - if (!coveragesInitialized_) - model->getSurfaceModel()->initializeCoverages(points.size()); - auto coverages = model->getSurfaceModel()->getCoverages(); - std::ofstream covMetricFile; - if (coverages != nullptr) { - Timer timer; - useCoverages = true; - Logger::getInstance().addInfo("Using coverages.").print(); - - if (maxIterations == std::numeric_limits::max() && - coverageDeltaThreshold == 0.) { - maxIterations = 10; - Logger::getInstance() - .addWarning("No coverage initialization parameters set. Using " + - std::to_string(maxIterations) + - " initialization iterations.") - .print(); - } - - if (!coveragesInitialized_) { - timer.start(); - Logger::getInstance().addInfo("Initializing coverages ... ").print(); - // debug output - if (logLevel >= 5) - covMetricFile.open(name + "_covMetric.txt"); - - for (unsigned iteration = 0; iteration < maxIterations; ++iteration) { - // We need additional signal handling when running the C++ code from - // the Python bindings to allow interrupts in the Python scripts -#ifdef VIENNAPS_PYTHON_BUILD - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - // save current coverages to compare with the new ones - auto prevStepCoverages = - SmartPointer>::New(*coverages); - - auto fluxes = - calculateFluxes(rayTracer, useCoverages, useProcessParams); - - // update coverages in the model - model->getSurfaceModel()->updateCoverages(fluxes, materialIds); - - // metric to check for convergence - coverages = model->getSurfaceModel()->getCoverages(); - auto metric = - calculateCoverageDeltaMetric(coverages, prevStepCoverages); - - if (logLevel >= 3) { - mergeScalarData(mesh->getCellData(), coverages); - mergeScalarData(mesh->getCellData(), fluxes); - auto surfaceData = model->getSurfaceModel()->getSurfaceData(); - if (surfaceData) - mergeScalarData(mesh->getCellData(), surfaceData); - printDiskMesh(mesh, name + "_covInit_" + std::to_string(iteration) + - ".vtp"); - - Logger::getInstance() - .addInfo("Iteration: " + std::to_string(iteration + 1)) - .print(); - - std::stringstream stream; - stream << std::setprecision(4) << std::fixed; - stream << "Coverage delta metric: "; - for (int i = 0; i < coverages->getScalarDataSize(); i++) { - stream << coverages->getScalarDataLabel(i) << ": " << metric[i] - << "\t"; - } - Logger::getInstance().addInfo(stream.str()).print(); - - // log metric to file - if (logLevel >= 5) { - for (auto val : metric) { - covMetricFile << val << ";"; - } - covMetricFile << "\n"; - } - } - - if (checkCoveragesConvergence(metric)) { - Logger::getInstance() - .addInfo("Coverages converged after " + - std::to_string(iteration + 1) + " iterations.") - .print(); - break; - } - } // end coverage initialization loop - coveragesInitialized_ = true; - - if (logLevel >= 5) - covMetricFile.close(); - - timer.finish(); - Logger::getInstance() - .addTiming("Coverage initialization", timer) - .print(); - } - } // end coverage initialization - - auto fluxes = calculateFluxes(rayTracer, useCoverages, useProcessParams); - mergeScalarData(mesh->getCellData(), fluxes); - auto surfaceData = model->getSurfaceModel()->getSurfaceData(); - if (surfaceData) - mergeScalarData(mesh->getCellData(), surfaceData); - - return mesh; - } - - // Run the process. - void apply() { - /* ---------- Check input --------- */ - if (!domain) { - Logger::getInstance().addWarning("No domain passed to Process.").print(); - return; - } - - if (domain->getLevelSets().empty()) { - Logger::getInstance().addWarning("No level sets in domain.").print(); - return; - } - - if (!model) { - Logger::getInstance() - .addWarning("No process model passed to Process.") - .print(); - return; - } - model->initialize(domain, processDuration); - const auto name = model->getProcessName().value_or("default"); - - if (model->getGeometricModel()) { - model->getGeometricModel()->setDomain(domain); - Logger::getInstance().addInfo("Applying geometric model...").print(); - model->getGeometricModel()->apply(); - return; - } - - if (processDuration == 0.) { - // apply only advection callback - if (model->getAdvectionCallback()) { - model->getAdvectionCallback()->setDomain(domain); - model->getAdvectionCallback()->applyPreAdvect(0); - } else { - Logger::getInstance() - .addWarning("No advection callback passed to Process.") - .print(); - } - return; - } - - if (!model->getSurfaceModel()) { - Logger::getInstance() - .addWarning("No surface model passed to Process.") - .print(); - return; - } - - if (!model->getVelocityField()) { - Logger::getInstance() - .addWarning("No velocity field passed to Process.") - .print(); - return; - } - - /* ------ Process Setup ------ */ - const unsigned int logLevel = Logger::getLogLevel(); - Timer processTimer; - processTimer.start(); - - double remainingTime = processDuration; - const NumericType gridDelta = domain->getGrid().getGridDelta(); - - auto diskMesh = SmartPointer>::New(); - auto translator = SmartPointer::New(); - viennals::ToDiskMesh meshConverter(diskMesh); - meshConverter.setTranslator(translator); - if (domain->getMaterialMap() && - domain->getMaterialMap()->size() == domain->getLevelSets().size()) { - meshConverter.setMaterialMap(domain->getMaterialMap()->getMaterialMap()); - } - - auto transField = SmartPointer>::New( - model->getVelocityField(), domain->getMaterialMap()); - transField->setTranslator(translator); - - viennals::Advect advectionKernel; - advectionKernel.setVelocityField(transField); - advectionKernel.setIntegrationScheme(advectionParams.integrationScheme); - advectionKernel.setTimeStepRatio(advectionParams.timeStepRatio); - advectionKernel.setSaveAdvectionVelocities(advectionParams.velocityOutput); - advectionKernel.setDissipationAlpha(advectionParams.dissipationAlpha); - advectionKernel.setIgnoreVoids(advectionParams.ignoreVoids); - advectionKernel.setCheckDissipation(advectionParams.checkDissipation); - // normals vectors are only necessary for analytical velocity fields - if (model->getVelocityField()->getTranslationFieldOptions() > 0) - advectionKernel.setCalculateNormalVectors(false); - - for (auto dom : domain->getLevelSets()) { - meshConverter.insertNextLevelSet(dom); - advectionKernel.insertNextLevelSet(dom); - } - - /* --------- Setup for ray tracing ----------- */ - const bool useRayTracing = !model->getParticleTypes().empty(); - - viennaray::BoundaryCondition rayBoundaryCondition[D]; - viennaray::Trace rayTracer; - - if (useRayTracing) { - initializeRayTracer(rayTracer); - - // initialize particle data logs - particleDataLogs.resize(model->getParticleTypes().size()); - for (std::size_t i = 0; i < model->getParticleTypes().size(); i++) { - int logSize = model->getParticleLogSize(i); - if (logSize > 0) { - particleDataLogs[i].data.resize(1); - particleDataLogs[i].data[0].resize(logSize); - } - } - } - - // Determine whether advection callback is used - const bool useAdvectionCallback = model->getAdvectionCallback() != nullptr; - if (useAdvectionCallback) { - model->getAdvectionCallback()->setDomain(domain); - } - - // Determine whether there are process parameters used in ray tracing - model->getSurfaceModel()->initializeProcessParameters(); - const bool useProcessParams = - model->getSurfaceModel()->getProcessParameters() != nullptr; - - if (useProcessParams) - Logger::getInstance().addInfo("Using process parameters.").print(); - if (useAdvectionCallback) - Logger::getInstance().addInfo("Using advection callback.").print(); - - bool useCoverages = false; - - // Initialize coverages - meshConverter.apply(); - auto numPoints = diskMesh->getNodes().size(); - model->getSurfaceModel()->initializeSurfaceData(numPoints); - if (!coveragesInitialized_) - model->getSurfaceModel()->initializeCoverages(numPoints); - auto coverages = model->getSurfaceModel()->getCoverages(); - std::ofstream covMetricFile; - if (coverages != nullptr) { - Timer timer; - useCoverages = true; - Logger::getInstance().addInfo("Using coverages.").print(); - - if (maxIterations == std::numeric_limits::max() && - coverageDeltaThreshold == 0.) { - maxIterations = 10; - Logger::getInstance() - .addWarning("No coverage initialization parameters set. Using " + - std::to_string(maxIterations) + - " initialization iterations.") - .print(); - } - - // debug output - if (logLevel >= 5) - covMetricFile.open(name + "_covMetric.txt"); - - if (!coveragesInitialized_) { - timer.start(); - Logger::getInstance().addInfo("Initializing coverages ... ").print(); - - auto points = diskMesh->getNodes(); - auto normals = *diskMesh->getCellData().getVectorData("Normals"); - auto materialIds = - *diskMesh->getCellData().getScalarData("MaterialIds"); - if (rayTracingParams.diskRadius == 0.) { - rayTracer.setGeometry(points, normals, gridDelta); - } else { - rayTracer.setGeometry(points, normals, gridDelta, - rayTracingParams.diskRadius); - } - rayTracer.setMaterialIds(materialIds); - - model->initialize(domain, -1.); - - for (unsigned iteration = 0; iteration < maxIterations; ++iteration) { - // We need additional signal handling when running the C++ code from - // the Python bindings to allow interrupts in the Python scripts -#ifdef VIENNAPS_PYTHON_BUILD - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - // save current coverages to compare with the new ones - auto prevStepCoverages = - SmartPointer>::New(*coverages); - - auto fluxes = - calculateFluxes(rayTracer, useCoverages, useProcessParams); - - // update coverages in the model - model->getSurfaceModel()->updateCoverages(fluxes, materialIds); - - // metric to check for convergence - coverages = model->getSurfaceModel()->getCoverages(); - auto metric = - calculateCoverageDeltaMetric(coverages, prevStepCoverages); - - if (logLevel >= 3) { - mergeScalarData(diskMesh->getCellData(), coverages); - mergeScalarData(diskMesh->getCellData(), fluxes); - auto surfaceData = model->getSurfaceModel()->getSurfaceData(); - if (surfaceData) - mergeScalarData(diskMesh->getCellData(), surfaceData); - printDiskMesh(diskMesh, name + "_covInit_" + - std::to_string(iteration) + ".vtp"); - Logger::getInstance() - .addInfo("Iteration: " + std::to_string(iteration + 1)) - .print(); - - std::stringstream stream; - stream << std::setprecision(4) << std::fixed; - stream << "Coverage delta metric: "; - for (int i = 0; i < coverages->getScalarDataSize(); i++) { - stream << coverages->getScalarDataLabel(i) << ": " << metric[i] - << "\t"; - } - Logger::getInstance().addInfo(stream.str()).print(); - - // log metric to file - if (logLevel >= 5) { - for (auto val : metric) { - covMetricFile << val << ";"; - } - covMetricFile << "\n"; - } - } - - if (checkCoveragesConvergence(metric)) { - Logger::getInstance() - .addInfo("Coverages converged after " + - std::to_string(iteration + 1) + " iterations.") - .print(); - break; - } - } // end coverage initialization loop - coveragesInitialized_ = true; - - timer.finish(); - Logger::getInstance() - .addTiming("Coverage initialization", timer) - .print(); - } - } // end coverage initialization - - double previousTimeStep = 0.; - size_t counter = 0; - unsigned lsVelCounter = 0; - Timer rtTimer; - Timer callbackTimer; - Timer advTimer; - while (remainingTime > 0.) { - // We need additional signal handling when running the C++ code from the - // Python bindings to allow interrupts in the Python scripts -#ifdef VIENNAPS_PYTHON_BUILD - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - // Expand LS based on the integration scheme - advectionKernel.prepareLS(); - model->initialize(domain, remainingTime); - - auto fluxes = SmartPointer>::New(); - meshConverter.apply(); - auto materialIds = *diskMesh->getCellData().getScalarData("MaterialIds"); - auto points = diskMesh->getNodes(); - - // rate calculation by top-down ray tracing - if (useRayTracing) { - rtTimer.start(); - auto normals = *diskMesh->getCellData().getVectorData("Normals"); - if (rayTracingParams.diskRadius == 0.) { - rayTracer.setGeometry(points, normals, gridDelta); - } else { - rayTracer.setGeometry(points, normals, gridDelta, - rayTracingParams.diskRadius); - } - rayTracer.setMaterialIds(materialIds); - - fluxes = calculateFluxes(rayTracer, useCoverages, useProcessParams); - - rtTimer.finish(); - Logger::getInstance() - .addTiming("Top-down flux calculation", rtTimer) - .print(); - } - - // update coverages and calculate coverage delta metric - if (useCoverages) { - coverages = model->getSurfaceModel()->getCoverages(); - auto prevStepCoverages = - SmartPointer>::New(*coverages); - - // update coverages in the model - model->getSurfaceModel()->updateCoverages(fluxes, materialIds); - - if (coverageDeltaThreshold > 0) { - auto metric = - calculateCoverageDeltaMetric(coverages, prevStepCoverages); - while (!checkCoveragesConvergence(metric)) { - Logger::getInstance() - .addInfo("Coverages did not converge. Repeating flux " - "calculation.") - .print(); - - prevStepCoverages = - SmartPointer>::New(*coverages); - - rtTimer.start(); - fluxes = calculateFluxes(rayTracer, useCoverages, useProcessParams); - rtTimer.finish(); - model->getSurfaceModel()->updateCoverages(fluxes, materialIds); - - coverages = model->getSurfaceModel()->getCoverages(); - metric = calculateCoverageDeltaMetric(coverages, prevStepCoverages); - if (logLevel >= 5) { - for (auto val : metric) { - covMetricFile << val << ";"; - } - covMetricFile << "\n"; - } - - Logger::getInstance() - .addTiming("Top-down flux calculation", rtTimer) - .print(); - } - } - } - - // calculate velocities from fluxes - auto velocities = model->getSurfaceModel()->calculateVelocities( - fluxes, points, materialIds); - - // prepare velocity field - model->getVelocityField()->prepare(domain, velocities, - processDuration - remainingTime); - if (model->getVelocityField()->getTranslationFieldOptions() == 2) - transField->buildKdTree(points); - - // print debug output - if (logLevel >= 4) { - if (velocities) - diskMesh->getCellData().insertNextScalarData(*velocities, - "velocities"); - if (useCoverages) { - auto coverages = model->getSurfaceModel()->getCoverages(); - mergeScalarData(diskMesh->getCellData(), coverages); - } - auto surfaceData = model->getSurfaceModel()->getSurfaceData(); - if (surfaceData) - mergeScalarData(diskMesh->getCellData(), surfaceData); - mergeScalarData(diskMesh->getCellData(), fluxes); - printDiskMesh(diskMesh, name + "_" + std::to_string(counter) + ".vtp"); - if (domain->getCellSet()) { - domain->getCellSet()->writeVTU(name + "_cellSet_" + - std::to_string(counter) + ".vtu"); - } - counter++; - } - - // apply advection callback - if (useAdvectionCallback) { - callbackTimer.start(); - bool continueProcess = model->getAdvectionCallback()->applyPreAdvect( - processDuration - remainingTime); - callbackTimer.finish(); - Logger::getInstance() - .addTiming("Advection callback pre-advect", callbackTimer) - .print(); - - if (!continueProcess) { - Logger::getInstance() - .addInfo("Process stopped early by AdvectionCallback during " - "`preAdvect`.") - .print(); - break; - } - } - - // adjust time step near end - if (remainingTime - previousTimeStep < 0.) { - advectionKernel.setAdvectionTime(remainingTime); - } - - // move coverages to LS, so they are moved with the advection step - if (useCoverages) - moveCoveragesToTopLS(translator, - model->getSurfaceModel()->getCoverages()); - advTimer.start(); - advectionKernel.apply(); - advTimer.finish(); - Logger::getInstance().addTiming("Surface advection", advTimer).print(); - - if (advectionParams.velocityOutput) { - auto lsMesh = SmartPointer>::New(); - viennals::ToMesh(domain->getLevelSets().back(), lsMesh) - .apply(); - viennals::VTKWriter( - lsMesh, "ls_velocities_" + std::to_string(lsVelCounter++) + ".vtp") - .apply(); - } - - // update the translator to retrieve the correct coverages from the LS - meshConverter.apply(); - if (useCoverages) - updateCoveragesFromAdvectedSurface( - translator, model->getSurfaceModel()->getCoverages()); - - // apply advection callback - if (useAdvectionCallback) { - callbackTimer.start(); - bool continueProcess = model->getAdvectionCallback()->applyPostAdvect( - advectionKernel.getAdvectedTime()); - callbackTimer.finish(); - Logger::getInstance() - .addTiming("Advection callback post-advect", callbackTimer) - .print(); - if (!continueProcess) { - Logger::getInstance() - .addInfo("Process stopped early by AdvectionCallback during " - "`postAdvect`.") - .print(); - break; - } - } - - previousTimeStep = advectionKernel.getAdvectedTime(); - if (previousTimeStep == std::numeric_limits::max()) { - Logger::getInstance() - .addInfo("Process halted: Surface velocities are zero across the " - "entire surface.") - .print(); - remainingTime = 0.; - break; - } - remainingTime -= previousTimeStep; - - if (Logger::getLogLevel() >= 2) { - std::stringstream stream; - stream << std::fixed << std::setprecision(4) - << "Process time: " << processDuration - remainingTime << " / " - << processDuration << " " - << units::Time::getInstance().toShortString(); - Logger::getInstance().addInfo(stream.str()).print(); - } - } - - processTime = processDuration - remainingTime; - processTimer.finish(); - - Logger::getInstance() - .addTiming("\nProcess " + name, processTimer) - .addTiming("Surface advection total time", - advTimer.totalDuration * 1e-9, - processTimer.totalDuration * 1e-9) - .print(); - if (useRayTracing) { - Logger::getInstance() - .addTiming("Top-down flux calculation total time", - rtTimer.totalDuration * 1e-9, - processTimer.totalDuration * 1e-9) - .print(); - } - if (useAdvectionCallback) { - Logger::getInstance() - .addTiming("Advection callback total time", - callbackTimer.totalDuration * 1e-9, - processTimer.totalDuration * 1e-9) - .print(); - } - model->reset(); - if (logLevel >= 5) - covMetricFile.close(); - } - - void writeParticleDataLogs(std::string fileName) { + void writeParticleDataLogs(const std::string &fileName) { std::ofstream file(fileName.c_str()); for (std::size_t i = 0; i < particleDataLogs.size(); i++) { @@ -790,22 +67,7 @@ template class Process { file.close(); } -private: - static void printDiskMesh(SmartPointer> mesh, - std::string name) { - viennals::VTKWriter(mesh, std::move(name)).apply(); - } - - static void - mergeScalarData(viennals::PointData &scalarData, - SmartPointer> dataToInsert) { - int numScalarData = dataToInsert->getScalarDataSize(); - for (int i = 0; i < numScalarData; i++) { - scalarData.insertReplaceScalarData(*dataToInsert->getScalarData(i), - dataToInsert->getScalarDataLabel(i)); - } - } - +protected: static viennaray::TracingData movePointDataToRayData( SmartPointer> pointData) { viennaray::TracingData rayData; @@ -830,120 +92,77 @@ template class Process { rayData.getVectorDataLabel(i)); } - void moveCoveragesToTopLS( - SmartPointer translator, - SmartPointer> coverages) { - auto topLS = domain->getLevelSets().back(); - for (size_t i = 0; i < coverages->getScalarDataSize(); i++) { - auto covName = coverages->getScalarDataLabel(i); - std::vector levelSetData(topLS->getNumberOfPoints(), 0); - auto cov = coverages->getScalarData(covName); - for (const auto iter : *translator.get()) { - levelSetData[iter.first] = cov->at(iter.second); - } - if (auto data = topLS->getPointData().getScalarData(covName, true); - data != nullptr) { - *data = std::move(levelSetData); - } else { - topLS->getPointData().insertNextScalarData(std::move(levelSetData), - covName); - } - } - } - - void updateCoveragesFromAdvectedSurface( - SmartPointer translator, - SmartPointer> coverages) const { - auto topLS = domain->getLevelSets().back(); - for (size_t i = 0; i < coverages->getScalarDataSize(); i++) { - auto covName = coverages->getScalarDataLabel(i); - auto levelSetData = topLS->getPointData().getScalarData(covName); - auto covData = coverages->getScalarData(covName); - covData->resize(translator->size()); - for (const auto it : *translator.get()) { - covData->at(it.second) = levelSetData->at(it.first); - } - } - } - - static std::vector calculateCoverageDeltaMetric( - SmartPointer> updated, - SmartPointer> previous) { - - assert(updated->getScalarDataSize() == previous->getScalarDataSize()); - std::vector delta(updated->getScalarDataSize(), 0.); - -#pragma omp parallel for - for (int i = 0; i < updated->getScalarDataSize(); i++) { - auto label = updated->getScalarDataLabel(i); - auto updatedData = updated->getScalarData(label); - auto previousData = previous->getScalarData(label); - for (size_t j = 0; j < updatedData->size(); j++) { - auto diff = updatedData->at(j) - previousData->at(j); - delta[i] += diff * diff; - } + bool checkInput() override { return true; } - delta[i] /= updatedData->size(); - } - - return delta; - } - - bool - checkCoveragesConvergence(const std::vector &deltaMetric) const { - for (auto val : deltaMetric) { - if (val > coverageDeltaThreshold) - return false; - } - return true; - } - - void initializeRayTracer(viennaray::Trace &tracer) const { + void initFluxEngine() override { // Map the domain boundary to the ray tracing boundaries viennaray::BoundaryCondition rayBoundaryCondition[D]; - if (rayTracingParams.ignoreFluxBoundaries) { + if (rayTracingParams_.ignoreFluxBoundaries) { for (unsigned i = 0; i < D; ++i) rayBoundaryCondition[i] = viennaray::BoundaryCondition::IGNORE; } else { for (unsigned i = 0; i < D; ++i) rayBoundaryCondition[i] = utils::convertBoundaryCondition( - domain->getGrid().getBoundaryConditions(i)); + domain_->getGrid().getBoundaryConditions(i)); } - tracer.setBoundaryConditions(rayBoundaryCondition); - tracer.setSourceDirection(rayTracingParams.sourceDirection); - tracer.setNumberOfRaysPerPoint(rayTracingParams.raysPerPoint); - tracer.setUseRandomSeeds(rayTracingParams.useRandomSeeds); - tracer.setCalculateFlux(false); - - auto source = model->getSource(); - if (source) { - tracer.setSource(source); + rayTracer.setBoundaryConditions(rayBoundaryCondition); + rayTracer.setSourceDirection(rayTracingParams_.sourceDirection); + rayTracer.setNumberOfRaysPerPoint(rayTracingParams_.raysPerPoint); + rayTracer.setUseRandomSeeds(rayTracingParams_.useRandomSeeds); + rayTracer.setCalculateFlux(false); + + if (auto source = processModel_->getSource()) { + rayTracer.setSource(source); Logger::getInstance().addInfo("Using custom source.").print(); } - auto primaryDirection = model->getPrimaryDirection(); - if (primaryDirection) { + if (auto primaryDirection = processModel_->getPrimaryDirection()) { Logger::getInstance() .addInfo("Using primary direction: " + utils::arrayToString(primaryDirection.value())) .print(); - tracer.setPrimaryDirection(primaryDirection.value()); + rayTracer.setPrimaryDirection(primaryDirection.value()); + } + + // initialize particle data logs + particleDataLogs.resize(processModel_->getParticleTypes().size()); + for (std::size_t i = 0; i < processModel_->getParticleTypes().size(); i++) { + if (int logSize = processModel_->getParticleLogSize(i); logSize > 0) { + particleDataLogs[i].data.resize(1); + particleDataLogs[i].data[0].resize(logSize); + } } } - auto calculateFluxes(viennaray::Trace &rayTracer, - const bool useCoverages, const bool useProcessParams) { + void setFluxEngineGeometry() override { + assert(diskMesh_ != nullptr); + auto points = diskMesh_->getNodes(); + auto normals = *diskMesh_->getCellData().getVectorData("Normals"); + auto materialIds = *diskMesh_->getCellData().getScalarData("MaterialIds"); + + if (rayTracingParams_.diskRadius == 0.) { + rayTracer.setGeometry(points, normals, domain_->getGrid().getGridDelta()); + } else { + rayTracer.setGeometry(points, normals, domain_->getGrid().getGridDelta(), + rayTracingParams_.diskRadius); + } + rayTracer.setMaterialIds(materialIds); + } + + SmartPointer> + calculateFluxes(const bool useCoverages, + const bool useProcessParams) override { viennaray::TracingData rayTracingData; + auto surfaceModel = this->model_->getSurfaceModel(); // move coverages to the ray tracer if (useCoverages) { - rayTracingData = - movePointDataToRayData(model->getSurfaceModel()->getCoverages()); + rayTracingData = movePointDataToRayData(surfaceModel->getCoverages()); } if (useProcessParams) { // store scalars in addition to coverages - auto processParams = model->getSurfaceModel()->getProcessParameters(); + auto processParams = surfaceModel->getProcessParameters(); NumericType numParams = processParams->getScalarData().size(); rayTracingData.setNumberOfScalarData(numParams); for (size_t i = 0; i < numParams; ++i) { @@ -955,21 +174,21 @@ template class Process { if (useCoverages || useProcessParams) rayTracer.setGlobalData(rayTracingData); - auto fluxes = runRayTracer(rayTracer); + auto fluxes = runRayTracer(); // move coverages back in the model if (useCoverages) - moveRayDataToPointData(model->getSurfaceModel()->getCoverages(), - rayTracingData); + moveRayDataToPointData(surfaceModel->getCoverages(), rayTracingData); return fluxes; } - auto runRayTracer(viennaray::Trace &rayTracer) { +private: + SmartPointer> runRayTracer() { auto fluxes = SmartPointer>::New(); unsigned particleIdx = 0; - for (auto &particle : model->getParticleTypes()) { - int dataLogSize = model->getParticleLogSize(particleIdx); + for (auto &particle : processModel_->getParticleTypes()) { + int dataLogSize = processModel_->getParticleLogSize(particleIdx); if (dataLogSize > 0) { rayTracer.getDataLog().data.resize(1); rayTracer.getDataLog().data[0].resize(dataLogSize, 0.); @@ -984,9 +203,9 @@ template class Process { auto flux = std::move(localData.getVectorData(i)); // normalize and smooth - rayTracer.normalizeFlux(flux, rayTracingParams.normalizationType); - if (rayTracingParams.smoothingNeighbors > 0) - rayTracer.smoothFlux(flux, rayTracingParams.smoothingNeighbors); + rayTracer.normalizeFlux(flux, rayTracingParams_.normalizationType); + if (rayTracingParams_.smoothingNeighbors > 0) + rayTracer.smoothFlux(flux, rayTracingParams_.smoothingNeighbors); fluxes->insertNextScalarData(std::move(flux), localData.getVectorDataLabel(i)); @@ -1001,17 +220,14 @@ template class Process { } private: - psDomainType domain; - SmartPointer> model; + // Members + using ProcessBase::domain_; + using ProcessBase::rayTracingParams_; + using ProcessBase::diskMesh_; + + SmartPointer> processModel_; + viennaray::Trace rayTracer; std::vector> particleDataLogs; - NumericType processDuration; - unsigned maxIterations = std::numeric_limits::max(); - NumericType coverageDeltaThreshold = 0.; - bool coveragesInitialized_ = false; - NumericType processTime = 0.; - - AdvectionParameters advectionParams; - RayTracingParameters rayTracingParams; }; -} // namespace viennaps +} // namespace viennaps \ No newline at end of file diff --git a/include/viennaps/psProcessBase.hpp b/include/viennaps/psProcessBase.hpp new file mode 100644 index 00000000..a7162d0c --- /dev/null +++ b/include/viennaps/psProcessBase.hpp @@ -0,0 +1,780 @@ +#pragma once + +#include "psProcessModelBase.hpp" +#include "psProcessParams.hpp" +#include "psTranslationField.hpp" +#include "psUnits.hpp" + +#include +#include +#include +#include +#include + +namespace viennaps { + +using namespace viennacore; + +/// This class server as the main process tool, applying a user- or pre-defined +/// process model to a domain. Depending on the user inputs surface advection, a +/// single callback function or a geometric advection is applied. +template class ProcessBase { + using TranslatorType = std::unordered_map; + using DomainType = SmartPointer>; + +public: + ProcessBase() = default; + ProcessBase(DomainType domain) : domain_(domain) {} + // Constructor for a process with a pre-configured process model. + ProcessBase(DomainType domain, + SmartPointer> processModel, + const NumericType duration = 0.) + : domain_(domain), model_(processModel), processDuration_(duration) {} + virtual ~ProcessBase() { assert(!covMetricFile.is_open()); } + + // Set the process domain. + void setDomain(DomainType domain) { domain_ = domain; } + + /* ----- Process parameters ----- */ + + // Set the duration of the process. + void setProcessDuration(NumericType duration) { processDuration_ = duration; } + + // Returns the duration of the recently run process. This duration can + // sometimes slightly vary from the set process duration, due to the maximum + // time step according to the CFL condition. + NumericType getProcessDuration() const { return processTime_; } + + /* ----- Ray tracing parameters ----- */ + + // Set the number of iterations to initialize the coverages. + void setMaxCoverageInitIterations(unsigned maxIt) { maxIterations_ = maxIt; } + + // Set the threshold for the coverage delta metric to reach convergence. + void setCoverageDeltaThreshold(NumericType threshold) { + coverageDeltaThreshold_ = threshold; + } + + // Set the source direction, where the rays should be traced from. + void setSourceDirection(viennaray::TraceDirection passedDirection) { + rayTracingParams_.sourceDirection = passedDirection; + } + + // Specify the number of rays to be traced for each particle throughout the + // process. The total count of rays is the product of this number and the + // number of points in the process geometry. + void setNumberOfRaysPerPoint(unsigned numRays) { + rayTracingParams_.raysPerPoint = numRays; + } + + // Enable flux smoothing. The flux at each surface point, calculated + // by the ray tracer, is averaged over the surface point neighbors. + void enableFluxSmoothing() { rayTracingParams_.smoothingNeighbors = 1; } + + // Disable flux smoothing. + void disableFluxSmoothing() { rayTracingParams_.smoothingNeighbors = 0; } + + void enableFluxBoundaries() { + rayTracingParams_.ignoreFluxBoundaries = false; + } + + // Ignore boundary conditions during the flux calculation. + void disableFluxBoundaries() { + rayTracingParams_.ignoreFluxBoundaries = true; + } + + // Enable the use of random seeds for ray tracing. This is useful to + // prevent the formation of artifacts in the flux calculation. + void enableRandomSeeds() { rayTracingParams_.useRandomSeeds = true; } + + // Disable the use of random seeds for ray tracing. + void disableRandomSeeds() { rayTracingParams_.useRandomSeeds = false; } + + void setRayTracingParameters( + const RayTracingParameters &passedRayTracingParams) { + rayTracingParams_ = passedRayTracingParams; + } + + auto &getRayTracingParameters() { return rayTracingParams_; } + + /* ----- Advection parameters ----- */ + + // Set the integration scheme for solving the level-set equation. + // Possible integration schemes are specified in + // viennals::IntegrationSchemeEnum. + void setIntegrationScheme(IntegrationScheme passedIntegrationScheme) { + advectionParams_.integrationScheme = passedIntegrationScheme; + } + + // Enable the output of the advection velocities on the level-set mesh. + void enableAdvectionVelocityOutput() { + advectionParams_.velocityOutput = true; + } + + // Disable the output of the advection velocities on the level-set mesh. + void disableAdvectionVelocityOutput() { + advectionParams_.velocityOutput = false; + } + + // Set the CFL (Courant-Friedrichs-Levy) condition to use during surface + // advection in the level-set. The CFL condition defines the maximum distance + // a surface is allowed to move in a single advection step. It MUST be below + // 0.5 to guarantee numerical stability. Defaults to 0.4999. + void setTimeStepRatio(NumericType cfl) { + advectionParams_.timeStepRatio = cfl; + } + + void setAdvectionParameters( + const AdvectionParameters &passedAdvectionParams) { + advectionParams_ = passedAdvectionParams; + } + + auto &getAdvectionParameters() { return advectionParams_; } + + /* ----- Process execution ----- */ + + // A single flux calculation is performed on the domain surface. The result is + // stored as point data on the nodes of the mesh. + SmartPointer> calculateFlux() { + + if (!checkModelAndDomain()) + return nullptr; + + model_->initialize(domain_, 0.); + const auto name = model_->getProcessName().value_or("default"); + + if (!model_->useFluxEngine()) { + Logger::getInstance() + .addWarning("Process model '" + name + "' does not use flux engine.") + .print(); + return nullptr; + } + + if (!checkInput()) + return nullptr; + + initFluxEngine(); + const auto logLevel = Logger::getLogLevel(); + + // Generate disk mesh from domain + diskMesh_ = SmartPointer>::New(); + viennals::ToDiskMesh meshGenerator(diskMesh_); + if (domain_->getMaterialMap() && + domain_->getMaterialMap()->size() == domain_->getLevelSets().size()) { + meshGenerator.setMaterialMap(domain_->getMaterialMap()->getMaterialMap()); + } + for (auto ls : domain_->getLevelSets()) { + meshGenerator.insertNextLevelSet(ls); + } + meshGenerator.apply(); + auto numPoints = diskMesh_->getNodes().size(); + const auto &materialIds = + *diskMesh_->getCellData().getScalarData("MaterialIds"); + + setFluxEngineGeometry(); + + const bool useProcessParams = + model_->getSurfaceModel()->getProcessParameters() != nullptr; + bool useCoverages = false; + + model_->getSurfaceModel()->initializeSurfaceData(numPoints); + // Check if coverages are used + if (!coveragesInitialized_) + model_->getSurfaceModel()->initializeCoverages(numPoints); + auto coverages = model_->getSurfaceModel()->getCoverages(); + + if (coverages != nullptr) { + useCoverages = true; + Logger::getInstance().addInfo("Using coverages.").print(); + + // debug output + if (logLevel >= 5) + covMetricFile.open(name + "_covMetric.txt"); + + coverageInitIterations(useProcessParams); + + if (logLevel >= 5) + covMetricFile.close(); + } + + auto fluxes = calculateFluxes(useCoverages, useProcessParams); + mergeScalarData(diskMesh_->getCellData(), fluxes); + if (auto surfaceData = model_->getSurfaceModel()->getSurfaceData()) + mergeScalarData(diskMesh_->getCellData(), surfaceData); + + return diskMesh_; + } + + // Run the process. + void apply() { + /* ---------- Check input --------- */ + if (!checkModelAndDomain()) + return; + + model_->initialize(domain_, processDuration_); + const auto name = model_->getProcessName().value_or("default"); + + if (model_->getGeometricModel()) { + model_->getGeometricModel()->setDomain(domain_); + Logger::getInstance().addInfo("Applying geometric model...").print(); + model_->getGeometricModel()->apply(); + return; + } + + if (processDuration_ == 0.) { + // apply only advection callback + if (model_->getAdvectionCallback()) { + model_->getAdvectionCallback()->setDomain(domain_); + model_->getAdvectionCallback()->applyPreAdvect(0); + } else { + Logger::getInstance() + .addWarning("No advection callback passed to Process.") + .print(); + } + return; + } + + if (!model_->getSurfaceModel()) { + Logger::getInstance() + .addWarning("No surface model passed to Process.") + .print(); + return; + } + + if (!model_->getVelocityField()) { + Logger::getInstance() + .addWarning("No velocity field passed to Process.") + .print(); + return; + } + + // check implementation specific input + if (!checkInput()) + return; + + /* ------ Process Setup ------ */ + const unsigned int logLevel = Logger::getLogLevel(); + Timer processTimer; + processTimer.start(); + + auto surfaceModel = model_->getSurfaceModel(); + auto velocityField = model_->getVelocityField(); + auto advectionCallback = model_->getAdvectionCallback(); + double remainingTime = processDuration_; + const NumericType gridDelta = domain_->getGrid().getGridDelta(); + + diskMesh_ = SmartPointer>::New(); + auto translator = SmartPointer::New(); + viennals::ToDiskMesh meshGenerator(diskMesh_); + meshGenerator.setTranslator(translator); + if (domain_->getMaterialMap() && + domain_->getMaterialMap()->size() == domain_->getLevelSets().size()) { + meshGenerator.setMaterialMap(domain_->getMaterialMap()->getMaterialMap()); + } + + translationField_ = SmartPointer>::New( + velocityField, domain_->getMaterialMap()); + translationField_->setTranslator(translator); + + viennals::Advect advectionKernel; + advectionKernel.setVelocityField(translationField_); + advectionKernel.setIntegrationScheme(advectionParams_.integrationScheme); + advectionKernel.setTimeStepRatio(advectionParams_.timeStepRatio); + advectionKernel.setSaveAdvectionVelocities(advectionParams_.velocityOutput); + advectionKernel.setDissipationAlpha(advectionParams_.dissipationAlpha); + advectionKernel.setIgnoreVoids(advectionParams_.ignoreVoids); + advectionKernel.setCheckDissipation(advectionParams_.checkDissipation); + // normals vectors are only necessary for analytical velocity fields + if (velocityField->getTranslationFieldOptions() > 0) + advectionKernel.setCalculateNormalVectors(false); + + for (auto dom : domain_->getLevelSets()) { + meshGenerator.insertNextLevelSet(dom); + advectionKernel.insertNextLevelSet(dom); + } + + // Check if the process model uses ray tracing + const bool useFluxEngine = model_->useFluxEngine(); + if (useFluxEngine) { + initFluxEngine(); + } + + // Determine whether advection callback is used + const bool useAdvectionCallback = advectionCallback != nullptr; + if (useAdvectionCallback) { + advectionCallback->setDomain(domain_); + } + + // Determine whether there are process parameters used in ray tracing + surfaceModel->initializeProcessParameters(); + const bool useProcessParams = + surfaceModel->getProcessParameters() != nullptr; + + if (useProcessParams) + Logger::getInstance().addInfo("Using process parameters.").print(); + if (useAdvectionCallback) + Logger::getInstance().addInfo("Using advection callback.").print(); + + bool useCoverages = false; + + // Initialize coverages + meshGenerator.apply(); + auto numPoints = diskMesh_->getNodes().size(); + surfaceModel->initializeSurfaceData(numPoints); + if (!coveragesInitialized_) + surfaceModel->initializeCoverages(numPoints); + auto coverages = surfaceModel->getCoverages(); + if (coverages != nullptr) { + Logger::getInstance().addInfo("Using coverages.").print(); + useCoverages = true; + + // debug output + if (logLevel >= 5) + covMetricFile.open(name + "_covMetric.txt"); + + setFluxEngineGeometry(); + coverageInitIterations(useProcessParams); + } + + double previousTimeStep = 0.; + size_t counter = 0; + unsigned lsVelCounter = 0; + Timer rtTimer; + Timer callbackTimer; + Timer advTimer; + while (remainingTime > 0.) { + // We need additional signal handling when running the C++ code from the + // Python bindings to allow interrupts in the Python scripts +#ifdef VIENNAPS_PYTHON_BUILD + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); +#endif + // Expand LS based on the integration scheme + advectionKernel.prepareLS(); + model_->initialize(domain_, remainingTime); + + auto fluxes = SmartPointer>::New(); + meshGenerator.apply(); + auto materialIds = *diskMesh_->getCellData().getScalarData("MaterialIds"); + auto points = diskMesh_->getNodes(); + + // rate calculation by implementation specific flux engine + if (useFluxEngine) { + rtTimer.start(); + + setFluxEngineGeometry(); + fluxes = calculateFluxes(useCoverages, useProcessParams); + + rtTimer.finish(); + Logger::getInstance().addTiming("Flux calculation", rtTimer).print(); + } + + // update coverages and calculate coverage delta metric + if (useCoverages) { + coverages = surfaceModel->getCoverages(); + auto prevStepCoverages = + SmartPointer>::New(*coverages); + + // update coverages in the model + surfaceModel->updateCoverages(fluxes, materialIds); + + if (coverageDeltaThreshold_ > 0) { + auto metric = + this->calculateCoverageDeltaMetric(coverages, prevStepCoverages); + while (!this->checkCoveragesConvergence(metric)) { + Logger::getInstance() + .addInfo("Coverages did not converge. Repeating flux " + "calculation.") + .print(); + + prevStepCoverages = + SmartPointer>::New(*coverages); + + rtTimer.start(); + fluxes = calculateFluxes(useCoverages, useProcessParams); + rtTimer.finish(); + surfaceModel->updateCoverages(fluxes, materialIds); + + coverages = surfaceModel->getCoverages(); + metric = this->calculateCoverageDeltaMetric(coverages, + prevStepCoverages); + + Logger::getInstance() + .addTiming("Top-down flux calculation", rtTimer) + .print(); + } + } + } + + // calculate velocities from fluxes + auto velocities = + surfaceModel->calculateVelocities(fluxes, points, materialIds); + + // prepare velocity field + velocityField->prepare(domain_, velocities, + processDuration_ - remainingTime); + if (velocityField->getTranslationFieldOptions() == 2) + translationField_->buildKdTree(points); + + // print debug output + if (logLevel >= 4) { + if (velocities) + diskMesh_->getCellData().insertNextScalarData(*velocities, + "velocities"); + if (useCoverages) { + mergeScalarData(diskMesh_->getCellData(), + surfaceModel->getCoverages()); + } + if (auto surfaceData = surfaceModel->getSurfaceData()) + mergeScalarData(diskMesh_->getCellData(), surfaceData); + mergeScalarData(diskMesh_->getCellData(), fluxes); + printDiskMesh(diskMesh_, name + "_" + std::to_string(counter) + ".vtp"); + if (domain_->getCellSet()) { + domain_->getCellSet()->writeVTU(name + "_cellSet_" + + std::to_string(counter) + ".vtu"); + } + counter++; + } + + // apply advection callback + if (useAdvectionCallback) { + callbackTimer.start(); + bool continueProcess = + advectionCallback->applyPreAdvect(processDuration_ - remainingTime); + callbackTimer.finish(); + Logger::getInstance() + .addTiming("Advection callback pre-advect", callbackTimer) + .print(); + + if (!continueProcess) { + Logger::getInstance() + .addInfo("Process stopped early by AdvectionCallback during " + "`preAdvect`.") + .print(); + break; + } + } + + // adjust time step near end + if (remainingTime - previousTimeStep < 0.) { + advectionKernel.setAdvectionTime(remainingTime); + } + + // move coverages to LS, so they are moved with the advection step + if (useCoverages) + this->moveCoveragesToTopLS(translator, surfaceModel->getCoverages()); + advTimer.start(); + advectionKernel.apply(); + advTimer.finish(); + Logger::getInstance().addTiming("Surface advection", advTimer).print(); + + if (advectionParams_.velocityOutput) { + auto lsMesh = SmartPointer>::New(); + viennals::ToMesh(domain_->getLevelSets().back(), lsMesh) + .apply(); + viennals::VTKWriter( + lsMesh, "ls_velocities_" + std::to_string(lsVelCounter++) + ".vtp") + .apply(); + } + + // update the translator to retrieve the correct coverages from the LS + meshGenerator.apply(); + if (useCoverages) + this->updateCoveragesFromAdvectedSurface(translator, + surfaceModel->getCoverages()); + + // apply advection callback + if (useAdvectionCallback) { + callbackTimer.start(); + bool continueProcess = advectionCallback->applyPostAdvect( + advectionKernel.getAdvectedTime()); + callbackTimer.finish(); + Logger::getInstance() + .addTiming("Advection callback post-advect", callbackTimer) + .print(); + if (!continueProcess) { + Logger::getInstance() + .addInfo("Process stopped early by AdvectionCallback during " + "`postAdvect`.") + .print(); + break; + } + } + + previousTimeStep = advectionKernel.getAdvectedTime(); + if (previousTimeStep == std::numeric_limits::max()) { + Logger::getInstance() + .addInfo("Process halted: Surface velocities are zero across the " + "entire surface.") + .print(); + remainingTime = 0.; + break; + } + remainingTime -= previousTimeStep; + + if (Logger::getLogLevel() >= 2) { + std::stringstream stream; + stream << std::fixed << std::setprecision(4) + << "Process time: " << processDuration_ - remainingTime << " / " + << processDuration_ << " " + << units::Time::getInstance().toShortString(); + Logger::getInstance().addInfo(stream.str()).print(); + } + } + + processTime_ = processDuration_ - remainingTime; + processTimer.finish(); + + Logger::getInstance() + .addTiming("\nProcess " + name, processTimer) + .addTiming("Surface advection total time", + advTimer.totalDuration * 1e-9, + processTimer.totalDuration * 1e-9) + .print(); + if (useFluxEngine) { + Logger::getInstance() + .addTiming("Flux calculation total time", + rtTimer.totalDuration * 1e-9, + processTimer.totalDuration * 1e-9) + .print(); + } + if (useAdvectionCallback) { + Logger::getInstance() + .addTiming("Advection callback total time", + callbackTimer.totalDuration * 1e-9, + processTimer.totalDuration * 1e-9) + .print(); + } + model_->reset(); + if (useCoverages && logLevel >= 5) + covMetricFile.close(); + } + +protected: + bool checkModelAndDomain() const { + if (!domain_) { + Logger::getInstance().addWarning("No domain passed to Process.").print(); + return false; + } + + if (domain_->getLevelSets().empty()) { + Logger::getInstance().addWarning("No level sets in domain.").print(); + return false; + } + + if (!model_) { + Logger::getInstance() + .addWarning("No process model passed to Process.") + .print(); + return false; + } + + return true; + } + + static void printDiskMesh(SmartPointer> mesh, + std::string name) { + viennals::VTKWriter(mesh, std::move(name)).apply(); + } + + static void + mergeScalarData(viennals::PointData &scalarData, + SmartPointer> dataToInsert) { + int numScalarData = dataToInsert->getScalarDataSize(); + for (int i = 0; i < numScalarData; i++) { + scalarData.insertReplaceScalarData(*dataToInsert->getScalarData(i), + dataToInsert->getScalarDataLabel(i)); + } + } + + void moveCoveragesToTopLS( + SmartPointer const &translator, + SmartPointer> coverages) { + auto topLS = domain_->getLevelSets().back(); + for (size_t i = 0; i < coverages->getScalarDataSize(); i++) { + auto covName = coverages->getScalarDataLabel(i); + std::vector levelSetData(topLS->getNumberOfPoints(), 0); + auto cov = coverages->getScalarData(covName); + for (const auto iter : *translator.get()) { + levelSetData[iter.first] = cov->at(iter.second); + } + if (auto data = topLS->getPointData().getScalarData(covName, true); + data != nullptr) { + *data = std::move(levelSetData); + } else { + topLS->getPointData().insertNextScalarData(std::move(levelSetData), + covName); + } + } + } + + void updateCoveragesFromAdvectedSurface( + SmartPointer const &translator, + SmartPointer> coverages) const { + auto topLS = domain_->getLevelSets().back(); + for (size_t i = 0; i < coverages->getScalarDataSize(); i++) { + auto covName = coverages->getScalarDataLabel(i); + auto levelSetData = topLS->getPointData().getScalarData(covName); + auto covData = coverages->getScalarData(covName); + covData->resize(translator->size()); + for (const auto it : *translator.get()) { + covData->at(it.second) = levelSetData->at(it.first); + } + } + } + + std::vector calculateCoverageDeltaMetric( + SmartPointer> updated, + SmartPointer> previous) { + + assert(updated->getScalarDataSize() == previous->getScalarDataSize()); + std::vector delta(updated->getScalarDataSize(), 0.); + +#pragma omp parallel for + for (int i = 0; i < updated->getScalarDataSize(); i++) { + auto label = updated->getScalarDataLabel(i); + auto updatedData = updated->getScalarData(label); + auto previousData = previous->getScalarData(label); + for (size_t j = 0; j < updatedData->size(); j++) { + auto diff = updatedData->at(j) - previousData->at(j); + delta[i] += diff * diff; + } + + delta[i] /= updatedData->size(); + } + + logMetric(delta); + return delta; + } + + bool + checkCoveragesConvergence(const std::vector &deltaMetric) const { + for (auto val : deltaMetric) { + if (val > coverageDeltaThreshold_) + return false; + } + return true; + } + + void coverageInitIterations(const bool useProcessParams) { + auto coverages = model_->getSurfaceModel()->getCoverages(); + const auto name = model_->getProcessName().value_or("default"); + const auto logLevel = Logger::getLogLevel(); + + if (maxIterations_ == std::numeric_limits::max() && + coverageDeltaThreshold_ == 0.) { + maxIterations_ = 10; + Logger::getInstance() + .addWarning("No coverage initialization parameters set. Using " + + std::to_string(maxIterations_) + + " initialization iterations.") + .print(); + } + + if (!coveragesInitialized_) { + Timer timer; + timer.start(); + Logger::getInstance().addInfo("Initializing coverages ... ").print(); + + for (unsigned iteration = 0; iteration < maxIterations_; ++iteration) { + // We need additional signal handling when running the C++ code from + // the Python bindings to allow interrupts in the Python scripts +#ifdef VIENNAPS_PYTHON_BUILD + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); +#endif + // save current coverages to compare with the new ones + auto prevStepCoverages = + SmartPointer>::New(*coverages); + + auto fluxes = calculateFluxes(true, useProcessParams); + + // update coverages in the model + assert(diskMesh_ != nullptr); + auto const &materialIds = + *diskMesh_->getCellData().getScalarData("MaterialIds"); + model_->getSurfaceModel()->updateCoverages(fluxes, materialIds); + + // metric to check for convergence + coverages = model_->getSurfaceModel()->getCoverages(); + auto metric = + calculateCoverageDeltaMetric(coverages, prevStepCoverages); + + if (logLevel >= 3) { + mergeScalarData(diskMesh_->getCellData(), coverages); + mergeScalarData(diskMesh_->getCellData(), fluxes); + if (auto surfaceData = model_->getSurfaceModel()->getSurfaceData()) + mergeScalarData(diskMesh_->getCellData(), surfaceData); + printDiskMesh(diskMesh_, name + "_covInit_" + + std::to_string(iteration) + ".vtp"); + + Logger::getInstance() + .addInfo("Iteration: " + std::to_string(iteration + 1)) + .print(); + + std::stringstream stream; + stream << std::setprecision(4) << std::fixed; + stream << "Coverage delta metric: "; + for (int i = 0; i < coverages->getScalarDataSize(); i++) { + stream << coverages->getScalarDataLabel(i) << ": " << metric[i] + << "\t"; + } + Logger::getInstance().addInfo(stream.str()).print(); + } + + if (checkCoveragesConvergence(metric)) { + Logger::getInstance() + .addInfo("Coverages converged after " + + std::to_string(iteration + 1) + " iterations.") + .print(); + break; + } + } // end coverage initialization loop + coveragesInitialized_ = true; + + timer.finish(); + Logger::getInstance().addTiming("Coverage initialization", timer).print(); + } else { + Logger::getInstance().addInfo("Coverages already initialized.").print(); + } + } + + void logMetric(const std::vector &metric) { + if (Logger::getLogLevel() < 5) + return; + + assert(covMetricFile.is_open()); + for (auto val : metric) { + covMetricFile << val << ";"; + } + covMetricFile << "\n"; + } + + // Implementation specific functions + virtual bool checkInput() = 0; + + virtual void initFluxEngine() = 0; + + virtual void setFluxEngineGeometry() = 0; + + virtual SmartPointer> + calculateFluxes(const bool useCoverages, const bool useProcessParams) = 0; + +protected: + DomainType domain_; + SmartPointer> model_; + NumericType processDuration_; + NumericType processTime_ = 0.; + unsigned maxIterations_ = std::numeric_limits::max(); + NumericType coverageDeltaThreshold_ = 0.; + bool coveragesInitialized_ = false; + + AdvectionParameters advectionParams_; + RayTracingParameters rayTracingParams_; + + SmartPointer> diskMesh_ = nullptr; + SmartPointer> translationField_ = nullptr; + std::ofstream covMetricFile; +}; + +} // namespace viennaps \ No newline at end of file diff --git a/include/viennaps/psProcessModel.hpp b/include/viennaps/psProcessModel.hpp index 9e7cf773..1a3fa002 100644 --- a/include/viennaps/psProcessModel.hpp +++ b/include/viennaps/psProcessModel.hpp @@ -1,52 +1,31 @@ #pragma once -#include "psAdvectionCallback.hpp" -#include "psDomain.hpp" -#include "psGeometricModel.hpp" -#include "psSurfaceModel.hpp" -#include "psVelocityField.hpp" +#include "psProcessModelBase.hpp" #include -#include #include #include -#include - namespace viennaps { using namespace viennacore; /// The process model combines all models (particle types, surface model, /// geometric model, advection callback) -template class ProcessModel { +template +class ProcessModel : public ProcessModelBase { protected: std::vector>> particles; SmartPointer> source = nullptr; std::vector particleLogSize; - SmartPointer> surfaceModel = nullptr; - SmartPointer> advectionCallback = nullptr; - SmartPointer> geometricModel = nullptr; - SmartPointer> velocityField = nullptr; - std::optional processName = std::nullopt; std::optional> primaryDirection = std::nullopt; public: - virtual ~ProcessModel() = default; - - virtual void initialize(SmartPointer> domain, - const NumericType processDuration) {} - virtual void reset() {} - auto &getParticleTypes() { return particles; } - auto getSurfaceModel() const { return surfaceModel; } - auto getAdvectionCallback() const { return advectionCallback; } - auto getGeometricModel() const { return geometricModel; } - auto getVelocityField() const { return velocityField; } auto getSource() { return source; } - auto getProcessName() const { return processName; } + bool useFluxEngine() override { return particles.size() > 0; } /// Set a primary direction for the source distribution (tilted distribution). virtual std::optional> @@ -58,8 +37,6 @@ template class ProcessModel { return particleLogSize[particleIdx]; } - void setProcessName(std::string name) { processName = std::move(name); } - virtual void setPrimaryDirection(const std::array passedPrimaryDirection) { primaryDirection = Normalize(passedPrimaryDirection); @@ -77,26 +54,6 @@ template class ProcessModel { void setSource(SmartPointer> passedSource) { source = passedSource; } - - void - setSurfaceModel(SmartPointer> passedSurfaceModel) { - surfaceModel = passedSurfaceModel; - } - - void setAdvectionCallback( - SmartPointer> passedAdvectionCallback) { - advectionCallback = passedAdvectionCallback; - } - - void setGeometricModel( - SmartPointer> passedGeometricModel) { - geometricModel = passedGeometricModel; - } - - void setVelocityField( - SmartPointer> passedVelocityField) { - velocityField = passedVelocityField; - } }; } // namespace viennaps diff --git a/include/viennaps/psProcessModelBase.hpp b/include/viennaps/psProcessModelBase.hpp new file mode 100644 index 00000000..7d49c86e --- /dev/null +++ b/include/viennaps/psProcessModelBase.hpp @@ -0,0 +1,60 @@ +#pragma once + +#include "psAdvectionCallback.hpp" +#include "psDomain.hpp" +#include "psGeometricModel.hpp" +#include "psSurfaceModel.hpp" +#include "psVelocityField.hpp" + +namespace viennaps { + +using namespace viennacore; + +/// The process model combines all models (particle types, surface model, +/// geometric model, advection callback) +template class ProcessModelBase { +protected: + SmartPointer> surfaceModel = nullptr; + SmartPointer> advectionCallback = nullptr; + SmartPointer> geometricModel = nullptr; + SmartPointer> velocityField = nullptr; + std::optional processName = std::nullopt; + +public: + virtual ~ProcessModelBase() = default; + + virtual void initialize(SmartPointer> domain, + const NumericType processDuration) {} + virtual void reset() {} + virtual bool useFluxEngine() { return false; } + + auto getSurfaceModel() const { return surfaceModel; } + auto getAdvectionCallback() const { return advectionCallback; } + auto getGeometricModel() const { return geometricModel; } + auto getVelocityField() const { return velocityField; } + auto getProcessName() const { return processName; } + + void setProcessName(const std::string &name) { processName = name; } + + void + setSurfaceModel(SmartPointer> passedSurfaceModel) { + surfaceModel = passedSurfaceModel; + } + + void setAdvectionCallback( + SmartPointer> passedAdvectionCallback) { + advectionCallback = passedAdvectionCallback; + } + + void setGeometricModel( + SmartPointer> passedGeometricModel) { + geometricModel = passedGeometricModel; + } + + void setVelocityField( + SmartPointer> passedVelocityField) { + velocityField = passedVelocityField; + } +}; + +} // namespace viennaps diff --git a/python/pyWrap.cpp b/python/pyWrap.cpp index cffb6275..b5e1abcb 100644 --- a/python/pyWrap.cpp +++ b/python/pyWrap.cpp @@ -341,10 +341,11 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) { .def_static("getInstance", &Logger::getInstance, pybind11::return_value_policy::reference) .def("addDebug", &Logger::addDebug) - .def("addTiming", - (Logger & (Logger::*)(std::string, double)) & Logger::addTiming) - .def("addTiming", (Logger & (Logger::*)(std::string, double, double)) & + .def("addTiming", (Logger & (Logger::*)(const std::string &, double)) & Logger::addTiming) + .def("addTiming", + (Logger & (Logger::*)(const std::string &, double, double)) & + Logger::addTiming) .def("addInfo", &Logger::addInfo) .def("addWarning", &Logger::addWarning) .def("addError", &Logger::addError, pybind11::arg("s"), diff --git a/tests/particleSampling/CMakeLists.txt b/tests/particleSampling/CMakeLists.txt new file mode 100644 index 00000000..a2a34b76 --- /dev/null +++ b/tests/particleSampling/CMakeLists.txt @@ -0,0 +1,7 @@ +project(particleSampling LANGUAGES CXX) + +add_executable(${PROJECT_NAME} "${PROJECT_NAME}.cpp") +target_link_libraries(${PROJECT_NAME} PRIVATE ViennaPS) + +add_dependencies(ViennaPS_Tests ${PROJECT_NAME}) +add_test(NAME ${PROJECT_NAME} COMMAND $) diff --git a/tests/particleSampling/particleSampling.cpp b/tests/particleSampling/particleSampling.cpp new file mode 100644 index 00000000..a954d8ba --- /dev/null +++ b/tests/particleSampling/particleSampling.cpp @@ -0,0 +1,209 @@ +#include +#include + +using namespace viennacore; + +template struct UnivariateDistribution { + std::vector pdf; + std::vector support; // (x-values) +}; + +template struct BivariateDistribution { + std::vector> pdf; // 2D grid of values + std::vector support_x; // (x-values) + std::vector support_y; // (y-values) +}; + +template +class UnivariateDistributionParticle final + : public viennaray::Particle, + NumericType> { + NumericType energy_; + Sampling + directionSampling_; // uni-variate (1D) sampling instance + // using inverse-transform sampling + Sampling + energySampling_; // uni-variate (1D) sampling instance + // using Alias sampling + +public: + UnivariateDistributionParticle( + const UnivariateDistribution &directionDist, + const UnivariateDistribution &energyDist) { + // initialize sampling instances with custom distribution + directionSampling_.setPDF(directionDist.pdf, directionDist.support); + energySampling_.setPDF(energyDist.pdf, energyDist.support); + } + + // This new initialize function can be overridden to generate a custom + // direction for the particle + Vec3D initNewWithDirection(RNG &rngState) override { + + ////// Custom Energy Sampling + auto energySample = energySampling_.sample(rngState); + // the returned sample is a 1D array + energy_ = energySample[0]; + + ////// Custom Direction Sampling + Vec3D direction = {0., 0., 0.}; + auto directionSample = directionSampling_.sample(rngState); + // the returned sample is a 1D array + NumericType theta = directionSample[0]; + if constexpr (D == 2) { + direction[0] = std::sin(theta); + direction[1] = -std::cos(theta); + } else { + std::uniform_real_distribution uniform_dist(0, 2 * M_PI); + auto phi = uniform_dist(rngState); + direction[0] = std::sin(theta) * std::cos(phi); + direction[1] = std::sin(theta) * std::sin(phi); + direction[2] = -std::cos(theta); + } + + std::cout << "Energy is " << energy_ << std::endl; + std::cout << "Direction is " << direction << std::endl; + + // THE RETURNED DIRECTION HAS TO BE NORMALIZED + // if this returns a 0-vector the default direction from + // the power cosine distribution is used + return direction; + } + + // override all other functions as usual +}; + +template +class BivariateDistributionParticle final + : public viennaray::Particle, + NumericType> { + NumericType energy_; + Sampling + directionNEnergySampling_; // bivariate (2D) sampling instance using + // accept-reject sampling + +public: + explicit BivariateDistributionParticle( + const BivariateDistribution &angleEnergyDistribution) { + // initialize sampling instances with custom distribution + directionNEnergySampling_.setPDF(angleEnergyDistribution.pdf, + angleEnergyDistribution.support_x, + angleEnergyDistribution.support_y); + } + + // This new initialize function can be overridden to generate a custom + // direction for the particle + Vec3D initNewWithDirection(RNG &rngState) override { + + // Generate sample + auto sample = directionNEnergySampling_.sample(rngState); + // sample is 2D array + + NumericType theta = sample[0]; + energy_ = sample[1]; + + Vec3D direction = {0., 0., 0.}; + if constexpr (D == 2) { + direction[0] = std::sin(theta); + direction[1] = -std::cos(theta); + } else { + std::uniform_real_distribution uniform_dist(0, 2 * M_PI); + auto phi = uniform_dist(rngState); + direction[0] = std::sin(theta) * std::cos(phi); + direction[1] = std::sin(theta) * std::sin(phi); + direction[2] = -std::cos(theta); + } + + std::cout << "Energy is " << energy_ << std::endl; + std::cout << "Direction is " << direction << std::endl; + + // THE RETURNED DIRECTION HAS TO BE NORMALIZED + // if this returns a 0-vector the default direction from + // the power cosine distribution is used + return direction; + } + + // override all other functions as usual +}; + +int main() { + using NumericType = double; + constexpr int D = 2; + Logger::setLogLevel(LogLevel::DEBUG); + + constexpr int N = 1000; // pdf resolution + + // ----- Setup with 2 uni-variate distributions ----- + { + UnivariateDistribution angularDistribution; + angularDistribution.pdf.resize(N); + angularDistribution.support.resize(N); + for (int i = 0; i < N; ++i) { + const NumericType theta = i * M_PI_2 / N; + angularDistribution.support[i] = theta; + angularDistribution.pdf[i] = + std::exp(-theta * theta / 0.5) * std::sin(theta); + } + + constexpr NumericType maxEnergy = 200.; + constexpr NumericType meanEnergy = 100.; + UnivariateDistribution energyDistribution; + energyDistribution.pdf.resize(N); + energyDistribution.support.resize(N); + for (int i = 0; i < N; ++i) { + const NumericType energy = i * maxEnergy / N; + energyDistribution.support[i] = energy; + energyDistribution.pdf[i] = + std::exp(-(energy - meanEnergy) * (energy - meanEnergy) / 20); + } + + auto testParticle = + std::make_unique>( + angularDistribution, energyDistribution); + // check if particle can be copied correctly with sampling instances + auto testCopy = testParticle->clone(); + + RNG rngState(236521); + for (int i = 0; i < 10; ++i) { + auto direction = testCopy->initNewWithDirection(rngState); + assert(IsNormalized(direction)); + } + } + + // ----- Setup with 1 bi-variate distributions ----- + { + BivariateDistribution angularEnergyDistribution; + angularEnergyDistribution.support_x.resize(N); + angularEnergyDistribution.support_y.resize(N); + angularEnergyDistribution.pdf.resize(N); + + constexpr NumericType maxEnergy = 200.; + constexpr NumericType meanEnergy = 100.; + for (int i = 0; i < N; ++i) { + const NumericType theta = i * M_PI_2 / N; + NumericType energy = i * maxEnergy / N; + angularEnergyDistribution.support_x[i] = theta; + angularEnergyDistribution.support_y[i] = energy; + + angularEnergyDistribution.pdf[i].resize(N); + for (int j = 0; j < N; ++j) { + energy = j * maxEnergy / N; + // some made-up PDF + const NumericType pdfVal = + std::exp(-theta * theta / 0.5) * std::sin(theta) * + std::exp(-(energy - meanEnergy) * (energy - meanEnergy) / 20); + angularEnergyDistribution.pdf[i][j] = pdfVal; + } + } + + auto testParticle = + std::make_unique>( + angularEnergyDistribution); + auto testCopy = testParticle->clone(); + + RNG rngState(235123612); + for (int i = 0; i < 10; ++i) { + auto direction = testCopy->initNewWithDirection(rngState); + assert(IsNormalized(direction)); + } + } +} \ No newline at end of file