Skip to content

Commit 6061e52

Browse files
Adaptive time stepping when approaching interfaces in lsAdvect (#136)
* Limit time step at material interfaces and test with python SquareEtch example * Update logger * Reverted material interface behavior and added removeLastMaterial functionality * format * Fix material interface stability with adaptive sub-stepping of thin layers (soft landing) * Added WENO integration scheme * Add option to enable adaptive time stepping during advection * Add python bindings * Bump version * Add option to set the adaptive time stepping threshold * Fix MultiSurfaceMesh ctors --------- Co-authored-by: filipovic <filipovic@iue.tuwien.ac.at>
1 parent 96ab9eb commit 6061e52

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

44 files changed

+813
-330
lines changed

CMakeLists.txt

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
22
project(
33
ViennaLS
44
LANGUAGES CXX
5-
VERSION 5.2.1)
5+
VERSION 5.3.0)
66

77
# --------------------------------------------------------------------------------------------------------
88
# Library options
@@ -111,9 +111,9 @@ include(cmake/vtk.cmake)
111111

112112
CPMAddPackage(
113113
NAME ViennaCore
114-
VERSION 1.7.0
114+
VERSION 1.8.0
115115
GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore"
116-
OPTIONS "VIENNACORE_FORMAT_EXCLUDE docs/ build/"
116+
OPTIONS "VIENNACORE_FORMAT_EXCLUDE build/"
117117
EXCLUDE_FROM_ALL ${VIENNALS_BUILD_PYTHON})
118118

119119
CPMAddPackage(

README.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -157,7 +157,7 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum
157157

158158
* Installation with CPM
159159
```cmake
160-
CPMAddPackage("gh:viennatools/viennals@5.2.1")
160+
CPMAddPackage("gh:viennatools/viennals@5.3.0")
161161
```
162162

163163
* With a local installation

examples/Deposition/Deposition.cpp

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -110,10 +110,11 @@ int main() {
110110
// the other could be taken as a mask layer for advection
111111
advectionKernel.insertNextLevelSet(substrate);
112112
advectionKernel.insertNextLevelSet(newLayer);
113-
114113
advectionKernel.setVelocityField(velocities);
115114
// advectionKernel.setAdvectionTime(4.);
116115
unsigned counter = 1;
116+
advectionKernel.setIntegrationScheme(
117+
ls::IntegrationSchemeEnum::WENO_5TH_ORDER);
117118
for (NumericType time = 0; time < 4.;
118119
time += advectionKernel.getAdvectedTime()) {
119120
advectionKernel.apply();

examples/SquareEtch/SquareEtch.cpp

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -80,7 +80,7 @@ class analyticalField : public ls::VelocityField<double> {
8080
int main() {
8181

8282
constexpr int D = 2;
83-
omp_set_num_threads(1);
83+
omp_set_num_threads(8);
8484

8585
// Change this to use the analytical velocity field
8686
const bool useAnalyticalVelocity = false;
@@ -173,7 +173,8 @@ int main() {
173173
// Analytical velocity fields and dissipation coefficients
174174
// can only be used with this integration scheme
175175
advectionKernel.setIntegrationScheme(
176-
ls::IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
176+
// ls::IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
177+
ls::IntegrationSchemeEnum::WENO_5TH_ORDER);
177178
} else {
178179
// for numerical velocities, just use the default
179180
// integration scheme, which is not accurate for certain

examples/SquareEtch/SquareEtch.py

Lines changed: 186 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
import viennals as vls
2+
import argparse
3+
4+
vls.Logger.setLogLevel(vls.LogLevel.INFO)
5+
6+
# 1. Define the Velocity Field
7+
# We inherit from viennals.VelocityField to define custom logic in Python
8+
class EtchingField(vls.VelocityField):
9+
def getScalarVelocity(self, coordinate, material, normalVector, pointId):
10+
if material == 2:
11+
return -0.1
12+
if material == 1:
13+
return -1.
14+
return 0.0
15+
16+
def getVectorVelocity(self, coordinate, material, normalVector, pointId):
17+
return [0.0] * len(coordinate)
18+
19+
# 1. Define the Velocity Field
20+
# We inherit from viennals.VelocityField to define custom logic in Python
21+
class DepositionField(vls.VelocityField):
22+
def getScalarVelocity(self, coordinate, material, normalVector, pointId):
23+
return 0.75
24+
25+
def getVectorVelocity(self, coordinate, material, normalVector, pointId):
26+
return [0.0] * len(coordinate)
27+
28+
def main():
29+
# 1. Parse Arguments
30+
parser = argparse.ArgumentParser(description="Run Square Etch simulation in 2D or 3D.")
31+
parser.add_argument(
32+
"-D", "--dim",
33+
type=int,
34+
default=2,
35+
choices=[2, 3],
36+
help="Dimension of the simulation (2 or 3). Default is 2."
37+
)
38+
args = parser.parse_args()
39+
40+
DIMENSION = args.dim
41+
vls.setDimension(DIMENSION)
42+
vls.setNumThreads(8)
43+
44+
extent = 30.0
45+
gridDelta = 0.47
46+
47+
# Define bounds and boundary conditions
48+
trenchBottom = -2.0
49+
bounds = [-extent, extent, -extent, extent]
50+
origin = [0.0, 0.0]
51+
if DIMENSION == 3:
52+
bounds.append(-extent)
53+
bounds.append(extent)
54+
origin.append(0.0)
55+
boundaryCons = []
56+
for i in range(DIMENSION - 1):
57+
boundaryCons.append(vls.BoundaryConditionEnum.REFLECTIVE_BOUNDARY)
58+
boundaryCons.append(vls.BoundaryConditionEnum.INFINITE_BOUNDARY)
59+
60+
planeNormal = list(origin)
61+
downNormal = list(origin)
62+
planeNormal[DIMENSION - 1] = 1.0
63+
downNormal[DIMENSION - 1] = -1.0
64+
65+
# Create the Substrate Domain
66+
substrate = vls.Domain(bounds, boundaryCons, gridDelta)
67+
68+
# Create initial flat geometry (Plane at y/z=0)
69+
plane = vls.Plane(origin, planeNormal)
70+
vls.MakeGeometry(substrate, plane).apply()
71+
72+
# --------------------------------------
73+
# 3. Trench Geometry
74+
# --------------------------------------
75+
trench = vls.Domain(bounds, boundaryCons, gridDelta)
76+
77+
# Define Box Corners based on dimension
78+
minCorner = list(origin)
79+
maxCorner = list(origin)
80+
originMask = list (origin)
81+
for i in range(DIMENSION - 1):
82+
minCorner[i] = -extent / 1.5
83+
maxCorner[i] = extent / 1.5
84+
minCorner[DIMENSION - 1] = trenchBottom
85+
maxCorner[DIMENSION - 1] = 1.0
86+
originMask[DIMENSION - 1] = trenchBottom + 1e-9
87+
88+
box = vls.Box(minCorner, maxCorner)
89+
vls.MakeGeometry(trench, box).apply()
90+
91+
# Subtract trench from substrate (Relative Complement)
92+
vls.BooleanOperation(substrate, trench, vls.BooleanOperationEnum.RELATIVE_COMPLEMENT).apply()
93+
94+
# 4. Create the Mask Layer
95+
# The mask prevents etching at the bottom of the trench in specific configurations,
96+
# or acts as a hard stop/masking material.
97+
mask = vls.Domain(bounds, boundaryCons, gridDelta)
98+
99+
vls.MakeGeometry(mask, vls.Plane(originMask, downNormal)).apply()
100+
101+
# Intersect with substrate geometry
102+
vls.BooleanOperation(mask, substrate, vls.BooleanOperationEnum.INTERSECT).apply()
103+
104+
# 5. Export Initial State
105+
print(f"Running in {DIMENSION}D.")
106+
print("Extracting initial meshes...")
107+
mesh = vls.Mesh()
108+
109+
# Save substrate
110+
vls.ToSurfaceMesh(substrate, mesh).apply()
111+
vls.VTKWriter(mesh, f"substrate-{DIMENSION}D.vtp").apply()
112+
113+
# Save mask
114+
vls.ToSurfaceMesh(mask, mesh).apply()
115+
vls.VTKWriter(mesh, f"mask-{DIMENSION}D.vtp").apply()
116+
117+
print("Creating new layer...")
118+
polymer = vls.Domain(substrate)
119+
120+
121+
# 6. Setup Advection
122+
deposition = DepositionField()
123+
etching = EtchingField()
124+
125+
print("Advecting")
126+
advectionKernel = vls.Advect()
127+
128+
# the level set to be advected has to be inserted last
129+
# the other could be taken as a mask layer for advection
130+
advectionKernel.insertNextLevelSet(mask)
131+
advectionKernel.insertNextLevelSet(substrate)
132+
advectionKernel.insertNextLevelSet(polymer)
133+
134+
advectionKernel.setSaveAdvectionVelocities(True)
135+
# advectionKernel.setVelocityField(etching)
136+
137+
advectionKernel.setVelocityField(deposition)
138+
advectionKernel.apply()
139+
140+
vls.ToSurfaceMesh(polymer, mesh).apply()
141+
vls.VTKWriter(mesh, f"newLayer-{DIMENSION}D.vtp").apply()
142+
143+
advectionKernel.setSaveAdvectionVelocities(True)
144+
advectionKernel.setVelocityField(etching)
145+
146+
# Use default integration scheme (Lax Friedrichs 1st order) as in the C++ else branch
147+
# advectionKernel.setIntegrationScheme(viennals.IntegrationSchemeEnum.LAX_FRIEDRICHS_1ST_ORDER)
148+
149+
# 7. Time Loop
150+
finalTime = 10.0
151+
counter = 1
152+
currentTime = 0.0
153+
154+
# The advect kernel calculates stable time steps automatically.
155+
# We call apply() repeatedly until we reach finalTime.
156+
157+
# Note: In the C++ example, the loop structure is:
158+
# for (double time = 0.; time < finalTime; time += advectionKernel.getAdvectedTime())
159+
# This implies one step is taken, then time is updated.
160+
161+
# Save initial state
162+
vls.ToSurfaceMesh(polymer, mesh).apply()
163+
vls.VTKWriter(mesh, f"numerical-{DIMENSION}D-0.vtp").apply()
164+
165+
while currentTime < finalTime:
166+
advectionKernel.apply()
167+
168+
stepTime = advectionKernel.getAdvectedTime()
169+
currentTime += stepTime
170+
171+
print(f"Advection step: {counter}, time: {stepTime:.4f} (Total: {currentTime:.4f})")
172+
173+
# Export result
174+
vls.ToSurfaceMesh(polymer, mesh).apply()
175+
vls.VTKWriter(mesh, f"numerical-{DIMENSION}D-{counter}.vtp").apply()
176+
177+
counter += 1
178+
179+
print(f"\nNumber of Advection steps taken: {counter}")
180+
181+
# Final export
182+
vls.ToSurfaceMesh(polymer, mesh).apply()
183+
vls.VTKWriter(mesh, f"final-{DIMENSION}D.vtp").apply()
184+
185+
if __name__ == "__main__":
186+
main()

0 commit comments

Comments
 (0)