From b0dee567f831a95f67d6bbcb3d6fc6957c73a7e8 Mon Sep 17 00:00:00 2001 From: Tobias Reiter <44025882+tobre1@users.noreply.github.com> Date: Mon, 24 Feb 2025 13:29:31 +0100 Subject: [PATCH] Add threshold value for coverage delta metric (#113) * Refactor coverge initialization and test for convergence * Add info message when repeating flux calculation * Add timer for additional flux calculation --- examples/holeEtching/holeEtching.cpp | 2 +- examples/holeEtching/testFluxes.py | 20 +- .../viennaps/models/psFluorocarbonEtching.hpp | 1 - include/viennaps/models/psSF6O2Etching.hpp | 1 - include/viennaps/models/psTEOSDeposition.hpp | 1 - include/viennaps/psProcess.hpp | 366 ++++++++++-------- pyproject.toml | 2 +- python/pyWrap.cpp | 4 + python/viennaps2d/viennaps2d.pyi | 1 + python/viennaps3d/viennaps3d.pyi | 1 + 10 files changed, 219 insertions(+), 180 deletions(-) diff --git a/examples/holeEtching/holeEtching.cpp b/examples/holeEtching/holeEtching.cpp index f8a4cfad..7693f11d 100644 --- a/examples/holeEtching/holeEtching.cpp +++ b/examples/holeEtching/holeEtching.cpp @@ -52,7 +52,7 @@ int main(int argc, char *argv[]) { Process process; process.setDomain(geometry); process.setProcessModel(model); - process.setMaxCoverageInitIterations(20); + process.setCoverageDeltaThreshold(1e-4); process.setNumberOfRaysPerPoint(params.get("raysPerPoint")); process.setProcessDuration(params.get("processTime")); process.setIntegrationScheme( diff --git a/examples/holeEtching/testFluxes.py b/examples/holeEtching/testFluxes.py index 9b37b1f0..eb13d5a4 100644 --- a/examples/holeEtching/testFluxes.py +++ b/examples/holeEtching/testFluxes.py @@ -23,9 +23,8 @@ vps.Time.setUnit("min") # hole geometry parameters -gridDelta = 0.025 # um -xExtent = 1.0 -yExtent = 1.0 +gridDelta = 0.03 # um +extent = 3.0 holeRadius = 0.175 maskHeight = 1.2 taperAngle = 1.193 @@ -62,22 +61,23 @@ # geometry setup, all units in um geometry = vps.Domain( - gridDelta=params["gridDelta"], - xExtent=params["xExtent"], - yExtent=params["yExtent"], + gridDelta=gridDelta, + xExtent=extent, + yExtent=extent, ) vps.MakeHole( domain=geometry, - holeRadius=params["holeRadius"], + holeRadius=holeRadius, holeDepth=0.0, - maskHeight=params["maskHeight"], - maskTaperAngle=params["taperAngle"], - holeShape=vps.HoleShape.Half, + maskHeight=maskHeight, + maskTaperAngle=taperAngle, + # holeShape=vps.HoleShape.Half, ).apply() process = vps.Process() process.setDomain(geometry) process.setMaxCoverageInitIterations(20) + process.setCoverageDeltaThreshold(1e-4) process.setProcessDuration(processDuration) process.setIntegrationScheme(integrationScheme) process.setNumberOfRaysPerPoint(numberOfRaysPerPoint) diff --git a/include/viennaps/models/psFluorocarbonEtching.hpp b/include/viennaps/models/psFluorocarbonEtching.hpp index e6a12a8e..791a093a 100644 --- a/include/viennaps/models/psFluorocarbonEtching.hpp +++ b/include/viennaps/models/psFluorocarbonEtching.hpp @@ -169,7 +169,6 @@ class FluorocarbonSurfaceModel : public SurfaceModel { calculateVelocities(SmartPointer> fluxes, const std::vector> &coordinates, const std::vector &materialIds) override { - updateCoverages(fluxes, materialIds); const auto numPoints = materialIds.size(); std::vector etchRate(numPoints, 0.); diff --git a/include/viennaps/models/psSF6O2Etching.hpp b/include/viennaps/models/psSF6O2Etching.hpp index cf9dec3e..f720634c 100644 --- a/include/viennaps/models/psSF6O2Etching.hpp +++ b/include/viennaps/models/psSF6O2Etching.hpp @@ -58,7 +58,6 @@ class SF6O2SurfaceModel : public SurfaceModel { calculateVelocities(SmartPointer> fluxes, const std::vector> &coordinates, const std::vector &materialIds) override { - updateCoverages(fluxes, materialIds); const auto numPoints = fluxes->getScalarData(0)->size(); std::vector etchRate(numPoints, 0.); diff --git a/include/viennaps/models/psTEOSDeposition.hpp b/include/viennaps/models/psTEOSDeposition.hpp index 6f5955fa..cef00f15 100644 --- a/include/viennaps/models/psTEOSDeposition.hpp +++ b/include/viennaps/models/psTEOSDeposition.hpp @@ -24,7 +24,6 @@ class SingleTEOSSurfaceModel : public SurfaceModel { SmartPointer> rates, const std::vector> &coordinates, const std::vector &materialIDs) override { - updateCoverages(rates, materialIDs); // define the surface reaction here auto particleFlux = rates->getScalarData("particleFlux"); std::vector velocity(particleFlux->size(), 0.); diff --git a/include/viennaps/psProcess.hpp b/include/viennaps/psProcess.hpp index c6e6e10c..3578a89a 100644 --- a/include/viennaps/psProcess.hpp +++ b/include/viennaps/psProcess.hpp @@ -65,6 +65,11 @@ template class Process { // 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) { + coverageDeltaTreshold = treshold; + } + // Set the source direction, where the rays should be traced from. void setSourceDirection(viennaray::TraceDirection passedDirection) { rayTracingParams.sourceDirection = passedDirection; @@ -350,6 +355,17 @@ template class Process { Timer timer; useCoverages = true; Logger::getInstance().addInfo("Using coverages.").print(); + + if (maxIterations == std::numeric_limits::max() && + coverageDeltaTreshold == 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"); @@ -357,6 +373,7 @@ template class Process { if (!coveragesInitialized_) { timer.start(); Logger::getInstance().addInfo("Initializing coverages ... ").print(); + auto points = diskMesh->getNodes(); auto normals = *diskMesh->getCellData().getVectorData("Normals"); auto materialIds = @@ -369,10 +386,11 @@ template class Process { } rayTracer.setMaterialIds(materialIds); - for (size_t iterations = 0; iterations < maxIterations; iterations++) { + 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 + // the Python bindings to allow interrupts in the Python scripts #ifdef VIENNAPS_PYTHON_BUILD if (PyErr_CheckSignals() != 0) throw pybind11::error_already_set(); @@ -380,90 +398,53 @@ template class Process { // save current coverages to compare with the new ones auto prevStepCoverages = SmartPointer>::New(*coverages); - // move coverages to the ray tracer - viennaray::TracingData rayTraceCoverages = - movePointDataToRayData(coverages); - if (useProcessParams) { - // store scalars in addition to coverages - auto processParams = - model->getSurfaceModel()->getProcessParameters(); - NumericType numParams = processParams->getScalarData().size(); - rayTraceCoverages.setNumberOfScalarData(numParams); - for (size_t i = 0; i < numParams; ++i) { - rayTraceCoverages.setScalarData( - i, processParams->getScalarData(i), - processParams->getScalarDataLabel(i)); - } - } - rayTracer.setGlobalData(rayTraceCoverages); - auto rates = SmartPointer>::New(); + auto fluxes = + calculateFluxes(rayTracer, useCoverages, useProcessParams); - std::size_t particleIdx = 0; - for (auto &particle : model->getParticleTypes()) { - int dataLogSize = model->getParticleLogSize(particleIdx); - if (dataLogSize > 0) { - rayTracer.getDataLog().data.resize(1); - rayTracer.getDataLog().data[0].resize(dataLogSize, 0.); - } - rayTracer.setParticleType(particle); - rayTracer.apply(); - - // fill up rates vector with rates from this particle type - auto &localData = rayTracer.getLocalData(); - int numRates = particle->getLocalDataLabels().size(); - for (int i = 0; i < numRates; ++i) { - auto rate = std::move(localData.getVectorData(i)); - - // normalize fluxes - rayTracer.normalizeFlux(rate, rayTracingParams.normalizationType); - if (rayTracingParams.smoothingNeighbors > 0) - rayTracer.smoothFlux(rate, rayTracingParams.smoothingNeighbors); - rates->insertNextScalarData(std::move(rate), - localData.getVectorDataLabel(i)); - } + // update coverages in the model + model->getSurfaceModel()->updateCoverages(fluxes, materialIds); - if (dataLogSize > 0) { - particleDataLogs[particleIdx].merge(rayTracer.getDataLog()); - } - ++particleIdx; - } + // metric to check for convergence + coverages = model->getSurfaceModel()->getCoverages(); + auto metric = + calculateCoverageDeltaMetric(coverages, prevStepCoverages); - // move coverages back in the model - moveRayDataToPointData(model->getSurfaceModel()->getCoverages(), - rayTraceCoverages); - model->getSurfaceModel()->updateCoverages(rates, materialIds); + if (logLevel >= 3) { + mergeScalarData(diskMesh->getCellData(), coverages); + mergeScalarData(diskMesh->getCellData(), fluxes); + printDiskMesh(diskMesh, name + "_covIinit_" + + std::to_string(iteration) + ".vtp"); + Logger::getInstance() + .addInfo("Iteration: " + std::to_string(iteration + 1)) + .print(); - // check for convergence - if (logLevel >= 5) { - coverages = model->getSurfaceModel()->getCoverages(); - auto metric = - calculateCoverageDeltaMetric(coverages, prevStepCoverages); - for (auto val : metric) { - covMetricFile << val << ";"; + 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"; } - covMetricFile << "\n"; } - if (logLevel >= 3) { - auto coverages = model->getSurfaceModel()->getCoverages(); - for (size_t idx = 0; idx < coverages->getScalarDataSize(); idx++) { - auto label = coverages->getScalarDataLabel(idx); - diskMesh->getCellData().insertNextScalarData( - *coverages->getScalarData(idx), label); - } - for (size_t idx = 0; idx < rates->getScalarDataSize(); idx++) { - auto label = rates->getScalarDataLabel(idx); - diskMesh->getCellData().insertNextScalarData( - *rates->getScalarData(idx), label); - } - printDiskMesh(diskMesh, name + "_covIinit_" + - std::to_string(iterations) + ".vtp"); + if (checkCoveragesConvergence(metric)) { Logger::getInstance() - .addInfo("Iteration: " + std::to_string(iterations)) + .addInfo("Coverages converged after " + + std::to_string(iteration + 1) + " iterations.") .print(); + break; } - } + } // end coverage initialization loop coveragesInitialized_ = true; timer.finish(); @@ -490,7 +471,7 @@ template class Process { advectionKernel.prepareLS(); model->initialize(domain, remainingTime); - auto rates = SmartPointer>::New(); + auto fluxes = SmartPointer>::New(); meshConverter.apply(); auto materialIds = *diskMesh->getCellData().getScalarData("MaterialIds"); auto points = diskMesh->getNodes(); @@ -507,88 +488,60 @@ template class Process { } rayTracer.setMaterialIds(materialIds); - // move coverages to ray tracer - viennaray::TracingData rayTraceCoverages; - if (useCoverages) { - rayTraceCoverages = - movePointDataToRayData(model->getSurfaceModel()->getCoverages()); - if (useProcessParams) { - // store scalars in addition to coverages - auto processParams = - model->getSurfaceModel()->getProcessParameters(); - NumericType numParams = processParams->getScalarData().size(); - rayTraceCoverages.setNumberOfScalarData(numParams); - for (size_t i = 0; i < numParams; ++i) { - rayTraceCoverages.setScalarData( - i, processParams->getScalarData(i), - processParams->getScalarDataLabel(i)); - } - } - rayTracer.setGlobalData(rayTraceCoverages); - } - - std::size_t particleIdx = 0; - for (auto &particle : model->getParticleTypes()) { - int dataLogSize = model->getParticleLogSize(particleIdx); - if (dataLogSize > 0) { - rayTracer.getDataLog().data.resize(1); - rayTracer.getDataLog().data[0].resize(dataLogSize, 0.); - } - rayTracer.setParticleType(particle); - rayTracer.apply(); - - // fill up rates vector with rates from this particle type - auto numRates = particle->getLocalDataLabels().size(); - auto &localData = rayTracer.getLocalData(); - for (int i = 0; i < numRates; ++i) { - auto rate = std::move(localData.getVectorData(i)); - - // normalize rates - rayTracer.normalizeFlux(rate, rayTracingParams.normalizationType); - if (rayTracingParams.smoothingNeighbors > 0) - rayTracer.smoothFlux(rate, rayTracingParams.smoothingNeighbors); - rates->insertNextScalarData(std::move(rate), - localData.getVectorDataLabel(i)); - } - - if (dataLogSize > 0) { - particleDataLogs[particleIdx].merge(rayTracer.getDataLog()); - } - ++particleIdx; - } + fluxes = calculateFluxes(rayTracer, useCoverages, useProcessParams); - // move coverages back to model - if (useCoverages) - moveRayDataToPointData(model->getSurfaceModel()->getCoverages(), - rayTraceCoverages); rtTimer.finish(); Logger::getInstance() .addTiming("Top-down flux calculation", rtTimer) .print(); } - // save coverages for comparison - SmartPointer> prevStepCoverages; - if (useCoverages && logLevel >= 5) { + // update coverages and calculate coverage delta metric + if (useCoverages) { coverages = model->getSurfaceModel()->getCoverages(); - prevStepCoverages = + auto prevStepCoverages = SmartPointer>::New(*coverages); - } - // get velocities from rates - auto velocities = model->getSurfaceModel()->calculateVelocities( - rates, points, materialIds); - - // calculate coverage metric - if (useCoverages && logLevel >= 5) { - auto metric = - calculateCoverageDeltaMetric(coverages, prevStepCoverages); - for (auto val : metric) { - covMetricFile << val << ";"; + // update coverages in the model + model->getSurfaceModel()->updateCoverages(fluxes, materialIds); + + if (coverageDeltaTreshold > 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(); + } } - covMetricFile << "\n"; } + // calculate velocities from fluxes + auto velocities = model->getSurfaceModel()->calculateVelocities( + fluxes, points, materialIds); + // prepare velocity field model->getVelocityField()->prepare(domain, velocities, processDuration - remainingTime); @@ -602,17 +555,9 @@ template class Process { "velocities"); if (useCoverages) { auto coverages = model->getSurfaceModel()->getCoverages(); - for (size_t idx = 0; idx < coverages->getScalarDataSize(); idx++) { - auto label = coverages->getScalarDataLabel(idx); - diskMesh->getCellData().insertNextScalarData( - *coverages->getScalarData(idx), label); - } - } - for (size_t idx = 0; idx < rates->getScalarDataSize(); idx++) { - auto label = rates->getScalarDataLabel(idx); - diskMesh->getCellData().insertNextScalarData( - *rates->getScalarData(idx), label); + mergeScalarData(diskMesh->getCellData(), coverages); } + mergeScalarData(diskMesh->getCellData(), fluxes); printDiskMesh(diskMesh, name + "_" + std::to_string(counter) + ".vtp"); if (domain->getCellSet()) { domain->getCellSet()->writeVTU(name + "_cellSet_" + @@ -753,13 +698,23 @@ template class Process { } private: - void printDiskMesh(SmartPointer> mesh, - std::string name) const { + static void printDiskMesh(SmartPointer> mesh, + std::string name) { viennals::VTKWriter(mesh, std::move(name)).apply(); } - viennaray::TracingData movePointDataToRayData( - SmartPointer> pointData) const { + 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)); + } + } + + static viennaray::TracingData movePointDataToRayData( + SmartPointer> pointData) { viennaray::TracingData rayData; const auto numData = pointData->getScalarDataSize(); rayData.setNumberOfVectorData(numData); @@ -772,9 +727,9 @@ template class Process { return std::move(rayData); } - void moveRayDataToPointData( + static void moveRayDataToPointData( SmartPointer> pointData, - viennaray::TracingData &rayData) const { + viennaray::TracingData &rayData) { pointData->clear(); const auto numData = rayData.getVectorData().size(); for (size_t i = 0; i < numData; ++i) @@ -818,9 +773,9 @@ template class Process { } } - std::vector calculateCoverageDeltaMetric( + static std::vector calculateCoverageDeltaMetric( SmartPointer> updated, - SmartPointer> previous) const { + SmartPointer> previous) { assert(updated->getScalarDataSize() == previous->getScalarDataSize()); std::vector delta(updated->getScalarDataSize(), 0.); @@ -841,6 +796,15 @@ template class Process { return delta; } + bool + checkCoveragesConvergence(const std::vector &deltaMetric) const { + for (auto val : deltaMetric) { + if (val > coverageDeltaTreshold) + return false; + } + return true; + } + void initializeRayTracer(viennaray::Trace &tracer) const { // Map the domain boundary to the ray tracing boundaries viennaray::BoundaryCondition rayBoundaryCondition[D]; @@ -873,11 +837,83 @@ template class Process { } } + auto calculateFluxes(viennaray::Trace &rayTracer, + const bool useCoverages, const bool useProcessParams) { + + viennaray::TracingData rayTracingData; + + // move coverages to the ray tracer + if (useCoverages) { + rayTracingData = + movePointDataToRayData(model->getSurfaceModel()->getCoverages()); + } + + if (useProcessParams) { + // store scalars in addition to coverages + auto processParams = model->getSurfaceModel()->getProcessParameters(); + NumericType numParams = processParams->getScalarData().size(); + rayTracingData.setNumberOfScalarData(numParams); + for (size_t i = 0; i < numParams; ++i) { + rayTracingData.setScalarData(i, processParams->getScalarData(i), + processParams->getScalarDataLabel(i)); + } + } + + if (useCoverages || useProcessParams) + rayTracer.setGlobalData(rayTracingData); + + auto fluxes = runRayTracer(rayTracer); + + // move coverages back in the model + if (useCoverages) + moveRayDataToPointData(model->getSurfaceModel()->getCoverages(), + rayTracingData); + + return fluxes; + } + + auto runRayTracer(viennaray::Trace &rayTracer) { + auto fluxes = SmartPointer>::New(); + unsigned particleIdx = 0; + for (auto &particle : model->getParticleTypes()) { + int dataLogSize = model->getParticleLogSize(particleIdx); + if (dataLogSize > 0) { + rayTracer.getDataLog().data.resize(1); + rayTracer.getDataLog().data[0].resize(dataLogSize, 0.); + } + rayTracer.setParticleType(particle); + rayTracer.apply(); + + // fill up fluxes vector with fluxes from this particle type + auto &localData = rayTracer.getLocalData(); + int numFluxes = particle->getLocalDataLabels().size(); + for (int i = 0; i < numFluxes; ++i) { + auto flux = std::move(localData.getVectorData(i)); + + // normalize and smooth + rayTracer.normalizeFlux(flux, rayTracingParams.normalizationType); + if (rayTracingParams.smoothingNeighbors > 0) + rayTracer.smoothFlux(flux, rayTracingParams.smoothingNeighbors); + + fluxes->insertNextScalarData(std::move(flux), + localData.getVectorDataLabel(i)); + } + + if (dataLogSize > 0) { + particleDataLogs[particleIdx].merge(rayTracer.getDataLog()); + } + ++particleIdx; + } + return fluxes; + } + +private: psDomainType domain; SmartPointer> model; std::vector> particleDataLogs; NumericType processDuration; - unsigned maxIterations = 20; + unsigned maxIterations = std::numeric_limits::max(); + NumericType coverageDeltaTreshold = 0.; bool coveragesInitialized_ = false; NumericType processTime = 0.; diff --git a/pyproject.toml b/pyproject.toml index 6a86907e..c3ce6591 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,7 +8,7 @@ name = "ViennaPS" readme = "README.md" license = {file = "LICENSE"} description = "Semiconductor fabrication process simulation library" -dependencies = ['ViennaLS'] +dependencies = ["ViennaLS>=4.2.1"] [project.urls] Homepage = "https://viennatools.github.io/ViennaPS/" diff --git a/python/pyWrap.cpp b/python/pyWrap.cpp index 14da064c..8de2223c 100644 --- a/python/pyWrap.cpp +++ b/python/pyWrap.cpp @@ -1300,6 +1300,10 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) { .def("setMaxCoverageInitIterations", &Process::setMaxCoverageInitIterations, "Set the number of iterations to initialize the coverages.") + .def("setCoverageDeltaThreshold", + &Process::setCoverageDeltaThreshold, + "Set the threshold for the coverage delta metric to reach " + "convergence.") .def("setIntegrationScheme", &Process::setIntegrationScheme, "Set the integration scheme for solving the level-set equation. " "Possible integration schemes are specified in " diff --git a/python/viennaps2d/viennaps2d.pyi b/python/viennaps2d/viennaps2d.pyi index 7791234d..1cfd0d4d 100644 --- a/python/viennaps2d/viennaps2d.pyi +++ b/python/viennaps2d/viennaps2d.pyi @@ -752,6 +752,7 @@ class Process: def getProcessDuration(self) -> float: ... def getRayTracingParameters(self) -> RayTracingParameters: ... def setAdvectionParameters(self, arg0: AdvectionParameters) -> None: ... + def setCoverageDeltaThreshold(self, arg0: float) -> None: ... def setDomain(self, *args, **kwargs): ... def setIntegrationScheme( self, arg0: viennals2d.viennals2d.IntegrationSchemeEnum diff --git a/python/viennaps3d/viennaps3d.pyi b/python/viennaps3d/viennaps3d.pyi index 79be024f..252cd35f 100644 --- a/python/viennaps3d/viennaps3d.pyi +++ b/python/viennaps3d/viennaps3d.pyi @@ -709,6 +709,7 @@ class Process: def getProcessDuration(self) -> float: ... def getRayTracingParameters(self) -> RayTracingParameters: ... def setAdvectionParameters(self, arg0: AdvectionParameters) -> None: ... + def setCoverageDeltaThreshold(self, arg0: float) -> None: ... def setDomain(self, *args, **kwargs): ... def setIntegrationScheme( self, arg0: viennals3d.viennals3d.IntegrationSchemeEnum