diff --git a/Examples/ExampleProcess/SurfaceModel.hpp b/Examples/ExampleProcess/SurfaceModel.hpp index abc88a7f..c70b1326 100644 --- a/Examples/ExampleProcess/SurfaceModel.hpp +++ b/Examples/ExampleProcess/SurfaceModel.hpp @@ -20,9 +20,10 @@ class SurfaceModel : public psSurfaceModel { processParams->insertNextScalar(0., "processParameter"); } - psSmartPointer> - calculateVelocities(psSmartPointer> Rates, - const std::vector &materialIds) override { + psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { // use coverages and rates here to calculate the velocity here return psSmartPointer>::New( *Rates->getScalarData("particleRate")); diff --git a/Examples/Redeposition/CMakeLists.txt b/Examples/Redeposition/CMakeLists.txt new file mode 100644 index 00000000..e3fecfbc --- /dev/null +++ b/Examples/Redeposition/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required(VERSION 3.4) + +project("Redeposition") + +add_executable(${PROJECT_NAME} ${PROJECT_NAME}.cpp) +target_include_directories(${PROJECT_NAME} PUBLIC ${VIENNAPS_INCLUDE_DIRS}) +target_link_libraries(${PROJECT_NAME} PRIVATE ${VIENNAPS_LIBRARIES}) +configure_file(config.txt ${CMAKE_CURRENT_BINARY_DIR}/config.txt COPYONLY) + +add_dependencies(buildExamples ${PROJECT_NAME}) diff --git a/Examples/Redeposition/Parameters.hpp b/Examples/Redeposition/Parameters.hpp new file mode 100644 index 00000000..d1536525 --- /dev/null +++ b/Examples/Redeposition/Parameters.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include +#include + +#include + +template struct Parameters { + // Domain + T gridDelta = 0.2; + T xExtent = 10.0; + + // Geometry + int numLayers = 26; + T layerHeight = 2; + T substrateHeight = 4; + T holeRadius = 2; + + // Process + T diffusionCoefficient = 10.; // diffusion cofficient + T sink = 0.001; // sink strength + // convection velocity in the scallops towards the center + T scallopStreamVelocity = 10.; + // convection velocity in the center towards the sink on the top + T holeStreamVelocity = 10.; + + Parameters() {} + + void fromMap(std::unordered_map &m) { + psUtils::AssignItems( // + m, // + psUtils::Item{"gridDelta", gridDelta}, // + psUtils::Item{"xExtent", xExtent}, // + psUtils::Item{"numLayers", numLayers}, // + psUtils::Item{"layerHeight", layerHeight}, // + psUtils::Item{"substrateHeight", substrateHeight}, // + psUtils::Item{"holeRadius", holeRadius}, // + psUtils::Item{"diffusionCoefficient", diffusionCoefficient}, // + psUtils::Item{"sink", sink}, // + psUtils::Item{"scallopStreamVelocity", scallopStreamVelocity}, // + psUtils::Item{"holeStreamVelocity", holeStreamVelocity} // + ); + } +}; diff --git a/Examples/Redeposition/Redeposition.cpp b/Examples/Redeposition/Redeposition.cpp new file mode 100644 index 00000000..91bf9b00 --- /dev/null +++ b/Examples/Redeposition/Redeposition.cpp @@ -0,0 +1,97 @@ +#include +#include + +#include +#include +#include +#include + +#include + +#include "Parameters.hpp" + +int main(int argc, char **argv) { + using NumericType = double; + + omp_set_num_threads(12); + // Attention: + // This model/example currently only works in 2D mode + constexpr int D = 2; + + Parameters params; + if (argc > 1) { + auto config = psUtils::readConfigFile(argv[1]); + if (config.empty()) { + std::cerr << "Empty config provided" << std::endl; + return -1; + } + params.fromMap(config); + } + + auto domain = psSmartPointer>::New(); + psMakeStack builder( + domain, params.gridDelta, params.xExtent, 0., params.numLayers, + params.layerHeight, params.substrateHeight, params.holeRadius, false); + builder.apply(); + // copy top layer for deposition + auto depoLayer = psSmartPointer>::New( + domain->getLevelSets()->back()); + domain->insertNextLevelSet(depoLayer); + domain->generateCellSet(2., true /* true means cell set above surface */); + auto &cellSet = domain->getCellSet(); + cellSet->writeVTU("initial.vtu"); + // we need neighborhood information for solving the + // convection-diffusion equation on the cell set + cellSet->buildNeighborhood(); + + // The redeposition model captures byproducts from the selective etching + // process in the cell set. The byproducts are then distributed by solving a + // convection-diffusion equation on the cell set. + auto redepoModel = psSmartPointer>::New( + domain, params.diffusionCoefficient, params.sink, + params.scallopStreamVelocity, params.holeStreamVelocity, + builder.getTopLayer(), builder.getHeight(), builder.getHoleRadius()); + auto velField = psSmartPointer>::New(); + auto etchSurfModel = + psSmartPointer>::New(); + auto depoSurfModel = + psSmartPointer>::New( + cellSet, 1., builder.getHeight(), builder.getTopLayer()); + + // run the selective etching process + { + auto etchModel = psSmartPointer>::New(); + etchModel->setSurfaceModel(etchSurfModel); + etchModel->setVelocityField(velField); + etchModel->setAdvectionCallback(redepoModel); + etchModel->setProcessName("SelectiveEtching"); + + psProcess processEtch; + processEtch.setDomain(domain); + processEtch.setProcessModel(etchModel); + processEtch.setProcessDuration(50.); + + processEtch.apply(); + } + + // run the redeposition based on the captured byproducts + { + auto depoModel = psSmartPointer>::New(); + depoModel->setSurfaceModel(depoSurfModel); + depoModel->setVelocityField(velField); + depoModel->setProcessName("Redeposition"); + + psProcess processRedepo; + processRedepo.setDomain(domain); + processRedepo.setProcessModel(depoModel); + processRedepo.setProcessDuration(0.2); + + processRedepo.apply(); + cellSet->updateMaterials(); + cellSet->writeVTU("RedepositionCellSet.vtu"); + } + + psWriteVisualizationMesh(domain, "FinalStack").apply(); + + return 0; +} \ No newline at end of file diff --git a/Examples/Redeposition/config.txt b/Examples/Redeposition/config.txt new file mode 100644 index 00000000..01b1bff2 --- /dev/null +++ b/Examples/Redeposition/config.txt @@ -0,0 +1,16 @@ +# Config file for a redeposition example process +# Domain +gridDelta=0.1 +xExtent=10.0 + +# Geometry +int numLayers=26 +layerHeight=2 +substrateHeight=4 +holeRadius=2 + +# Process +diffusionCoefficient=10. +sink=0.001 +scallopStreamVelocity=2. +holeStreamVelocity=2. \ No newline at end of file diff --git a/README.md b/README.md index ca0e6dad..5fafe340 100644 --- a/README.md +++ b/README.md @@ -97,6 +97,13 @@ By changing the dimension of the hole etching example (_D = 2_), we can easily s +### Redeposition During Selective Etching + +This example demonstrates capturing etching byproducts and the subsequent redeposition during a selective etching process in a Si3N4/SiO2 stack. The etching byproducts are captured in a cell set description of the etching plasma. To model the dynamics of these etching byproducts, a convection-diffusion equation is solved on the cell set using finite differences. The redeposition is then captured by adding up the byproducts in every step and using this information to generate a velocity field on the etched surface. +
+ +
+ ## Application It is also possible to build an application which can parse a custom configuration file and execute pre-defined processes. The application can be built using CMake (make sure all dependencies are installed/ have been built previously): diff --git a/data/images/redeposition.gif b/data/images/redeposition.gif new file mode 100644 index 00000000..e18fb724 Binary files /dev/null and b/data/images/redeposition.gif differ diff --git a/include/CellSet/csDenseCellSet.hpp b/include/CellSet/csDenseCellSet.hpp index f4eccd57..1faa7cd8 100644 --- a/include/CellSet/csDenseCellSet.hpp +++ b/include/CellSet/csDenseCellSet.hpp @@ -14,20 +14,23 @@ #include +#include + /** This class represents a cell-based voxel implementation of a volume. The depth of the cell set in z-direction can be specified. */ template class csDenseCellSet { private: - using gridType = lsSmartPointer>; + using gridType = psSmartPointer>; using levelSetsType = - lsSmartPointer>>>; + psSmartPointer>>>; levelSetsType levelSets = nullptr; gridType cellGrid = nullptr; - lsSmartPointer> surface = nullptr; - lsSmartPointer> BVH = nullptr; + psSmartPointer> surface = nullptr; + psSmartPointer> BVH = nullptr; + std::vector> neighborhood; T gridDelta; size_t numberOfCells; T depth = 0.; @@ -49,10 +52,10 @@ template class csDenseCellSet { levelSets = passedLevelSets; if (cellGrid == nullptr) - cellGrid = lsSmartPointer>::New(); + cellGrid = psSmartPointer>::New(); if (surface == nullptr) - surface = lsSmartPointer>::New(levelSets->back()); + surface = psSmartPointer>::New(levelSets->back()); else surface->deepCopy(levelSets->back()); @@ -83,14 +86,14 @@ template class csDenseCellSet { } lsToVoxelMesh voxelConverter(cellGrid); - auto plane = lsSmartPointer>::New(surface->getGrid()); + auto plane = psSmartPointer>::New(surface->getGrid()); if (depth > 0.) { T origin[D] = {0.}; T normal[D] = {0.}; origin[D - 1] = depthPlanePos; normal[D - 1] = 1.; lsMakeGeometry(plane, - lsSmartPointer>::New(origin, normal)) + psSmartPointer>::New(origin, normal)) .apply(); } if (!cellSetAboveSurface && depth > 0.) @@ -113,7 +116,7 @@ template class csDenseCellSet { fillingFractions = cellGrid->getCellData().getScalarData("fillingFraction"); calculateBounds(minBounds, maxBounds); - BVH = lsSmartPointer>::New(getBoundingBox(), BVHlayers); + BVH = psSmartPointer>::New(getBoundingBox(), BVHlayers); buildBVH(); } @@ -134,7 +137,7 @@ template class csDenseCellSet { gridType getCellGrid() { return cellGrid; } - lsSmartPointer> getBVH() const { return BVH; } + psSmartPointer> getBVH() const { return BVH; } T getDepth() const { return depth; } @@ -144,11 +147,11 @@ template class csDenseCellSet { return cellGrid->getNodes(); } - std::vector> getElements() { + std::vector> &getElements() { return cellGrid->template getElements<(1 << D)>(); } - lsSmartPointer> getSurface() { return surface; } + psSmartPointer> getSurface() { return surface; } levelSetsType getLevelSets() const { return levelSets; } @@ -250,14 +253,14 @@ template class csDenseCellSet { lsToVoxelMesh voxelConverter(cellGrid); auto plane = - lsSmartPointer>::New(levelSets->back()->getGrid()); + psSmartPointer>::New(levelSets->back()->getGrid()); if (depth > 0.) { T origin[D] = {0.}; T normal[D] = {0.}; origin[D - 1] = depthPlanePos; normal[D - 1] = 1.; lsMakeGeometry(plane, - lsSmartPointer>::New(origin, normal)) + psSmartPointer>::New(origin, normal)) .apply(); } if (!cellSetAboveSurface && depth > 0.) @@ -287,18 +290,18 @@ template class csDenseCellSet { // Updates the surface of the cell set. The new surface should be below the // old surface as this function can only remove cells from the cell set. void updateSurface() { - auto updateCellGrid = lsSmartPointer>::New(); + auto updateCellGrid = psSmartPointer>::New(); lsToVoxelMesh voxelConverter(updateCellGrid); if (depth != 0.) { - auto plane = lsSmartPointer>::New(surface->getGrid()); + auto plane = psSmartPointer>::New(surface->getGrid()); T origin[D] = {0.}; T normal[D] = {0.}; origin[D - 1] = depthPlanePos; normal[D - 1] = 1.; lsMakeGeometry(plane, - lsSmartPointer>::New(origin, normal)) + psSmartPointer>::New(origin, normal)) .apply(); voxelConverter.insertNextLevelSet(plane); } @@ -346,6 +349,42 @@ template class csDenseCellSet { } } + void buildNeighborhood() { + const auto &cells = cellGrid->template getElements<(1 << D)>(); + unsigned const numNodes = cellGrid->getNodes().size(); + unsigned const numCells = cells.size(); + + neighborhood.clear(); + neighborhood.resize(numCells); + + std::vector> nodeCellConnections(numNodes); + + // for each node, store which cells are connected with the node + for (unsigned cellIdx = 0; cellIdx < numCells; cellIdx++) { + for (unsigned cellNodeIdx = 0; cellNodeIdx < (1 << D); cellNodeIdx++) { + nodeCellConnections[cells[cellIdx][cellNodeIdx]].push_back(cellIdx); + } + } + + for (unsigned cellIdx = 0; cellIdx < numCells; cellIdx++) { + for (unsigned cellNodeIdx = 0; cellNodeIdx < (1 << D); cellNodeIdx++) { + auto &cellsAtNode = nodeCellConnections[cells[cellIdx][cellNodeIdx]]; + for (const auto &neighborCell : cellsAtNode) { + if (neighborCell != cellIdx) { + neighborhood[cellIdx].insert(neighborCell); + } + } + } + } + } + + std::set &getNeighbors(size_t cellIdx) { + assert(!neighborhood.empty() && + "Querying neighbors without creating neighborhood structure"); + assert(cellIdx < numberOfCells && "Cell idx out of bounds"); + return neighborhood[cellIdx]; + } + private: int findIndex(const csTriple &point) { const auto &elems = cellGrid->template getElements<(1 << D)>(); @@ -453,4 +492,4 @@ template class csDenseCellSet { } }; -#endif +#endif \ No newline at end of file diff --git a/include/Geometries/psMakeFin.hpp b/include/Geometries/psMakeFin.hpp index 096db440..2fbbf71d 100644 --- a/include/Geometries/psMakeFin.hpp +++ b/include/Geometries/psMakeFin.hpp @@ -5,6 +5,9 @@ #include #include +/** + * Creates a fin gemeotry in in z(3D)/y(2D) direction. + */ template class psMakeFin { using LSPtrType = psSmartPointer>; using PSPtrType = psSmartPointer>; diff --git a/include/Geometries/psMakeHole.hpp b/include/Geometries/psMakeHole.hpp index aa328b01..fa0b7846 100644 --- a/include/Geometries/psMakeHole.hpp +++ b/include/Geometries/psMakeHole.hpp @@ -8,9 +8,9 @@ #include /** - Creates a hole geometry in z direction. For 2D geometries a regular trench is - created. -*/ + * Creates a hole geometry in z direction. For 2D geometries a regular trench is + * created. + */ template class psMakeHole { using LSPtrType = psSmartPointer>; using PSPtrType = psSmartPointer>; diff --git a/include/Geometries/psMakePlane.hpp b/include/Geometries/psMakePlane.hpp index 82325419..f67a7f3a 100644 --- a/include/Geometries/psMakePlane.hpp +++ b/include/Geometries/psMakePlane.hpp @@ -8,8 +8,8 @@ #include /** - Creates a plane in z(3D)/y(2D) direction. -*/ + * Creates a plane in z(3D)/y(2D) direction. + */ template class psMakePlane { using LSPtrType = psSmartPointer>; using PSPtrType = psSmartPointer>; diff --git a/include/Geometries/psMakeStack.hpp b/include/Geometries/psMakeStack.hpp new file mode 100644 index 00000000..35a40f8b --- /dev/null +++ b/include/Geometries/psMakeStack.hpp @@ -0,0 +1,170 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include + +/** + * Creates a stack of alternating material layers with an etched + * hole(3D)/trench(2D) in the middle. + */ +template class psMakeStack { + using PSPtrType = psSmartPointer>; + using LSPtrType = psSmartPointer>; + + const T gridDelta = 0.1; + const T xExtent = 10; + const T yExtent = 10; + T bounds[2 * D]; + T normal[D]; + T origin[D] = {0.}; + + const int numLayers = 11; + const T layerHeight = 2; + const T substrateHeight = 4; + const T holeRadius = 2; + const bool periodicBoundary = false; + + typename lsDomain::BoundaryType boundaryConds[D]; + + PSPtrType geometry; + +public: + psMakeStack(PSPtrType passedDomain, const int passedNumLayers, + const T passedGridDelta) + : numLayers(passedNumLayers), gridDelta(passedGridDelta), + geometry(passedDomain) { + init(); + } + + psMakeStack(PSPtrType passedDomain, const T passedGridDelta, + const T passedXExtent, const T passedYExtent, + const int passedNumLayers, const T passedLayerHeight, + const T passedSubstrateHeight, const T passedHoleRadius, + const bool periodic = false) + : geometry(passedDomain), gridDelta(passedGridDelta), + xExtent(passedXExtent), yExtent(passedYExtent), + numLayers(passedNumLayers), layerHeight(passedLayerHeight), + substrateHeight(passedSubstrateHeight), holeRadius(passedHoleRadius), + periodicBoundary(periodic) { + init(); + } + + void apply() { + if constexpr (D == 2) { + create2DGeometry(); + } else { + create3DGeometry(); + } + } + + int getTopLayer() const { return numLayers + 1; } + + T getHeight() const { return substrateHeight + numLayers * layerHeight; } + + T getHoleRadius() const { return holeRadius; } + +private: + void create2DGeometry() { + geometry->clear(); + + // Silicon substrate + auto substrate = LSPtrType::New(bounds, boundaryConds, gridDelta); + origin[D - 1] = substrateHeight; + auto plane = lsSmartPointer>::New(origin, normal); + lsMakeGeometry(substrate, plane).apply(); + + // Si3N4/SiO2 layers + T current = substrateHeight + layerHeight; + for (int i = 0; i < numLayers; ++i) { + auto ls = LSPtrType::New(bounds, boundaryConds, gridDelta); + origin[D - 1] = substrateHeight + layerHeight * (i + 1); + auto plane = lsSmartPointer>::New(origin, normal); + lsMakeGeometry(ls, plane).apply(); + geometry->insertNextLevelSet(ls); + } + + // cut out middle + auto cutOut = LSPtrType::New(bounds, boundaryConds, gridDelta); + T minPoint[D] = {-holeRadius, substrateHeight * 0.8}; + T maxPoint[D] = {holeRadius, + substrateHeight + layerHeight * numLayers + gridDelta}; + lsMakeGeometry(cutOut, + lsSmartPointer>::New(minPoint, maxPoint)) + .apply(); + + for (auto layer : *geometry->getLevelSets()) { + lsBooleanOperation(layer, cutOut, + lsBooleanOperationEnum::RELATIVE_COMPLEMENT) + .apply(); + } + } + + void create3DGeometry() { + geometry->clear(); + + // cut out middle + auto cutOut = LSPtrType::New(bounds, boundaryConds, gridDelta); + origin[D - 1] = substrateHeight * 0.8; + lsMakeGeometry( + cutOut, lsSmartPointer>::New( + origin, normal, (numLayers + 1) * layerHeight, holeRadius)) + .apply(); + + // Silicon substrate + auto substrate = LSPtrType::New(bounds, boundaryConds, gridDelta); + origin[D - 1] = substrateHeight; + auto plane = lsSmartPointer>::New(origin, normal); + lsMakeGeometry(substrate, plane).apply(); + lsBooleanOperation(substrate, cutOut, + lsBooleanOperationEnum::RELATIVE_COMPLEMENT) + .apply(); + + // Si3N4/SiO2 layers + for (int i = 0; i < numLayers; ++i) { + auto ls = LSPtrType::New(bounds, boundaryConds, gridDelta); + origin[D - 1] = substrateHeight + layerHeight * (i + 1); + auto plane = lsSmartPointer>::New(origin, normal); + lsMakeGeometry(ls, plane).apply(); + lsBooleanOperation(ls, cutOut, + lsBooleanOperationEnum::RELATIVE_COMPLEMENT) + .apply(); + geometry->insertNextLevelSet(ls); + } + } + + void init() { + bounds[0] = -xExtent; + bounds[1] = xExtent; + normal[0] = 0.; + if (periodicBoundary) + boundaryConds[0] = lsDomain::BoundaryType::PERIODIC_BOUNDARY; + else + boundaryConds[0] = lsDomain::BoundaryType::REFLECTIVE_BOUNDARY; + + if constexpr (D == 2) { + normal[1] = 1.; + bounds[2] = 0; + bounds[3] = layerHeight * numLayers + gridDelta; + boundaryConds[1] = lsDomain::BoundaryType::INFINITE_BOUNDARY; + } else { + normal[1] = 0.; + normal[2] = 1.; + bounds[2] = -yExtent; + bounds[3] = yExtent; + bounds[4] = 0; + bounds[5] = layerHeight * numLayers + gridDelta; + if (periodicBoundary) + boundaryConds[0] = lsDomain::BoundaryType::PERIODIC_BOUNDARY; + else + boundaryConds[0] = lsDomain::BoundaryType::REFLECTIVE_BOUNDARY; + boundaryConds[2] = lsDomain::BoundaryType::INFINITE_BOUNDARY; + } + } +}; \ No newline at end of file diff --git a/include/Geometries/psMakeTrench.hpp b/include/Geometries/psMakeTrench.hpp index f1933133..621e3eb1 100644 --- a/include/Geometries/psMakeTrench.hpp +++ b/include/Geometries/psMakeTrench.hpp @@ -8,8 +8,8 @@ #include /** - Creates a trench geometry in z(3D)/y(2D) direction. -*/ + * Creates a trench geometry in z(3D)/y(2D) direction. + */ template class psMakeTrench { using LSPtrType = psSmartPointer>; using PSPtrType = psSmartPointer>; diff --git a/include/Models/DirectionalEtching.hpp b/include/Models/DirectionalEtching.hpp index ea6fc97e..ee0e41ca 100644 --- a/include/Models/DirectionalEtching.hpp +++ b/include/Models/DirectionalEtching.hpp @@ -47,9 +47,10 @@ class DirectionalEtchVelocityField : public psVelocityField { template class DirectionalEtchingSurfaceModel : public psSurfaceModel { public: - psSmartPointer> - calculateVelocities(psSmartPointer> Rates, - const std::vector &materialIds) override { + psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { return nullptr; } }; diff --git a/include/Models/SF6O2Etching.hpp b/include/Models/SF6O2Etching.hpp index 7db0c3f2..2dc30063 100644 --- a/include/Models/SF6O2Etching.hpp +++ b/include/Models/SF6O2Etching.hpp @@ -41,9 +41,10 @@ class SF6O2SurfaceModel : public psSurfaceModel { Coverages->insertNextScalarData(cov, "oCoverage"); } - psSmartPointer> - calculateVelocities(psSmartPointer> Rates, - const std::vector &materialIds) override { + psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { updateCoverages(Rates); const auto numPoints = Rates->getScalarData(0)->size(); std::vector etchRate(numPoints, 0.); diff --git a/include/Models/SimpleDeposition.hpp b/include/Models/SimpleDeposition.hpp index 65347e74..48350fe8 100644 --- a/include/Models/SimpleDeposition.hpp +++ b/include/Models/SimpleDeposition.hpp @@ -12,9 +12,10 @@ template class SimpleDepositionSurfaceModel : public psSurfaceModel { public: - psSmartPointer> - calculateVelocities(psSmartPointer> Rates, - const std::vector &materialIds) override { + psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { const auto depoRate = Rates->getScalarData("depoRate"); return psSmartPointer>::New(*depoRate); diff --git a/include/Models/StackRedeposition.hpp b/include/Models/StackRedeposition.hpp new file mode 100644 index 00000000..0412e6ef --- /dev/null +++ b/include/Models/StackRedeposition.hpp @@ -0,0 +1,271 @@ +#pragma once + +#include + +#include + +#include +#include + +// The selective etching model works in accordance with the geometry generated +// by psMakeStack +template +class SelectiveEtchingSurfaceModel : public psSurfaceModel { +public: + psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { + + const auto numPoints = materialIds.size(); + std::vector etchRate(numPoints, 0.); + + for (size_t i = 0; i < numPoints; ++i) { + int matId = static_cast(materialIds[i]); + if (matId == 0) { + etchRate[i] = 0; + continue; + } + + if (matId % 2 == 0) { + etchRate[i] = -0.1; + } else { + etchRate[i] = 0; + } + } + + return psSmartPointer>::New(etchRate); + } +}; + +template +class RedepositionDynamics : public psAdvectionCalback { + using psAdvectionCalback::domain; + + psSmartPointer> cellSet = nullptr; + int plasmaMaterial = 0; + T diffusionCoefficient = 1.; + T sink = 1; + T scallopStreamVel = 1.; + T holeStreamVel = 1.; + T top; + T holeRadius; + SelectiveEtchingSurfaceModel surfModel; + lsSmartPointer> velocities; + std::vector> nodes; + static constexpr T eps = 1e-4; + bool redepoRun = false; + +public: + RedepositionDynamics(psSmartPointer> passedDomain, + const T passedDiffCoeff, const T passedSink, + const T passedScallopVel, const T passedHoleVel, + const int passedPlasmaMat, const T passedTop, + const T passedRadius) + : diffusionCoefficient(passedDiffCoeff), sink(passedSink), + scallopStreamVel(passedScallopVel), holeStreamVel(passedHoleVel), + plasmaMaterial(passedPlasmaMat), top(passedTop), + holeRadius(passedRadius) { + domain = passedDomain; + assert(domain->getUseCellSet()); + cellSet = domain->getCellSet(); + cellSet->addScalarData("byproductSum", 0.); + } + + bool applyPreAdvect(const T processTime) override final { + prepare(); + return true; + } + + bool applyPostAdvect(const T advectionTime) override final { + cellSet->updateMaterials(); + addByproducts(advectionTime); + diffuseByproducts(advectionTime); + return true; + } + +private: + void diffuseByproducts(const T timeStep) { + auto data = cellSet->getFillingFractions(); + auto materialIds = cellSet->getScalarData("Material"); + auto elems = cellSet->getElements(); + auto nodes = cellSet->getNodes(); + const auto gridDelta = cellSet->getGridDelta(); + // calculate time discretization + const T dt = + std::min(gridDelta * gridDelta / diffusionCoefficient * 0.245, 1.); + const int numSteps = static_cast(timeStep / dt); + const T C = dt * diffusionCoefficient / (gridDelta * gridDelta); + + for (int iter = 0; iter < numSteps; iter++) { + std::vector solution(data->size(), 0.); + + // std::cout << (*std::max_element(materialIds->begin(), + // materialIds->end())) + // << std::endl; + +#pragma omp parallel for + for (size_t e = 0; e < data->size(); e++) { + if (materialIds->at(e) != plasmaMaterial) { + continue; + } + + int numNeighbors = 0; + auto coord = nodes[elems[e][0]]; + std::array, D> gridCoords = {}; + for (int i = 0; i < D; i++) { + coord[i] += gridDelta / 2.; + gridCoords[i][0] = -1; + gridCoords[i][1] = -1; + } + + for (const auto n : cellSet->getNeighbors(e)) { + if (materialIds->at(n) != plasmaMaterial) + continue; + + auto neighborCoord = nodes[elems[n][0]]; + for (int i = 0; i < D; i++) + neighborCoord[i] += gridDelta / 2.; + + if (csUtil::distance(coord, neighborCoord) < gridDelta + 1e-4) { + + for (int d = 0; d < D; d++) { + if (coord[d] - neighborCoord[d] > eps) { + gridCoords[d][0] = n; + } else if (coord[d] - neighborCoord[d] < -eps) { + gridCoords[d][1] = n; + } + } + + solution[e] += data->at(n); + numNeighbors++; + } + } + + // diffusion + solution[e] = + data->at(e) + + C * (solution[e] - static_cast(numNeighbors) * data->at(e)); + + // sink at the top + if (coord[1] > top - gridDelta) { + solution[e] -= sink; + solution[e] = std::max(solution[e], 0.); + continue; + } + + // convection + if (std::abs(coord[0]) < holeRadius) { + // in hole + assert(gridCoords[1][1] != -1 && "holeStream up neighbor wrong"); + solution[e] -= dt / gridDelta * holeStreamVel * + (data->at(gridCoords[1][1]) - data->at(e)); + } else { + if (coord[0] < 0) { + // left side scallop - use forward difference + assert(gridCoords[0][1] != -1 && + "scallopStream right neighbor wrong"); + solution[e] -= dt / gridDelta * scallopStreamVel * + (data->at(gridCoords[0][1]) - data->at(e)); + } else { + // right side scallop - use backward difference + assert(gridCoords[0][0] != -1 && + "scallopStream left neighbor wrong"); + solution[e] += dt / gridDelta * scallopStreamVel * + (data->at(e) - data->at(gridCoords[0][0])); + } + } + } + *data = std::move(solution); + } + auto sum = cellSet->getScalarData("byproductSum"); + +#pragma omp parallel for + for (size_t e = 0; e < data->size(); e++) { + if (materialIds->at(e) != plasmaMaterial) { + continue; + } + + sum->at(e) += data->at(e); + } + } + + void prepare() { + lsToDiskMesh meshConv; + for (auto ls : *domain->getLevelSets()) + meshConv.insertNextLevelSet(ls); + + auto mesh = lsSmartPointer>::New(); + meshConv.setMesh(mesh); + meshConv.apply(); + + nodes = mesh->getNodes(); + + velocities = surfModel.calculateVelocities( + nullptr, nodes, *mesh->getCellData().getScalarData("MaterialIds")); + } + + void addByproducts(const T timeStep) { + for (size_t j = 0; j < nodes.size(); j++) { + cellSet->addFillingFraction(nodes[j], -velocities->at(j) * timeStep); + } + } +}; + +template +class RedepositionSurfaceModel : public psSurfaceModel { +public: + RedepositionSurfaceModel( + psSmartPointer> passedCellSet, + const NumericType passedFactor, const NumericType passedTop, + const int passedPlasmaMat) + : cellSet(passedCellSet), top(passedTop), plasmaMaterial(passedPlasmaMat), + depoMaterial(passedPlasmaMat - 1), depositionFactor(passedFactor) {} + + psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { + + const auto numPoints = materialIds.size(); + std::vector depoRate(numPoints, 0.); + auto ff = cellSet->getScalarData("byproductSum"); + auto cellMatIds = cellSet->getScalarData("Material"); + + for (size_t i = 0; i < numPoints; ++i) { + int matId = static_cast(materialIds[i]); + if (matId == 0 || matId == plasmaMaterial) + continue; + + if (matId % 2 != 0 || matId == depoMaterial) { + auto node = coordinates[i]; + + if (node[D - 1] >= top) + continue; + + auto cellIdx = cellSet->getIndex(node); + int n = 0; + for (const auto ni : cellSet->getNeighbors(cellIdx)) { + const int neighborMatId = cellMatIds->at(ni); + + if (neighborMatId == plasmaMaterial) { + depoRate[i] += ff->at(ni); + n++; + } + } + if (n > 1) + depoRate[i] /= static_cast(n); + depoRate[i] *= depositionFactor; + } + } + + return psSmartPointer>::New(depoRate); + } + +private: + psSmartPointer> cellSet; + const NumericType top = 0.; + const int plasmaMaterial = 0; + const int depoMaterial = 0; + const NumericType depositionFactor = 1.; +}; \ No newline at end of file diff --git a/include/Models/WetEtching.hpp b/include/Models/WetEtching.hpp index e6e2eba3..9c32e5c4 100644 --- a/include/Models/WetEtching.hpp +++ b/include/Models/WetEtching.hpp @@ -91,9 +91,10 @@ class WetEtchingVelocityField : public psVelocityField { template class WetEtchingSurfaceModel : public psSurfaceModel { public: - psSmartPointer> - calculateVelocities(psSmartPointer> Rates, - const std::vector &materialIds) override { + psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { return nullptr; } }; diff --git a/include/psProcess.hpp b/include/psProcess.hpp index 8b11d46d..ec9f2f79 100644 --- a/include/psProcess.hpp +++ b/include/psProcess.hpp @@ -263,9 +263,9 @@ template class psProcess { auto Rates = psSmartPointer>::New(); meshConverter.apply(); auto materialIds = *diskMesh->getCellData().getScalarData("MaterialIds"); + auto points = diskMesh->getNodes(); if (useRayTracing) { - auto points = diskMesh->getNodes(); auto normals = *diskMesh->getCellData().getVectorData("Normals"); rayTrace.setGeometry(points, normals, gridDelta); rayTrace.setMaterialIds(materialIds); @@ -313,8 +313,8 @@ template class psProcess { rayTraceCoverages); } - auto velocitites = - model->getSurfaceModel()->calculateVelocities(Rates, materialIds); + auto velocitites = model->getSurfaceModel()->calculateVelocities( + Rates, points, materialIds); model->getVelocityField()->setVelocities(velocitites); #ifdef VIENNAPS_VERBOSE @@ -334,9 +334,13 @@ template class psProcess { diskMesh->getCellData().insertNextScalarData(*Rates->getScalarData(idx), label); } - if (printIntermediate) - printDiskMesh(diskMesh, - name + "_" + std::to_string(counter++) + ".vtp"); + if (printIntermediate) { + printDiskMesh(diskMesh, name + "_" + std::to_string(counter) + ".vtp"); + if (domain->getUseCellSet()) { + domain->getCellSet()->writeVTU(name + "_cellSet_" + + std::to_string(counter++) + ".vtu"); + } + } #endif // apply advection callback if (useAdvectionCallback) { diff --git a/include/psSurfaceModel.hpp b/include/psSurfaceModel.hpp index d5f252a9..cf1ebe73 100644 --- a/include/psSurfaceModel.hpp +++ b/include/psSurfaceModel.hpp @@ -26,9 +26,10 @@ template class psSurfaceModel { return processParams; } - virtual psSmartPointer> - calculateVelocities(psSmartPointer> Rates, - const std::vector &materialIDs) { + virtual psSmartPointer> calculateVelocities( + psSmartPointer> Rates, + const std::vector> &coordinates, + const std::vector &materialIDs) { return psSmartPointer>::New(); } diff --git a/include/psVelocityField.hpp b/include/psVelocityField.hpp index c76bf033..3fae96a9 100644 --- a/include/psVelocityField.hpp +++ b/include/psVelocityField.hpp @@ -34,4 +34,23 @@ template class psVelocityField { virtual bool useTranslationField() const { return true; } }; +template +class psDefaultVelocityField : public psVelocityField { +public: + virtual NumericType + getScalarVelocity(const std::array &coordinate, int material, + const std::array &normalVector, + unsigned long pointId) override { + return velocities->at(pointId); + } + + void setVelocities( + psSmartPointer> passedVelocities) override { + velocities = passedVelocities; + } + +private: + psSmartPointer> velocities; +}; + #endif \ No newline at end of file diff --git a/include/psWriteVisualizationMesh.hpp b/include/psWriteVisualizationMesh.hpp index dafd70bc..4cbfd3af 100644 --- a/include/psWriteVisualizationMesh.hpp +++ b/include/psWriteVisualizationMesh.hpp @@ -19,10 +19,6 @@ template class psWriteVisualizationMesh { int i = 0; for (auto ls : *domain->getLevelSets()) { visMesh.insertNextLevelSet(ls); - auto mesh = psSmartPointer>::New(); - lsToSurfaceMesh(ls, mesh).apply(); - psVTKWriter(mesh, "visMesh_" + std::to_string(i++) + ".vtp") - .apply(); } visMesh.apply(); }