Skip to content
Snippets Groups Projects
Commit 9de3da34 authored by Martin Bauer's avatar Martin Bauer
Browse files

pytest and coverage improvements 

- notebooks are considered in coverage analysis now
   -> custom notebook runner plugin adapted
- added images for documentation
parent ae4bbde0
Branches
Tags
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 647 additions and 280 deletions
doc/modules/lbmpy/materials/feature_optimization_overview.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/feature_optimization_overview.png

171 KiB

doc/modules/lbmpy/materials/lbm_01.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/lbm_01.png

29.4 KiB

doc/modules/lbmpy/materials/lbm_02.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/lbm_02.png

26.9 KiB

doc/modules/lbmpy/materials/lbm_03.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/lbm_03.png

29 KiB

doc/modules/lbmpy/materials/lbm_04.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/lbm_04.png

29.3 KiB

doc/modules/lbmpy/materials/lbm_05.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/lbm_05.png

29.1 KiB

doc/modules/lbmpy/materials/lbm_algorithm.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/lbm_algorithm.png

135 KiB

doc/modules/lbmpy/materials/simsetup.svg 0 → 100644
+ 425
− 0
View file @ 9de3da34
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:cc="http://creativecommons.org/ns#"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns="http://www.w3.org/2000/svg"
xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
width="818.26257"
height="382.85757"
id="svg2"
version="1.1"
inkscape:version="0.48.5 r10040"
sodipodi:docname="nondimensionalization.svg">
<defs
id="defs4">
<marker
style="overflow:visible"
id="DistanceEnd"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="DistanceEnd">
<g
id="g2301">
<path
style="fill:none;stroke:#ffffff;stroke-width:1.14999998;stroke-linecap:square"
d="M 0,0 -2,0"
id="path2316"
inkscape:connector-curvature="0" />
<path
style="fill:#000000;fill-rule:evenodd;stroke:none"
d="M 0,0 -13,4 -9,0 -13,-4 0,0 z"
id="path2312"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:square"
d="M 0,-4 0,40"
id="path2314"
inkscape:connector-curvature="0" />
</g>
</marker>
<marker
style="overflow:visible"
id="DistanceStart"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="DistanceStart">
<g
id="g2300">
<path
style="fill:none;stroke:#ffffff;stroke-width:1.14999998;stroke-linecap:square"
d="M 0,0 2,0"
id="path2306"
inkscape:connector-curvature="0" />
<path
style="fill:#000000;fill-rule:evenodd;stroke:none"
d="M 0,0 13,4 9,0 13,-4 0,0 z"
id="path2302"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:square"
d="M 0,-4 0,40"
id="path2304"
inkscape:connector-curvature="0" />
</g>
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend"
style="overflow:visible">
<path
id="path3882"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)"
inkscape:connector-curvature="0" />
</marker>
<marker
inkscape:stockid="Arrow1Lstart"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lstart"
style="overflow:visible">
<path
id="path3879"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(0.8,0,0,0.8,10,0)"
inkscape:connector-curvature="0" />
</marker>
<marker
inkscape:stockid="Arrow1Lstart"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lstart-3"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path3879-7"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(0.8,0,0,0.8,10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-3"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path3882-3"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lstart"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lstart-6"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path3879-3"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(0.8,0,0,0.8,10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-9"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path3882-35"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lstart"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lstart-7"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path3879-0"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(0.8,0,0,0.8,10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-96"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path3882-8"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
style="overflow:visible"
id="DistanceStart-3"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="DistanceStart">
<g
id="g2300-4">
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#ffffff;stroke-width:1.14999998;stroke-linecap:square"
d="M 0,0 2,0"
id="path2306-7" />
<path
inkscape:connector-curvature="0"
style="fill:#000000;fill-rule:evenodd;stroke:none"
d="M 0,0 13,4 9,0 13,-4 0,0 z"
id="path2302-0" />
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:square"
d="M 0,-4 0,40"
id="path2304-5" />
</g>
</marker>
<marker
style="overflow:visible"
id="DistanceEnd-3"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="DistanceEnd">
<g
id="g2301-2">
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#ffffff;stroke-width:1.14999998;stroke-linecap:square"
d="M 0,0 -2,0"
id="path2316-8" />
<path
inkscape:connector-curvature="0"
style="fill:#000000;fill-rule:evenodd;stroke:none"
d="M 0,0 -13,4 -9,0 -13,-4 0,0 z"
id="path2312-2" />
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:square"
d="M 0,-4 0,40"
id="path2314-7" />
</g>
</marker>
<marker
style="overflow:visible"
id="DistanceStart-0"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="DistanceStart">
<g
id="g2300-2">
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#ffffff;stroke-width:1.14999998;stroke-linecap:square"
d="M 0,0 2,0"
id="path2306-9" />
<path
inkscape:connector-curvature="0"
style="fill:#000000;fill-rule:evenodd;stroke:none"
d="M 0,0 13,4 9,0 13,-4 0,0 z"
id="path2302-9" />
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:square"
d="M 0,-4 0,40"
id="path2304-57" />
</g>
</marker>
<marker
style="overflow:visible"
id="DistanceEnd-5"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="DistanceEnd">
<g
id="g2301-22">
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#ffffff;stroke-width:1.14999998;stroke-linecap:square"
d="M 0,0 -2,0"
id="path2316-4" />
<path
inkscape:connector-curvature="0"
style="fill:#000000;fill-rule:evenodd;stroke:none"
d="M 0,0 -13,4 -9,0 -13,-4 0,0 z"
id="path2312-1" />
<path
inkscape:connector-curvature="0"
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:square"
d="M 0,-4 0,40"
id="path2314-2" />
</g>
</marker>
</defs>
<sodipodi:namedview
id="base"
pagecolor="#ffffff"
bordercolor="#666666"
borderopacity="1.0"
inkscape:pageopacity="0.0"
inkscape:pageshadow="2"
inkscape:zoom="0.98994949"
inkscape:cx="417.96939"
inkscape:cy="72.886099"
inkscape:document-units="px"
inkscape:current-layer="layer1"
showgrid="false"
inkscape:window-width="1600"
inkscape:window-height="1180"
inkscape:window-x="2557"
inkscape:window-y="-3"
inkscape:window-maximized="1"
fit-margin-top="20"
fit-margin-left="20"
fit-margin-bottom="20"
fit-margin-right="20">
<inkscape:grid
type="xygrid"
id="grid2985"
empspacing="5"
visible="true"
enabled="true"
snapvisiblegridlinesonly="true"
originx="92.030182px"
originy="-639.62961px" />
</sodipodi:namedview>
<metadata
id="metadata7">
<rdf:RDF>
<cc:Work
rdf:about="">
<dc:format>image/svg+xml</dc:format>
<dc:type
rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
<dc:title></dc:title>
</cc:Work>
</rdf:RDF>
</metadata>
<g
inkscape:label="Layer 1"
inkscape:groupmode="layer"
id="layer1"
transform="translate(92.030182,-29.875)">
<path
style="fill:none;stroke:#000000;stroke-width:5;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 0,52.362183 590,0"
id="path2989"
inkscape:connector-curvature="0"
sodipodi:nodetypes="cc" />
<path
style="fill:none;stroke:#000000;stroke-width:5;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 0,242.36218 590,0"
id="path2989-3"
inkscape:connector-curvature="0"
sodipodi:nodetypes="cc" />
<path
style="fill:none;stroke:#0000ff;stroke-width:5;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:0.59420286;stroke-dasharray:none"
d="M 2,52.362183 2,242.36218"
id="path3802"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#ff7f49;stroke-width:5;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:4;stroke-opacity:0.65700486;stroke-dasharray:none"
d="m 588,52.362186 0,189.999994"
id="path3802-0"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-start:url(#DistanceStart);marker-end:url(#DistanceEnd)"
d="m 0,298.36218 590,0"
id="path3873"
inkscape:connector-curvature="0" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="259.60919"
y="335.07266"
id="text4505"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan4507"
x="259.60919"
y="335.07266"
style="font-size:20px;font-style:italic;-inkscape-font-specification:Sans Italic">30 cm</tspan></text>
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-start:url(#DistanceStart);marker-end:url(#DistanceEnd)"
d="m 640,52.362183 0,189.999997"
id="path3873-1"
inkscape:connector-curvature="0"
sodipodi:nodetypes="cc" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="658"
y="154.36218"
id="text4505-6"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan4507-8"
x="658"
y="154.36218"
style="font-size:20px;font-style:italic;-inkscape-font-specification:Sans Italic">6 cm</tspan></text>
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-start:url(#Arrow1Lstart-7);marker-end:url(#Arrow1Lend-96)"
d="m -49.497475,328.08281 1.010153,45.45686 0,0 45.4568644,-1.01015"
id="path6047"
inkscape:connector-curvature="0" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="-10.101525"
y="392.73257"
id="text6727"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan6729"
x="-10.101525"
y="392.73257"
style="font-size:15px;font-weight:bold;-inkscape-font-specification:Sans Bold">x</tspan></text>
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="-72.213287"
y="339.9407"
id="text6727-1"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan6729-1"
x="-72.213287"
y="339.9407"
style="font-size:15px;font-weight:bold;-inkscape-font-specification:Sans Bold">y</tspan></text>
</g>
</svg>
doc/modules/lbmpy/materials/wing.png 0 → 100644
+ 0
− 0
View file @ 9de3da34
doc/modules/lbmpy/materials/wing.png

1.29 KiB

doc/modules/lbmpy/tutorial01_predefinedScenarios.ipynb
+ 18
− 18
View file @ 9de3da34
This diff is collapsed.
doc/modules/pyphasefield/Chapter4.ipynb
+ 19
− 50
View file @ 9de3da34
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## 4.4 Case Study-IV ## 4.4 Case Study-IV
- Phase-Field modeling of dendritic solidification - Phase-Field modeling of dendritic solidification
- Demonstrate the implementation of the anisotropy in the interfacial energy and the mobility that appears in phase-field models. - Demonstrate the implementation of the anisotropy in the interfacial energy and the mobility that appears in phase-field models.
- Two variables $\phi$ & $T$ - Two variables $\phi$ & $T$
- time: $t$, position $r$ - time: $t$, position $r$
- $F(\phi, m) = \int_V \frac{1}{2} \varepsilon^2 | \nabla \varphi |^2 + f(\varphi,m)dv $ Ginzburg-Landau type - $F(\phi, m) = \int_V \frac{1}{2} \varepsilon^2 | \nabla \varphi |^2 + f(\varphi,m)dv $ Ginzburg-Landau type
- Fist term: gradient energy - Fist term: gradient energy
- $\varepsilon$: anisotropic gradient energy coefficient (determines the thickness of interface layer) - $\varepsilon$: anisotropic gradient energy coefficient (determines the thickness of interface layer)
- Second term: local free energy (double well potential, with local minima at $\phi = 0, 1$ - Second term: local free energy (double well potential, with local minima at $\phi = 0, 1$
- $f(\phi, m) = \frac{1}{4}\varphi^4 - (\frac{1}{2}-\frac{1}{3}m)\varphi^3+(\frac{1}{4} - \frac{1}{2}m)\varphi^2$: functional form for free energy - $f(\phi, m) = \frac{1}{4}\varphi^4 - (\frac{1}{2}-\frac{1}{3}m)\varphi^3+(\frac{1}{4} - \frac{1}{2}m)\varphi^2$: functional form for free energy
- Interfacial anisotropy introduced by $\varepsilon$, which depends on direction of the outer normal vector - Interfacial anisotropy introduced by $\varepsilon$, which depends on direction of the outer normal vector
- $\varepsilon = \bar{\varepsilon} \sigma(\theta)$ - $\varepsilon = \bar{\varepsilon} \sigma(\theta)$
- $\bar{\varepsilon}$: mean value of $\varepsilon$ - $\bar{\varepsilon}$: mean value of $\varepsilon$
- $\sigma(\theta) = 1 + \delta \cos(j(\theta - \theta_0))$: anisotropy - $\sigma(\theta) = 1 + \delta \cos(j(\theta - \theta_0))$: anisotropy
- $\delta$: strength of anisotropy - $\delta$: strength of anisotropy
- $j$: mode number - $j$: mode number
- $\theta_0$: initial offset - $\theta_0$: initial offset
- $\theta = \tan^{-1} ( \frac{\partial \varphi/ \partial y}{\partial \varphi / \partial x} )$ - $\theta = \tan^{-1} ( \frac{\partial \varphi/ \partial y}{\partial \varphi / \partial x} )$
- $m(T) = \frac{\alpha}{\pi} \tan^{-1} (\gamma (T_{eq} - T))$ - $m(T) = \frac{\alpha}{\pi} \tan^{-1} (\gamma (T_{eq} - T))$
- $\alpha$: positive constant - $\alpha$: positive constant
- $T_{eq}$: equilibrium temperature - $T_{eq}$: equilibrium temperature
- $\tau \frac{\partial \varphi}{\partial t} = - \frac{\delta F}{\delta \Phi}$ - $\tau \frac{\partial \varphi}{\partial t} = - \frac{\delta F}{\delta \Phi}$
- $\Rightarrow \tau \frac{\partial \varphi}{\partial t} = \frac{\partial}{\partial y} (\varepsilon \frac{\partial \varepsilon}{\partial \theta} \frac{\partial \varphi}{\partial x}) - \frac{\partial}{\partial x} (\varepsilon \frac{\partial \varepsilon}{\partial \theta} \frac{\partial \varphi}{\partial y}) + \nabla \cdot (\varepsilon^2 \nabla \varphi) + \varphi(1-\varphi) ( \varphi - \frac{1}{2} + m)$ - $\Rightarrow \tau \frac{\partial \varphi}{\partial t} = \frac{\partial}{\partial y} (\varepsilon \frac{\partial \varepsilon}{\partial \theta} \frac{\partial \varphi}{\partial x}) - \frac{\partial}{\partial x} (\varepsilon \frac{\partial \varepsilon}{\partial \theta} \frac{\partial \varphi}{\partial y}) + \nabla \cdot (\varepsilon^2 \nabla \varphi) + \varphi(1-\varphi) ( \varphi - \frac{1}{2} + m)$
- Evolution of temperature field: - Evolution of temperature field:
- $\frac{\partial T}{\partial t} = \nabla^2 T+ \kappa \frac{\partial \varphi}{\partial t}$ - $\frac{\partial T}{\partial t} = \nabla^2 T+ \kappa \frac{\partial \varphi}{\partial t}$
- $T$ is nondimensionalized $\rightarrow$ cooling temperature $=0$, equilibirum $=1$ - $T$ is nondimensionalized $\rightarrow$ cooling temperature $=0$, equilibirum $=1$
- $\kappa$: dimensionless latent heat - $\kappa$: dimensionless latent heat
- (diffusion constant are identical in liquid and solid) - (diffusion constant are identical in liquid and solid)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#4.4.4 #4.4.4
import numpy as np import numpy as np
from timeit import default_timer as timer from timeit import default_timer as timer
#from pyevtk.hl import imageToVTK #from pyevtk.hl import imageToVTK
import math import math
def nucleus(Nx, Ny, seed): def nucleus(Nx, Ny, seed):
T = np.zeros((Nx+2, Ny+2)) T = np.zeros((Nx+2, Ny+2))
phi = np.zeros((Nx+2, Ny+2)) phi = np.zeros((Nx+2, Ny+2))
for i in range(0,Nx+2): for i in range(0,Nx+2):
for j in range(0,Ny+2): for j in range(0,Ny+2):
if (i-(Nx+2)/2)**2 + (j-(Ny+2)/2)**2 < seed: if (i-(Nx+2)/2)**2 + (j-(Ny+2)/2)**2 < seed:
phi[i,j] = 1.0 phi[i,j] = 1.0
return (phi, T) return (phi, T)
def laplacian(A, Nx, Ny, dx, dy): def laplacian(A, Nx, Ny, dx, dy):
out = np.zeros((Nx+2, Ny+2)) out = np.zeros((Nx+2, Ny+2))
out[1:Nx+1, 1:Ny+1] = ( A[2:Nx+2, 1:Ny+1] - 2. * A[1:Nx+1, 1:Ny+1] + A[0:Nx, 1:Ny+1] )/(dx*dx) + \ out[1:Nx+1, 1:Ny+1] = ( A[2:Nx+2, 1:Ny+1] - 2. * A[1:Nx+1, 1:Ny+1] + A[0:Nx, 1:Ny+1] )/(dx*dx) + \
( A[1:Nx+1, 2:Ny+2] - 2. * A[1:Nx+1, 1:Ny+1] + A[1:Nx+1, 0:Ny] )/(dy*dy) ( A[1:Nx+1, 2:Ny+2] - 2. * A[1:Nx+1, 1:Ny+1] + A[1:Nx+1, 0:Ny] )/(dy*dy)
return out return out
def gradient_x(A, Nx, Ny, dx): def gradient_x(A, Nx, Ny, dx):
grad_x = np.zeros((Nx+2, Ny+2)) grad_x = np.zeros((Nx+2, Ny+2))
grad_x[1:Nx+1, 1:Ny+1] = (A[2:Nx+2, 1:Ny+1] - A[0:Nx, 1:Ny+1])/dx grad_x[1:Nx+1, 1:Ny+1] = (A[2:Nx+2, 1:Ny+1] - A[0:Nx, 1:Ny+1])/dx
return grad_x return grad_x
def gradient_y(A, Nx, Ny, dy): def gradient_y(A, Nx, Ny, dy):
grad_y = np.zeros((Nx+2, Ny+2)) grad_y = np.zeros((Nx+2, Ny+2))
grad_y[1:Nx+1, 1:Ny+1] = (A[1:Nx+1, 2:Ny+2] - A[1:Nx+1, 0:Ny])/dy grad_y[1:Nx+1, 1:Ny+1] = (A[1:Nx+1, 2:Ny+2] - A[1:Nx+1, 0:Ny])/dy
return grad_y return grad_y
# Really periodic boundaries??? # Really periodic boundaries???
def periodic(A, Nx, Ny): def periodic(A, Nx, Ny):
A[0,:] = A[Nx, :] A[0,:] = A[Nx, :]
A[:,0] = A[:, Ny] A[:,0] = A[:, Ny]
A[Nx+1,:] = A[1, :] A[Nx+1,:] = A[1, :]
A[:, Ny+1] = A[:,1] A[:, Ny+1] = A[:,1]
# Simulation cell parameters # Simulation cell parameters
Nx, Ny = 300, 300 Nx, Ny = 300, 300
dx, dy = 0.03, 0.03 dx, dy = 0.03, 0.03
NxNy = Nx*Ny NxNy = Nx*Ny
# Time integration parameters # Time integration parameters
n_step = 2 n_step = 2
n_print = 50 n_print = 50
dtime = 1.0e-4 dtime = 1.0e-4
t = 0.0 t = 0.0
# Material specific parameters # Material specific parameters
tau = 0.0003 tau = 0.0003
epsilonb = 0.01 epsilonb = 0.01
kappa = 1.8 kappa = 1.8
delta = 0.02 delta = 0.02
aniso = 6.0 aniso = 6.0
alpha = 0.9 alpha = 0.9
gamma = 10.0 gamma = 10.0
t_eq = 1.0 t_eq = 1.0
theta_0 = 0.2 theta_0 = 0.2
seed = 5.0 seed = 5.0
# Prepare microstructure # Prepare microstructure
phi,T = nucleus(Nx, Ny, seed) phi,T = nucleus(Nx, Ny, seed)
#imageToVTK("dendritic-solidification_0", cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))}) #imageToVTK("dendritic-solidification_0", cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))})
# Simulate # Simulate
start = timer() start = timer()
for step in range(1, n_step): for step in range(1, n_step):
t += dtime t += dtime
# Boundary handling (periodic) # Boundary handling (periodic)
periodic(phi, Nx, Ny) periodic(phi, Nx, Ny)
periodic(T, Nx, Ny) periodic(T, Nx, Ny)
# Free energy derivative # Free energy derivative
#dfdcon = A * (2. * con * (1. - con)**2 - 2. * con**2 * (1. - con)) #dfdcon = A * (2. * con * (1. - con)**2 - 2. * con**2 * (1. - con))
# Laplacians # Laplacians
lap_phi = laplacian(phi, Nx, Ny, dx, dy) lap_phi = laplacian(phi, Nx, Ny, dx, dy)
lap_T = laplacian(T, Nx, Ny, dx, dy) lap_T = laplacian(T, Nx, Ny, dx, dy)
# Gradients # Gradients
phi_dx = gradient_x(phi, Nx, Ny, dx) phi_dx = gradient_x(phi, Nx, Ny, dx)
phi_dy = gradient_y(phi, Nx, Ny, dy) phi_dy = gradient_y(phi, Nx, Ny, dy)
# Angle # Angle
theta = np.arctan2(phi_dy, phi_dx) theta = np.arctan2(phi_dy, phi_dx)
# Epsilon and derivative # Epsilon and derivative
epsilon = epsilonb * (1.0 + delta * np.cos(aniso * (theta - theta_0))) epsilon = epsilonb * (1.0 + delta * np.cos(aniso * (theta - theta_0)))
epsilon_deriv = -epsilonb * aniso * delta * np.sin(aniso * (theta - theta_0)) epsilon_deriv = -epsilonb * aniso * delta * np.sin(aniso * (theta - theta_0))
first_term, second_term = (np.zeros((Nx+2, Ny+2)), np.zeros((Nx+2, Ny+2))) first_term, second_term = (np.zeros((Nx+2, Ny+2)), np.zeros((Nx+2, Ny+2)))
first_term[1:Nx+1, 1:Nx+1] = (epsilon[1:Nx+1, 2:Nx+2] * epsilon_deriv[1:Nx+1, 2:Nx+2] * phi_dx[1:Nx+1, 2:Nx+2]\ first_term[1:Nx+1, 1:Nx+1] = (epsilon[1:Nx+1, 2:Nx+2] * epsilon_deriv[1:Nx+1, 2:Nx+2] * phi_dx[1:Nx+1, 2:Nx+2]\
- epsilon[1:Nx+1, 0:Nx]*epsilon_deriv[1:Nx+1, 0:Nx]*phi_dx[1:Nx+1, 0:Nx])/dy - epsilon[1:Nx+1, 0:Nx]*epsilon_deriv[1:Nx+1, 0:Nx]*phi_dx[1:Nx+1, 0:Nx])/dy
second_term[1:Nx+1, 1:Nx+1] = (epsilon[2:Nx+2, 1:Nx+1] * epsilon_deriv[2:Nx+2, 1:Nx+1] * phi_dy[2:Nx+2, 1:Nx+1]\ second_term[1:Nx+1, 1:Nx+1] = (epsilon[2:Nx+2, 1:Nx+1] * epsilon_deriv[2:Nx+2, 1:Nx+1] * phi_dy[2:Nx+2, 1:Nx+1]\
- epsilon[0:Nx, 1:Nx+1]*epsilon_deriv[0:Nx, 1:Nx+1]*phi_dy[0:Nx, 1:Nx+1])/dx - epsilon[0:Nx, 1:Nx+1]*epsilon_deriv[0:Nx, 1:Nx+1]*phi_dy[0:Nx, 1:Nx+1])/dx
#first_term = gradient_y(epsilon*epsilon_deriv*phi_dx, Nx, Ny, dy) #first_term = gradient_y(epsilon*epsilon_deriv*phi_dx, Nx, Ny, dy)
#second_term = gradient_x(epsilon*epsilon_deriv*phi_dx, Nx, Ny, dx) #second_term = gradient_x(epsilon*epsilon_deriv*phi_dx, Nx, Ny, dx)
# Factor m # Factor m
m = alpha/math.pi * np.arctan(gamma*(t_eq - T)) m = alpha/math.pi * np.arctan(gamma*(t_eq - T))
# Time integration # Time integration
phi_old = phi.copy() phi_old = phi.copy()
phi = phi + (dtime/tau) * (first_term - second_term + epsilon**2 * lap_phi + phi*(1.0-phi)*(phi - 0.5 + m)) phi = phi + (dtime/tau) * (first_term - second_term + epsilon**2 * lap_phi + phi*(1.0-phi)*(phi - 0.5 + m))
T = T + dtime * lap_T + kappa*(phi - phi_old) T = T + dtime * lap_T + kappa*(phi - phi_old)
#if step % n_print == 0: #if step % n_print == 0:
#imageToVTK("dendritic-solidification_" + str(step), cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))}) #imageToVTK("dendritic-solidification_" + str(step), cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))})
end = timer() end = timer()
# Runtime is compared later # Runtime is compared later
runtime = end-start runtime = end-start
print("Done", runtime) print("Done", runtime)
``` ```
%% Output %% Output
Done 0.012862175000009302 Done 0.028609734028577805
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#4.4.4 #4.4.4
import numpy as np import numpy as np
#from pyevtk.hl import imageToVTK #from pyevtk.hl import imageToVTK
import math import math
import sympy as sp import sympy as sp
from pystencils import Field from pystencils import Field
from sympy import atan2, sin, cos, atan from sympy import atan2, sin, cos, atan
def laplacian_sympy(field, dx, dy): def laplacian_sympy(field, dx, dy):
return (field[1,0] - 2.0 * field[0,0] + field[-1,0])/(dx*dx) + (field[0,1] - 2.0*field[0,0] + field[0,-1])/(dy*dy) return (field[1,0] - 2.0 * field[0,0] + field[-1,0])/(dx*dx) + (field[0,1] - 2.0*field[0,0] + field[0,-1])/(dy*dy)
def nucleus(Nx, Ny, seed): def nucleus(Nx, Ny, seed):
T = np.zeros((Nx+2, Ny+2)) T = np.zeros((Nx+2, Ny+2))
phi = np.zeros((Nx+2, Ny+2)) phi = np.zeros((Nx+2, Ny+2))
for i in range(0,Nx+2): for i in range(0,Nx+2):
for j in range(0,Ny+2): for j in range(0,Ny+2):
if (i-(Nx+2)/2)**2 + (j-(Ny+2)/2)**2 < seed: if (i-(Nx+2)/2)**2 + (j-(Ny+2)/2)**2 < seed:
phi[i,j] = 1.0 phi[i,j] = 1.0
return (phi, T) return (phi, T)
def gradient_x_sympy(A, dx): def gradient_x_sympy(A, dx):
return (A[1,0] - A[-1,0])/dx return (A[1,0] - A[-1,0])/dx
def gradient_y_sympy(A, dy): def gradient_y_sympy(A, dy):
return (A[0,1] - A[0,-1])/dy return (A[0,1] - A[0,-1])/dy
# Really periodic boundaries??? # Really periodic boundaries???
def periodic(A, Nx, Ny): def periodic(A, Nx, Ny):
A[0,:] = A[Nx, :] A[0,:] = A[Nx, :]
A[:,0] = A[:, Ny] A[:,0] = A[:, Ny]
A[Nx+1,:] = A[1, :] A[Nx+1,:] = A[1, :]
A[:, Ny+1] = A[:,1] A[:, Ny+1] = A[:,1]
# Simulation cell parameters # Simulation cell parameters
Nx, Ny = 300, 300 Nx, Ny = 300, 300
dx, dy = 0.03, 0.03 dx, dy = 0.03, 0.03
NxNy = Nx*Ny NxNy = Nx*Ny
# Time integration parameters # Time integration parameters
n_step = 2 n_step = 2
n_print = 50 n_print = 50
dtime = 1.0e-4 dtime = 1.0e-4
t = 0.0 t = 0.0
# Material specific parameters # Material specific parameters
tau = 0.0003 tau = 0.0003
epsilonb = 0.01 epsilonb = 0.01
kappa = 1.8 kappa = 1.8
delta = 0.02 delta = 0.02
aniso = 6.0 aniso = 6.0
alpha = 0.9 alpha = 0.9
gamma = 10.0 gamma = 10.0
t_eq = 1.0 t_eq = 1.0
theta_0 = 0.2 theta_0 = 0.2
seed = 5.0 seed = 5.0
# Prepare microstructure # Prepare microstructure
phi,T = nucleus(Nx, Ny, seed) phi,T = nucleus(Nx, Ny, seed)
phi_old = phi.copy() phi_old = phi.copy()
#imageToVTK("dendritic-solidification_0", cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))}) #imageToVTK("dendritic-solidification_0", cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))})
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
size = (Nx, Ny) # domain size size = (Nx, Ny) # domain size
deriv = np.zeros((Nx+2, Ny+2)) deriv = np.zeros((Nx+2, Ny+2))
dst = np.zeros((Nx+2, Ny+2)) dst = np.zeros((Nx+2, Ny+2))
# Prepare microstructure # Prepare microstructure
tau_s, epsilonb_s, kappa_s, dx_s, dy_s, dt_s, theta_0_s, aniso_s, delta_s, alpha_s, gamma_s, t_eq_s \ tau_s, epsilonb_s, kappa_s, dx_s, dy_s, dt_s, theta_0_s, aniso_s, delta_s, alpha_s, gamma_s, t_eq_s \
= sp.symbols("tau, epsilonb, kappa, dx, dy, dt, theta_0, aniso, delta, alpha, gamma, t_eq") = sp.symbols("tau, epsilonb, kappa, dx, dy, dt, theta_0, aniso, delta, alpha, gamma, t_eq")
phi_s = Field.createFromNumpyArray("phi", phi) phi_s = Field.createFromNumpyArray("phi", phi)
phi_old_s = Field.createFromNumpyArray("phi_old", phi_old) phi_old_s = Field.createFromNumpyArray("phi_old", phi_old)
T_s = Field.createFromNumpyArray("T", T) T_s = Field.createFromNumpyArray("T", T)
eps = np.zeros((Nx+2, Ny+2)) eps = np.zeros((Nx+2, Ny+2))
eps_d = np.zeros((Nx+2, Ny+2)) eps_d = np.zeros((Nx+2, Ny+2))
phi_dx_array = np.zeros((Nx+2, Ny+2)) phi_dx_array = np.zeros((Nx+2, Ny+2))
phi_dy_array = np.zeros((Nx+2, Ny+2)) phi_dy_array = np.zeros((Nx+2, Ny+2))
theta_array = np.zeros((Nx+2, Ny+2)) theta_array = np.zeros((Nx+2, Ny+2))
epsilon = Field.createFromNumpyArray("epsilon", eps) epsilon = Field.createFromNumpyArray("epsilon", eps)
epsilon_deriv = Field.createFromNumpyArray("epsilon_deriv", eps_d) epsilon_deriv = Field.createFromNumpyArray("epsilon_deriv", eps_d)
phi_dx = Field.createFromNumpyArray("phi_dx", phi_dx_array) phi_dx = Field.createFromNumpyArray("phi_dx", phi_dx_array)
phi_dy = Field.createFromNumpyArray("phi_dy", phi_dy_array) phi_dy = Field.createFromNumpyArray("phi_dy", phi_dy_array)
theta = Field.createFromNumpyArray("theta", theta_array) theta = Field.createFromNumpyArray("theta", theta_array)
lap_phi = laplacian_sympy(phi_s, dx_s, dy_s) lap_phi = laplacian_sympy(phi_s, dx_s, dy_s)
lap_T = laplacian_sympy(T_s, dx_s, dy_s) lap_T = laplacian_sympy(T_s, dx_s, dy_s)
one = sp.Eq(phi_dx[1,1], gradient_x_sympy(phi_s, dx_s)) one = sp.Eq(phi_dx[1,1], gradient_x_sympy(phi_s, dx_s))
two = sp.Eq(phi_dy[1,1], gradient_y_sympy(phi_s, dy_s)) two = sp.Eq(phi_dy[1,1], gradient_y_sympy(phi_s, dy_s))
zero = sp.Eq(theta[0,0], atan2(phi_dy[0,0], phi_dx[0,0])) zero = sp.Eq(theta[0,0], atan2(phi_dy[0,0], phi_dx[0,0]))
three = sp.Eq(epsilon[0,0], epsilonb_s * (1.0 + delta * cos(aniso_s * (theta[0,0] - theta_0_s)))) three = sp.Eq(epsilon[0,0], epsilonb_s * (1.0 + delta * cos(aniso_s * (theta[0,0] - theta_0_s))))
four = sp.Eq(epsilon_deriv[0,0], -epsilonb_s * aniso_s * delta_s * sin(aniso_s * (theta[0,0] - theta_0_s))) four = sp.Eq(epsilon_deriv[0,0], -epsilonb_s * aniso_s * delta_s * sin(aniso_s * (theta[0,0] - theta_0_s)))
first_term = (epsilon[0,1] * epsilon_deriv[0,1] * phi_dx[0,1] - epsilon[0,-1] * epsilon_deriv[0,-1]*phi_dx[0,-1])/dy_s first_term = (epsilon[0,1] * epsilon_deriv[0,1] * phi_dx[0,1] - epsilon[0,-1] * epsilon_deriv[0,-1]*phi_dx[0,-1])/dy_s
second_term = (epsilon[1,0] * epsilon_deriv[1,0] * phi_dy[1,0] - epsilon[-1,0] * epsilon_deriv[-1,0]*phi_dy[-1,0])/dx_s second_term = (epsilon[1,0] * epsilon_deriv[1,0] * phi_dy[1,0] - epsilon[-1,0] * epsilon_deriv[-1,0]*phi_dy[-1,0])/dx_s
m = alpha_s/math.pi * atan(gamma_s * (t_eq_s - T[0,0])) m = alpha_s/math.pi * atan(gamma_s * (t_eq_s - T[0,0]))
five = sp.Eq(phi_old_s[0,0], phi_s[0,0]) five = sp.Eq(phi_old_s[0,0], phi_s[0,0])
six = sp.Eq(phi_s[0,0], phi_s[0,0] +(dt_s/tau_s) * (first_term - second_term + epsilon[0,0]**2 * lap_phi + phi[0,0] * (1.0 -phi[0,0])*(phi[0,0] - 0.5 + m))) six = sp.Eq(phi_s[0,0], phi_s[0,0] +(dt_s/tau_s) * (first_term - second_term + epsilon[0,0]**2 * lap_phi + phi[0,0] * (1.0 -phi[0,0])*(phi[0,0] - 0.5 + m)))
seven = sp.Eq(T_s[0,0], T_s[0,0] + dt_s * lap_T + kappa_s * (phi[0,0] - phi_old[0,0])) seven = sp.Eq(T_s[0,0], T_s[0,0] + dt_s * lap_T + kappa_s * (phi[0,0] - phi_old[0,0]))
parameters = [(tau_s, tau), (epsilonb_s, epsilonb), (kappa_s, kappa), (dx_s, dx), (dy_s, dy),\ parameters = [(tau_s, tau), (epsilonb_s, epsilonb), (kappa_s, kappa), (dx_s, dx), (dy_s, dy),\
(dt_s, dtime), (theta_0_s, theta_0), (aniso_s, aniso), (delta_s, delta), (alpha_s, alpha),\ (dt_s, dtime), (theta_0_s, theta_0), (aniso_s, aniso), (delta_s, delta), (alpha_s, alpha),\
(gamma_s, gamma), (t_eq_s, t_eq)] (gamma_s, gamma), (t_eq_s, t_eq)]
#one = one.subs(parameters) #one = one.subs(parameters)
#two = two.subs(parameters) #two = two.subs(parameters)
#zero = zero.subs(parameters) #zero = zero.subs(parameters)
#three = three.subs(parameters) #three = three.subs(parameters)
#four = four.subs(parameters) #four = four.subs(parameters)
#five = five.subs(parameters) #five = five.subs(parameters)
#six = six.subs(parameters) #six = six.subs(parameters)
#seven = seven.subs(parameters) #seven = seven.subs(parameters)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from pystencils.llvm import createKernel, make_python_function from pystencils.llvm import createKernel, makePythonFunction
ast = createKernel([one, two, zero, three, four, five, six, seven]) ast = createKernel([one, two, zero, three, four, five, six, seven])
print(ast) print(ast)
``` ```
%% Output %% Output
KernelFunction kernel([<double * RESTRICT fd_T>, <double * RESTRICT fd_epsilon>, <double * RESTRICT fd_epsilon_deriv>, <double * RESTRICT fd_phi>, <double * RESTRICT fd_phi_dx>, <double * RESTRICT fd_phi_dy>, <double * RESTRICT fd_phi_old>, <double * RESTRICT fd_theta>, <double aniso>, <double delta>, <double dt>, <double dx>, <double dy>, <double epsilonb>, <double tau>, <double theta_0>]) KernelFunction kernel([<double * RESTRICT fd_T>, <double * RESTRICT fd_epsilon>, <double * RESTRICT fd_epsilon_deriv>, <double * RESTRICT fd_phi>, <double * RESTRICT fd_phi_dx>, <double * RESTRICT fd_phi_dy>, <double * RESTRICT fd_phi_old>, <double * RESTRICT fd_theta>, <double aniso>, <double delta>, <double dt>, <double dx>, <double dy>, <double epsilonb>, <double tau>, <double theta_0>])
for(ctr_0=1; ctr_0<301; ctr_0+=1) Block for(ctr_0=1; ctr_0<301; ctr_0+=1)
fd_phi_dx_E = pointerArithmeticFunc(fd_phi_dx, 302*ctr_0 + 302) Block fd_phi_dx_E = pointerArithmeticFunc(fd_phi_dx, 302*ctr_0 + 302)
fd_phi_W = pointerArithmeticFunc(fd_phi, 302*ctr_0 - 302) fd_phi_W = pointerArithmeticFunc(fd_phi, 302*ctr_0 - 302)
fd_phi_E = pointerArithmeticFunc(fd_phi, 302*ctr_0 + 302) fd_phi_E = pointerArithmeticFunc(fd_phi, 302*ctr_0 + 302)
fd_phi_dy_E = pointerArithmeticFunc(fd_phi_dy, 302*ctr_0 + 302) fd_phi_dy_E = pointerArithmeticFunc(fd_phi_dy, 302*ctr_0 + 302)
fd_phi_C = pointerArithmeticFunc(fd_phi, 302*ctr_0) fd_phi_C = pointerArithmeticFunc(fd_phi, 302*ctr_0)
fd_theta_C = pointerArithmeticFunc(fd_theta, 302*ctr_0) fd_theta_C = pointerArithmeticFunc(fd_theta, 302*ctr_0)
fd_phi_dy_C = pointerArithmeticFunc(fd_phi_dy, 302*ctr_0) fd_phi_dy_C = pointerArithmeticFunc(fd_phi_dy, 302*ctr_0)
fd_phi_dx_C = pointerArithmeticFunc(fd_phi_dx, 302*ctr_0) fd_phi_dx_C = pointerArithmeticFunc(fd_phi_dx, 302*ctr_0)
fd_epsilon_C = pointerArithmeticFunc(fd_epsilon, 302*ctr_0) fd_epsilon_C = pointerArithmeticFunc(fd_epsilon, 302*ctr_0)
fd_epsilon_deriv_C = pointerArithmeticFunc(fd_epsilon_deriv, 302*ctr_0) fd_epsilon_deriv_C = pointerArithmeticFunc(fd_epsilon_deriv, 302*ctr_0)
fd_phi_old_C = pointerArithmeticFunc(fd_phi_old, 302*ctr_0) fd_phi_old_C = pointerArithmeticFunc(fd_phi_old, 302*ctr_0)
fd_epsilon_E = pointerArithmeticFunc(fd_epsilon, 302*ctr_0 + 302) fd_epsilon_E = pointerArithmeticFunc(fd_epsilon, 302*ctr_0 + 302)
fd_epsilon_deriv_E = pointerArithmeticFunc(fd_epsilon_deriv, 302*ctr_0 + 302) fd_epsilon_deriv_E = pointerArithmeticFunc(fd_epsilon_deriv, 302*ctr_0 + 302)
fd_epsilon_W = pointerArithmeticFunc(fd_epsilon, 302*ctr_0 - 302) fd_epsilon_W = pointerArithmeticFunc(fd_epsilon, 302*ctr_0 - 302)
fd_epsilon_deriv_W = pointerArithmeticFunc(fd_epsilon_deriv, 302*ctr_0 - 302) fd_epsilon_deriv_W = pointerArithmeticFunc(fd_epsilon_deriv, 302*ctr_0 - 302)
fd_phi_dy_W = pointerArithmeticFunc(fd_phi_dy, 302*ctr_0 - 302) fd_phi_dy_W = pointerArithmeticFunc(fd_phi_dy, 302*ctr_0 - 302)
fd_T_C = pointerArithmeticFunc(fd_T, 302*ctr_0) fd_T_C = pointerArithmeticFunc(fd_T, 302*ctr_0)
fd_T_E = pointerArithmeticFunc(fd_T, 302*ctr_0 + 302) fd_T_E = pointerArithmeticFunc(fd_T, 302*ctr_0 + 302)
fd_T_W = pointerArithmeticFunc(fd_T, 302*ctr_0 - 302) fd_T_W = pointerArithmeticFunc(fd_T, 302*ctr_0 - 302)
for(ctr_1=1; ctr_1<301; ctr_1+=1) for(ctr_1=1; ctr_1<301; ctr_1+=1)
fd_phi_dx_E[ctr_1 + 1] = dx**castFunc(-1, double)*(fd_phi_E[ctr_1] + fd_phi_W[ctr_1]*castFunc(-1, double)) Block fd_phi_dx_E[ctr_1 + 1] = dx**castFunc(-1, double)*(castFunc(-1, double)*fd_phi_W[ctr_1] + fd_phi_E[ctr_1])
fd_phi_dy_E[ctr_1 + 1] = dy**castFunc(-1, double)*(fd_phi_C[ctr_1 + 1] + fd_phi_C[ctr_1 - 1]*castFunc(-1, double)) fd_phi_dy_E[ctr_1 + 1] = dy**castFunc(-1, double)*(castFunc(-1, double)*fd_phi_C[ctr_1 - 1] + fd_phi_C[ctr_1 + 1])
fd_theta_C[ctr_1] = atan2(fd_phi_dy_C[ctr_1], fd_phi_dx_C[ctr_1]) fd_theta_C[ctr_1] = atan2(fd_phi_dy_C[ctr_1], fd_phi_dx_C[ctr_1])
fd_epsilon_C[ctr_1] = epsilonb*(1.0 + 0.02*cos(aniso*(theta_0*castFunc(-1, double) + fd_theta_C[ctr_1]))) fd_epsilon_C[ctr_1] = epsilonb*(0.02*cos(aniso*(theta_0*castFunc(-1, double) + fd_theta_C[ctr_1])) + 1.0)
fd_epsilon_deriv_C[ctr_1] = aniso*delta*epsilonb*castFunc(-1, double)*sin(aniso*(theta_0*castFunc(-1, double) + fd_theta_C[ctr_1])) fd_epsilon_deriv_C[ctr_1] = aniso*delta*epsilonb*castFunc(-1, double)*sin(aniso*(theta_0*castFunc(-1, double) + fd_theta_C[ctr_1]))
fd_phi_old_C[ctr_1] = fd_phi_C[ctr_1] fd_phi_old_C[ctr_1] = fd_phi_C[ctr_1]
fd_phi_C[ctr_1] = dt*tau**castFunc(-1, double)*(dy**castFunc(-1, double)*(fd_epsilon_C[ctr_1 + 1]*fd_epsilon_deriv_C[ctr_1 + 1]*fd_phi_dx_C[ctr_1 + 1] + fd_epsilon_C[ctr_1 - 1]*fd_epsilon_deriv_C[ctr_1 - 1]*fd_phi_dx_C[ctr_1 - 1]*castFunc(-1, double)) + castFunc(-1, double)*dx**castFunc(-1, double)*(fd_epsilon_E[ctr_1]*fd_epsilon_deriv_E[ctr_1]*fd_phi_dy_E[ctr_1] + fd_epsilon_W[ctr_1]*fd_epsilon_deriv_W[ctr_1]*fd_phi_dy_W[ctr_1]*castFunc(-1, double)) + fd_epsilon_C[ctr_1]**castFunc(2, double)*((fd_phi_C[ctr_1 + 1] + fd_phi_C[ctr_1 - 1] - 2.0*fd_phi_C[ctr_1])*dy**castFunc(-2, double) + (-2.0*fd_phi_C[ctr_1] + fd_phi_E[ctr_1] + fd_phi_W[ctr_1])*dx**castFunc(-2, double))) + fd_phi_C[ctr_1] fd_phi_C[ctr_1] = dt*tau**castFunc(-1, double)*(dx**castFunc(-1, double)*(castFunc(-1, double)*fd_epsilon_W[ctr_1]*fd_epsilon_deriv_W[ctr_1]*fd_phi_dy_W[ctr_1] + fd_epsilon_E[ctr_1]*fd_epsilon_deriv_E[ctr_1]*fd_phi_dy_E[ctr_1])*castFunc(-1, double) + dy**castFunc(-1, double)*(castFunc(-1, double)*fd_epsilon_C[ctr_1 - 1]*fd_epsilon_deriv_C[ctr_1 - 1]*fd_phi_dx_C[ctr_1 - 1] + fd_epsilon_C[ctr_1 + 1]*fd_epsilon_deriv_C[ctr_1 + 1]*fd_phi_dx_C[ctr_1 + 1]) + (dx**castFunc(-2, double)*(-2.0*fd_phi_C[ctr_1] + fd_phi_E[ctr_1] + fd_phi_W[ctr_1]) + dy**castFunc(-2, double)*(fd_phi_C[ctr_1 + 1] + fd_phi_C[ctr_1 - 1] - 2.0*fd_phi_C[ctr_1]))*fd_epsilon_C[ctr_1]**castFunc(2, double)) + fd_phi_C[ctr_1]
fd_T_C[ctr_1] = dt*((fd_T_C[ctr_1 + 1] + fd_T_C[ctr_1 - 1] - 2.0*fd_T_C[ctr_1])*dy**castFunc(-2, double) + (-2.0*fd_T_C[ctr_1] + fd_T_E[ctr_1] + fd_T_W[ctr_1])*dx**castFunc(-2, double)) + fd_T_C[ctr_1] fd_T_C[ctr_1] = dt*(dx**castFunc(-2, double)*(-2.0*fd_T_C[ctr_1] + fd_T_E[ctr_1] + fd_T_W[ctr_1]) + dy**castFunc(-2, double)*(fd_T_C[ctr_1 + 1] + fd_T_C[ctr_1 - 1] - 2.0*fd_T_C[ctr_1])) + fd_T_C[ctr_1]
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
test = sp.atan2(1,2) test = sp.atan2(1,2)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
kernel = make_python_function(ast) kernel = makePythonFunction(ast)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Simulate # Simulate
start = timer() start = timer()
for step in range(1, n_step): for step in range(1, n_step):
t += dtime t += dtime
# Boundary handling (periodic) # Boundary handling (periodic)
periodic(phi, Nx, Ny) periodic(phi, Nx, Ny)
periodic(T, Nx, Ny) periodic(T, Nx, Ny)
#kernel(T=T, phi=phi, phi_old=phi_old, epsilon=eps, epsilon_deriv = eps_d, phi_dx=phi_dx_array, phi_dy = phi_dy_array, theta = theta_array) #kernel(T=T, phi=phi, phi_old=phi_old, epsilon=eps, epsilon_deriv = eps_d, phi_dx=phi_dx_array, phi_dy = phi_dy_array, theta = theta_array)
kernel(T=T, phi=phi, phi_old=phi_old, epsilon=eps, epsilon_deriv = eps_d, phi_dx=phi_dx_array, phi_dy = phi_dy_array, theta = theta_array,\ kernel(T=T, phi=phi, phi_old=phi_old, epsilon=eps, epsilon_deriv = eps_d, phi_dx=phi_dx_array, phi_dy = phi_dy_array, theta = theta_array,\
tau=tau, epsilonb=epsilonb, kappa=kappa, dx=dx, dy=dy,\ tau=tau, epsilonb=epsilonb, kappa=kappa, dx=dx, dy=dy,\
dt=dtime, theta_0=theta_0, aniso=aniso, delta=delta, alpha=alpha,\ dt=dtime, theta_0=theta_0, aniso=aniso, delta=delta, alpha=alpha,\
gamma_s=gamma, t_eq_s=t_eq) gamma_s=gamma, t_eq_s=t_eq)
#if step % n_print == 0: #if step % n_print == 0:
#imageToVTK("dendritic-solidification_" + str(step), cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))}) #imageToVTK("dendritic-solidification_" + str(step), cellData = {"phi" : phi.reshape((Nx+2, Nx+2, 1)), "T": T.reshape((Nx+2, Nx+2, 1))})
end = timer() end = timer()
# Runtime is compared later # Runtime is compared later
runtime = end-start runtime = end-start
print("Done", runtime) print("Done", runtime)
``` ```
%% Output %% Output
Done 0.030210711999643536 Done 0.04964133305475116
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
......
doc/modules/pyphasefield/demo_heat_equation.ipynb
+ 34
− 22
View file @ 9de3da34
This diff is collapsed.
doc/modules/pystencils/demo_llvm_polly.ipynb
+ 47
− 51
View file @ 9de3da34
This diff is collapsed.
doc/modules/pystencils/kerncraft/3d-7pt.c deleted 100644 → 0
+ 0
− 10
View file @ ae4bbde0
double a[M][N][N];
double b[M][N][N];
double s;
for(int k=1; k<M-1; ++k)
for(int j=1; j<N-1; ++j)
for(int i=1; i<N-1; ++i)
b[k][j][i] = ( a[k][j][i-1] + a[k][j][i+1]
+ a[k][j-1][i] + a[k][j+1][i]
+ a[k-1][j][i] + a[k+1][j][i]) * s;
jupyterhub.Dockerfile
+ 1
− 0
View file @ 9de3da34
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
FROM i10git.cs.fau.de:5005/software/pystencils FROM i10git.cs.fau.de:5005/software/pystencils
RUN conda install --yes -c conda-forge nodejs configurable-http-proxy RUN conda install --yes -c conda-forge nodejs configurable-http-proxy
RUN pip install oauthenticator RUN pip install oauthenticator
ADD doc/jupyterhub/jupyterhub_config.py /etc/jupyterhub/jupyterhub_config.py ADD doc/jupyterhub/jupyterhub_config.py /etc/jupyterhub/jupyterhub_config.py
......
lbmpy/session.py
+ 3
− 10
View file @ 9de3da34
import sympy as sp import sympy as sp
import numpy as np import numpy as np
from pystencils.jupytersetup import *
from lbmpy.scenarios import * from lbmpy.scenarios import *
from lbmpy.creationfunctions import * from lbmpy.creationfunctions import *
from pystencils import makeSlice, showCode from pystencils import makeSlice, showCode
...@@ -7,13 +8,5 @@ from lbmpy.boundaries import * ...@@ -7,13 +8,5 @@ from lbmpy.boundaries import *
from lbmpy.postprocessing import * from lbmpy.postprocessing import *
from lbmpy.lbstep import LatticeBoltzmannStep from lbmpy.lbstep import LatticeBoltzmannStep
from lbmpy.geometry import * from lbmpy.geometry import *
from lbmpy.parameterization import ScalingWidget
try: import lbmpy.plot2d as plt
from IPython import get_ipython
ipython = get_ipython()
from pystencils.jupytersetup import *
from lbmpy.parameterization import ScalingWidget
import lbmpy.plot2d as plt
except ImportError:
pass
lbmpy_benchmarks/phasefield/test_analytical_NPhase_full.ipynb
+ 41
− 76
View file @ 9de3da34
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from lbmpy.session import * from lbmpy.session import *
from lbmpy.phasefield.analytical import * from lbmpy.phasefield.analytical import *
from lbmpy.chapman_enskog.derivative import * from pystencils.derivative import evaluateDiffs
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Full analytical checks for N phase model # Full analytical checks for N phase model
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
x = sp.Symbol("x") x = sp.Symbol("x")
c_a = analyticInterfaceProfile(x) c_a = analyticInterfaceProfile(x)
def gamma(i,j): def gamma(i,j):
if i == j: if i == j:
return 0 return 0
elif i < j: elif i < j:
return sp.Symbol("gamma_%d%d" % (i,j)) return sp.Symbol("gamma_%d%d" % (i,j))
else: else:
return sp.Symbol("gamma_%d%d" % (j,i)) return sp.Symbol("gamma_%d%d" % (j,i))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
The following checks if the 'tanh' shaped profile is a solution to $\mu_i = 0$ and if the excess free energy is the surface tension parameter. The following checks if the 'tanh' shaped profile is a solution to $\mu_i = 0$ and if the excess free energy is the surface tension parameter.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
numPhases = 5 numPhases = 3
for numPhases in range(2, 6):
print("Checking N =", numPhases)
c = sp.symbols("c_:%d" % (numPhases-1,))
F = freeEnergyFunctionalNPhases(orderParameters=c, surfaceTensions=gamma)
μ = chemicalPotentialsFromFreeEnergy(F, c)
lastPhaseIdx = numPhases-1
# Check all permutations of phases
for i in range(numPhases):
for j in range(numPhases):
if i == j:
continue
print(" -> Testing interface between", i, "and", j)
substitutions = {c_i: 0 for c_i in c}
if i != lastPhaseIdx:
substitutions[c[i]] = c_a
if j != lastPhaseIdx:
substitutions[c[j]] = 1 - c_a
for μ_i in μ:
res = evaluateDiffs(μ_i.subs(substitutions), x).expand()
assert res == 0, "Analytic interface profile wrong for phase between %d and %d" % (i,j)
twoPhaseFreeEnergy = F.subs(substitutions)
twoPhaseFreeEnergy = sp.simplify(evaluateDiffs(twoPhaseFreeEnergy, x))
result = coshIntegral(twoPhaseFreeEnergy, x)
assert result == gamma(i,j)
print("Checking N =", numPhases)
c = sp.symbols("c_:%d" % (numPhases-1,))
F = freeEnergyFunctionalNPhases(orderParameters=c, surfaceTensions=gamma)
μ = chemicalPotentialsFromFreeEnergy(F, c)
lastPhaseIdx = numPhases-1
# Check all permutations of phases
for i in range(numPhases):
for j in range(numPhases):
if i == j:
continue
print(" -> Testing interface between", i, "and", j)
substitutions = {c_i: 0 for c_i in c}
if i != lastPhaseIdx:
substitutions[c[i]] = c_a
if j != lastPhaseIdx:
substitutions[c[j]] = 1 - c_a
for μ_i in μ:
res = evaluateDiffs(μ_i.subs(substitutions), x).expand()
assert res == 0, "Analytic interface profile wrong for phase between %d and %d" % (i,j)
twoPhaseFreeEnergy = F.subs(substitutions)
twoPhaseFreeEnergy = sp.simplify(evaluateDiffs(twoPhaseFreeEnergy, x))
result = coshIntegral(twoPhaseFreeEnergy, x)
assert result == gamma(i,j)
``` ```
%% Output %% Output
Checking N = 2
-> Testing interface between 0 and 1
-> Testing interface between 1 and 0
Checking N = 3 Checking N = 3
-> Testing interface between 0 and 1 -> Testing interface between 0 and 1
-> Testing interface between 0 and 2 -> Testing interface between 0 and 2
-> Testing interface between 1 and 0 -> Testing interface between 1 and 0
-> Testing interface between 1 and 2 -> Testing interface between 1 and 2
-> Testing interface between 2 and 0 -> Testing interface between 2 and 0
-> Testing interface between 2 and 1 -> Testing interface between 2 and 1
Checking N = 4
-> Testing interface between 0 and 1 %% Cell type:code id: tags:
-> Testing interface between 0 and 2
-> Testing interface between 0 and 3 ``` python
-> Testing interface between 1 and 0 ```
-> Testing interface between 1 and 2
-> Testing interface between 1 and 3
-> Testing interface between 2 and 0
-> Testing interface between 2 and 1
-> Testing interface between 2 and 3
-> Testing interface between 3 and 0
-> Testing interface between 3 and 1
-> Testing interface between 3 and 2
Checking N = 5
-> Testing interface between 0 and 1
-> Testing interface between 0 and 2
-> Testing interface between 0 and 3
-> Testing interface between 0 and 4
-> Testing interface between 1 and 0
-> Testing interface between 1 and 2
-> Testing interface between 1 and 3
-> Testing interface between 1 and 4
-> Testing interface between 2 and 0
-> Testing interface between 2 and 1
-> Testing interface between 2 and 3
-> Testing interface between 2 and 4
-> Testing interface between 3 and 0
-> Testing interface between 3 and 1
-> Testing interface between 3 and 2
-> Testing interface between 3 and 4
-> Testing interface between 4 and 0
-> Testing interface between 4 and 1
-> Testing interface between 4 and 2
-> Testing interface between 4 and 3
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment