From cb8c344a68d0e48da17c64b1f74378ae6f0c4c46 Mon Sep 17 00:00:00 2001 From: nfarabullini Date: Wed, 13 Mar 2024 14:15:09 +0100 Subject: [PATCH] Removed md files --- .../exercises/1_simple_addition_solution.md | 85 ------ .../next/exercises/3_gradient_exercise.md | 120 -------- .../exercises/3_gradient_exercise_solution.md | 119 -------- docs/user/next/exercises/4_curl_exercise.md | 121 -------- .../exercises/4_curl_exercise_solution.md | 126 -------- .../exercises/5_vector_laplace_exercise.md | 264 ---------------- .../5_vector_laplace_exercise_solution.md | 286 ------------------ docs/user/next/exercises/6_where_domain.md | 163 ---------- .../exercises/6_where_domain_solutions.md | 164 ---------- docs/user/next/exercises/7_scan_operator.md | 189 ------------ .../exercises/7_scan_operator_solutions.md | 191 ------------ .../8_diffusion_exercise_solution.md | 159 ---------- 12 files changed, 1987 deletions(-) delete mode 100644 docs/user/next/exercises/1_simple_addition_solution.md delete mode 100644 docs/user/next/exercises/3_gradient_exercise.md delete mode 100644 docs/user/next/exercises/3_gradient_exercise_solution.md delete mode 100644 docs/user/next/exercises/4_curl_exercise.md delete mode 100644 docs/user/next/exercises/4_curl_exercise_solution.md delete mode 100644 docs/user/next/exercises/5_vector_laplace_exercise.md delete mode 100644 docs/user/next/exercises/5_vector_laplace_exercise_solution.md delete mode 100644 docs/user/next/exercises/6_where_domain.md delete mode 100644 docs/user/next/exercises/6_where_domain_solutions.md delete mode 100644 docs/user/next/exercises/7_scan_operator.md delete mode 100644 docs/user/next/exercises/7_scan_operator_solutions.md delete mode 100644 docs/user/next/exercises/8_diffusion_exercise_solution.md diff --git a/docs/user/next/exercises/1_simple_addition_solution.md b/docs/user/next/exercises/1_simple_addition_solution.md deleted file mode 100644 index ee3b0f9d17..0000000000 --- a/docs/user/next/exercises/1_simple_addition_solution.md +++ /dev/null @@ -1,85 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# 1. Simple Addition - -```{code-cell} ipython3 -import numpy as np -import matplotlib.pyplot as plt - -import gt4py.next as gtx -from gt4py.next import Dims -``` - -Next we implement the stencil and a numpy reference version, in order to verify them against each other. - -```{code-cell} ipython3 -I = gtx.Dimension("I") -J = gtx.Dimension("J") -size = 10 -``` - -```{code-cell} ipython3 -def addition_numpy(a: np.array, b: np.array) -> np.array: - c = a + b - return c -``` - -```{code-cell} ipython3 -@gtx.field_operator -def addition( - a: gtx.Field[Dims[I, J], float], b: gtx.Field[Dims[I, J], float] -) -> gtx.Field[Dims[I, J], float]: - return a + b -``` - -```{code-cell} ipython3 -def test_addition(backend=None): - domain = gtx.domain({I: size, J: size}) - - a_data = np.fromfunction(lambda xx, yy: xx, domain.shape, dtype=float) - a = gtx.as_field(domain, a_data, allocator=backend) - b_data = np.fromfunction(lambda xx, yy: yy, domain.shape, dtype=float) - b = gtx.as_field(domain, b_data, allocator=backend) - - c_numpy = addition_numpy(a.asnumpy(), b.asnumpy()) - - c = gtx.zeros(domain, allocator=backend) - - addition(a, b, out=c, offset_provider={}) - - assert np.allclose(c.asnumpy(), c_numpy) - - print("Result:") - print(c) - print(c.asnumpy()) - - # Plots - fig, ax = plt.subplot_mosaic([ - ['a', 'b', 'c'] - ]) - ax['a'].imshow(a.asnumpy()) - ax['b'].imshow(b.asnumpy()) - ax['c'].imshow(c.asnumpy()) - - print("\nTest successful!") -``` - -```{code-cell} ipython3 -test_addition() -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/3_gradient_exercise.md b/docs/user/next/exercises/3_gradient_exercise.md deleted file mode 100644 index 21abcf35a0..0000000000 --- a/docs/user/next/exercises/3_gradient_exercise.md +++ /dev/null @@ -1,120 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# 4. Gradient - -+++ - -Another example is the gradient defined at the center of a Cell $\mathbf{P}$ of a scalar function $f$. We approximate this by taking the sum over the three edges and multiplying $f(e)$ with the edge normal $\mathbf{n}_e$ and the edge length $L_e$ and dividing the resulting sum with the cell area $A_P$. -The result will be the two components of the gradient vector. - -![](../images/gradient_picture.png "Divergence") - -![](../images/gradient_formula.png "Divergence") - -```{code-cell} ipython3 -from helpers import * - - -import gt4py.next as gtx -``` - -```{code-cell} ipython3 -def gradient_numpy( - c2e: np.array, - f: np.array, - nx: np.array, - ny: np.array, - L: np.array, - A: np.array, - edge_orientation: np.array, -) -> gtx.tuple[np.array, np.array]: - # edge_orientation = np.expand_dims(edge_orientation, axis=-1) - f_x = np.sum(f[c2e] * nx[c2e] * L[c2e] * edge_orientation, axis=1) / A - f_y = np.sum(f[c2e] * ny[c2e] * L[c2e] * edge_orientation, axis=1) / A - return f_x, f_y -``` - -```{code-cell} ipython3 -@gtx.field_operator -def gradient( - f: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - L: gtx.Field[Dims[E], float], - A: gtx.Field[Dims[C], float], - edge_orientation: gtx.Field[Dims[C, C2EDim], float], -) -> gtx.tuple[gtx.Field[Dims[C], float], gtx.Field[Dims[C], float]]: - # TODO: fix components of gradient - f_x = A - f_y = A - return f_x, f_y -``` - -```{code-cell} ipython3 -def test_gradient(): - backend = None - # backend = gtfn_cpu - # backend = gtfn_gpu - - cell_domain = gtx.domain({C: n_cells}) - edge_domain = gtx.domain({E: n_edges}) - - f = random_field(edge_domain, allocator=backend) - nx = random_field(edge_domain, allocator=backend) - ny = random_field(edge_domain, allocator=backend) - L = random_field(edge_domain, allocator=backend) - A = random_field(cell_domain, allocator=backend) - edge_orientation = random_sign( - gtx.domain({C: n_cells, C2EDim: 3}), allocator=backend - ) - - gradient_numpy_x, gradient_numpy_y = gradient_numpy( - c2e_table, - f.asnumpy(), - nx.asnumpy(), - ny.asnumpy(), - L.asnumpy(), - A.asnumpy(), - edge_orientation.asnumpy(), - ) - - c2e_connectivity = gtx.NeighborTableOffsetProvider(c2e_table, C, E, 3, has_skip_values=False) - - gradient_gt4py_x = gtx.zeros(cell_domain, allocator=backend) - gradient_gt4py_y = gtx.zeros(cell_domain, allocator=backend) - - gradient( - f, - nx, - ny, - L, - A, - edge_orientation, - out=(gradient_gt4py_x, gradient_gt4py_y), - offset_provider={C2E.value: c2e_connectivity}, - ) - - assert np.allclose(gradient_gt4py_x.asnumpy(), gradient_numpy_x) - assert np.allclose(gradient_gt4py_y.asnumpy(), gradient_numpy_y) -``` - -```{code-cell} ipython3 -test_gradient() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/3_gradient_exercise_solution.md b/docs/user/next/exercises/3_gradient_exercise_solution.md deleted file mode 100644 index fb3080440e..0000000000 --- a/docs/user/next/exercises/3_gradient_exercise_solution.md +++ /dev/null @@ -1,119 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# 4. Gradient - -+++ - -Another example is the gradient defined at the center of a Cell $\mathbf{P}$ of a scalar function $f$. We approximate this by taking the sum over the three edges and multiplying $f(e)$ with the edge normal $\mathbf{n}_e$ and the edge length $L_e$ and dividing the resulting sum with the cell area $A_P$. -The result will be the two components of the gradient vector. - -![](../images/gradient_picture.png "Divergence") - -![](../images/gradient_formula.png "Divergence") - -```{code-cell} ipython3 -from helpers import * - - -import gt4py.next as gtx -``` - -```{code-cell} ipython3 -def gradient_numpy( - c2e: np.array, - f: np.array, - nx: np.array, - ny: np.array, - L: np.array, - A: np.array, - edge_orientation: np.array, -) -> gtx.tuple[np.array, np.array]: - # edge_orientation = np.expand_dims(edge_orientation, axis=-1) - f_x = np.sum(f[c2e] * nx[c2e] * L[c2e] * edge_orientation, axis=1) / A - f_y = np.sum(f[c2e] * ny[c2e] * L[c2e] * edge_orientation, axis=1) / A - return f_x, f_y -``` - -```{code-cell} ipython3 -@gtx.field_operator -def gradient( - f: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - L: gtx.Field[Dims[E], float], - A: gtx.Field[Dims[C], float], - edge_orientation: gtx.Field[Dims[C, C2EDim], float], -) -> gtx.tuple[gtx.Field[Dims[C], float], gtx.Field[Dims[C], float]]: - f_x = neighbor_sum(f(C2E) * nx(C2E) * L(C2E) * edge_orientation, axis=C2EDim) / A - f_y = neighbor_sum(f(C2E) * ny(C2E) * L(C2E) * edge_orientation, axis=C2EDim) / A - return f_x, f_y -``` - -```{code-cell} ipython3 -def test_gradient(): - backend = None - # backend = gtfn_cpu - # backend = gtfn_gpu - - cell_domain = gtx.domain({C: n_cells}) - edge_domain = gtx.domain({E: n_edges}) - - f = random_field(edge_domain, allocator=backend) - nx = random_field(edge_domain, allocator=backend) - ny = random_field(edge_domain, allocator=backend) - L = random_field(edge_domain, allocator=backend) - A = random_field(cell_domain, allocator=backend) - edge_orientation = random_sign( - gtx.domain({C: n_cells, C2EDim: 3}), allocator=backend - ) - - gradient_numpy_x, gradient_numpy_y = gradient_numpy( - c2e_table, - f.asnumpy(), - nx.asnumpy(), - ny.asnumpy(), - L.asnumpy(), - A.asnumpy(), - edge_orientation.asnumpy(), - ) - - c2e_connectivity = gtx.NeighborTableOffsetProvider(c2e_table, C, E, 3, has_skip_values=False) - - gradient_gt4py_x = gtx.zeros(cell_domain, allocator=backend) - gradient_gt4py_y = gtx.zeros(cell_domain, allocator=backend) - - gradient( - f, - nx, - ny, - L, - A, - edge_orientation, - out=(gradient_gt4py_x, gradient_gt4py_y), - offset_provider={C2E.value: c2e_connectivity}, - ) - - assert np.allclose(gradient_gt4py_x.asnumpy(), gradient_numpy_x) - assert np.allclose(gradient_gt4py_y.asnumpy(), gradient_numpy_y) -``` - -```{code-cell} ipython3 -test_gradient() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/4_curl_exercise.md b/docs/user/next/exercises/4_curl_exercise.md deleted file mode 100644 index b86d365255..0000000000 --- a/docs/user/next/exercises/4_curl_exercise.md +++ /dev/null @@ -1,121 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# 5. Curl - -+++ - -As the last example of the easier operations, we take a look at the curl of a vector field $\mathbf{v}$ defined at a vertex $\mathbf{N}$. -To approximate this, we once again iterate over all of the direct neighboring edges of the vertex in the center and for each edge take the dot product of the vector field $\mathbf{v}_e$ with the edge normals $\mathbf{n}_f$ and multiply that by the dual edge length $\hat{L}_e$. The resulting neighbor sum is then divided by the dual area $\hat{A}_N$, which is the area of the Voronoi cell around the Vertex $\mathbf{N}$. - - -![](../images/curl_picture.png "Divergence") - - -![](../images/curl_formula.png "Divergence") - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx -``` - -```{code-cell} ipython3 -def curl_numpy( - v2e: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - dualL: np.array, - dualA: np.array, - edge_orientation: np.array, -) -> np.array: - uv_curl = ( - np.sum( - (u[v2e] * nx[v2e] + v[v2e] * ny[v2e]) * dualL[v2e] * edge_orientation, - axis=1, - ) - / dualA - ) - - return uv_curl -``` - -```{code-cell} ipython3 -@gtx.field_operator -def curl( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - dualL: gtx.Field[Dims[E], float], - dualA: gtx.Field[Dims[V], float], - edge_orientation: gtx.Field[Dims[V, V2EDim], float], -) -> gtx.Field[Dims[V], float]: - # TODO: fix curl - uv_curl = dualA - - return uv_curl -``` - -```{code-cell} ipython3 -def test_curl(): - backend = None - # backend = gtfn_cpu - # backend = gtfn_gpu - - edge_domain = gtx.domain({E: n_edges}) - vertex_domain = gtx.domain({V: n_vertices}) - - u = random_field(edge_domain, allocator=backend) - v = random_field(edge_domain, allocator=backend) - nx = random_field(edge_domain, allocator=backend) - ny = random_field(edge_domain, allocator=backend) - dualL = random_field(edge_domain, allocator=backend) - dualA = random_field(vertex_domain, allocator=backend) - edge_orientation = random_sign( - gtx.domain({V: n_vertices, V2EDim: 6}), allocator=backend - ) - - divergence_ref = curl_numpy( - v2e_table, - u.asnumpy(), - v.asnumpy(), - nx.asnumpy(), - ny.asnumpy(), - dualL.asnumpy(), - dualA.asnumpy(), - edge_orientation.asnumpy(), - ) - - v2e_connectivity = gtx.NeighborTableOffsetProvider(v2e_table, V, E, 6, has_skip_values=False) - - curl_gt4py = gtx.zeros(vertex_domain, allocator=backend) - - curl( - u, v, nx, ny, dualL, dualA, edge_orientation, out = curl_gt4py, offset_provider = {V2E.value: v2e_connectivity} - ) - - assert np.allclose(curl_gt4py.asnumpy(), divergence_ref) -``` - -```{code-cell} ipython3 -test_curl() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/4_curl_exercise_solution.md b/docs/user/next/exercises/4_curl_exercise_solution.md deleted file mode 100644 index 72e3ccf728..0000000000 --- a/docs/user/next/exercises/4_curl_exercise_solution.md +++ /dev/null @@ -1,126 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# 5. Curl - -+++ - -As the last example of the easier operations, we take a look at the curl of a vector field $\mathbf{v}$ defined at a vertex $\mathbf{N}$. -To approximate this, we once again iterate over all of the direct neighboring edges of the vertex in the center and for each edge take the dot product of the vector field $\mathbf{v}_e$ with the edge normals $\mathbf{n}_f$ and multiply that by the dual edge length $\hat{L}_e$. The resulting neighbor sum is then divided by the dual area $\hat{A}_N$, which is the area of the Voronoi cell around the Vertex $\mathbf{N}$. - - -![](../images/curl_picture.png "Divergence") - - -![](../images/curl_formula.png "Divergence") - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx -``` - -```{code-cell} ipython3 -def curl_numpy( - v2e: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - dualL: np.array, - dualA: np.array, - edge_orientation: np.array, -) -> np.array: - uv_curl = ( - np.sum( - (u[v2e] * nx[v2e] + v[v2e] * ny[v2e]) * dualL[v2e] * edge_orientation, - axis=1, - ) - / dualA - ) - - return uv_curl -``` - -```{code-cell} ipython3 -@gtx.field_operator -def curl( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - dualL: gtx.Field[Dims[E], float], - dualA: gtx.Field[Dims[V], float], - edge_orientation: gtx.Field[Dims[V, V2EDim], float], -) -> gtx.Field[Dims[V], float]: - uv_curl = ( - neighbor_sum( - (u(V2E) * nx(V2E) + v(V2E) * ny(V2E)) * dualL(V2E) * edge_orientation, - axis=V2EDim, - ) - / dualA - ) - - return uv_curl -``` - -```{code-cell} ipython3 -def test_curl(): - backend = None - # backend = gtfn_cpu - # backend = gtfn_gpu - - edge_domain = gtx.domain({E: n_edges}) - vertex_domain = gtx.domain({V: n_vertices}) - - u = random_field(edge_domain, allocator=backend) - v = random_field(edge_domain, allocator=backend) - nx = random_field(edge_domain, allocator=backend) - ny = random_field(edge_domain, allocator=backend) - dualL = random_field(edge_domain, allocator=backend) - dualA = random_field(vertex_domain, allocator=backend) - edge_orientation = random_sign( - gtx.domain({V: n_vertices, V2EDim: 6}), allocator=backend - ) - - divergence_ref = curl_numpy( - v2e_table, - u.asnumpy(), - v.asnumpy(), - nx.asnumpy(), - ny.asnumpy(), - dualL.asnumpy(), - dualA.asnumpy(), - edge_orientation.asnumpy(), - ) - - v2e_connectivity = gtx.NeighborTableOffsetProvider(v2e_table, V, E, 6, has_skip_values=False) - - curl_gt4py = gtx.zeros(vertex_domain, allocator=backend) - - curl( - u, v, nx, ny, dualL, dualA, edge_orientation, out = curl_gt4py, offset_provider = {V2E.value: v2e_connectivity} - ) - - assert np.allclose(curl_gt4py.asnumpy(), divergence_ref) -``` - -```{code-cell} ipython3 -test_curl() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/5_vector_laplace_exercise.md b/docs/user/next/exercises/5_vector_laplace_exercise.md deleted file mode 100644 index f15c416b28..0000000000 --- a/docs/user/next/exercises/5_vector_laplace_exercise.md +++ /dev/null @@ -1,264 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# 5. Vector Laplacian - -+++ - -Starting off from laplacian formula by taking a scalar function (scalar field) and returning a scalar **p** that measures by how much the average value of f over a surface centered at p deviates from f(p). - -![](../images/laplacian.png "Laplacian") - -Look at vector Laplacian formula (relevant for Navier stokes viscous stress term) by taking a vector function (vector field) and returning a vector. - -![](../images/vector_laplacian.png "Vector Laplacian") - -Compute normal component of vector laplacian on finite volume meshes. - -![](../images/vector_laplacian_normal_component.png "Normal Component of Vector Laplacian") - -Can reuse divergence and curl as defined in previous exercises, however need two more directional gradients: - -![](../images/directional_gradient_n_picture.png "Gradient in n direction") - -![](../images/directional_gradient_n_formula.png "Gradient in n direction") - -![](../images/directional_gradient_tau_picture.png "Gradient in tau direction") - -![](../images/directional_gradient_tau_formula.png "Gradient in tau direction") - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx -``` - -```{code-cell} ipython3 -def divergence_numpy( - c2e: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - L: np.array, - A: np.array, - edge_orientation: np.array, -) -> np.array: - uv_div = ( - np.sum( - (u[c2e] * nx[c2e] + v[c2e] * ny[c2e]) * L[c2e] * edge_orientation, axis=1 - ) - / A - ) - return uv_div -``` - -```{code-cell} ipython3 -def curl_numpy( - v2e: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - dualL: np.array, - dualA: np.array, - edge_orientation: np.array, -) -> np.array: - uv_curl = ( - np.sum( - (u[v2e] * nx[v2e] + v[v2e] * ny[v2e]) * dualL[v2e] * edge_orientation, - axis=1, - ) - / dualA - ) - - return uv_curl -``` - -```{code-cell} ipython3 -def laplacian_numpy( - c2e: np.array, - v2e: np.array, - e2v: np.array, - e2c: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - L: np.array, - dualL: np.array, - tangent_orientation: np.array, - A: np.array, - dualA: np.array, - edge_orientation_vertex: np.array, - edge_orientation_cell: np.array, -) -> np.array: - # compute curl (on vertices) - uv_curl = curl_numpy(v2e, u, v, nx, ny, dualL, dualA, edge_orientation_vertex) - - # compute divergence (on cells) - uv_div = divergence_numpy(c2e, u, v, nx, ny, L, A, edge_orientation_cell) - - # first term of of nabla2 (gradient of curl) - grad_of_curl = (uv_curl[e2v[:, 1]] - uv_curl[e2v[:, 0]])*tangent_orientation/L - - # second term of of nabla2 (gradient of divergence) - grad_of_div = (uv_div[e2c[:, 1]] - uv_div[e2c[:, 0]])/dualL - - # finalize nabla2 (difference between the two gradients) - uv_nabla2 = grad_of_div - grad_of_curl - - return uv_nabla2 -``` - -```{code-cell} ipython3 -@gtx.field_operator -def divergence( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - L: gtx.Field[Dims[E], float], - A: gtx.Field[Dims[C], float], - edge_orientation: gtx.Field[Dims[C, C2EDim], float], -) -> gtx.Field[Dims[C], float]: - # compute divergence - uv_div = A - - return uv_div -``` - -```{code-cell} ipython3 -@gtx.field_operator -def curl( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - dualL: gtx.Field[Dims[E], float], - dualA: gtx.Field[Dims[V], float], - edge_orientation: gtx.Field[Dims[V, V2EDim], float], -) -> gtx.Field[Dims[V], float]: - # compute curl - uv_curl = dualA - - return uv_curl -``` - -```{code-cell} ipython3 -@gtx.field_operator -def laplacian_fvm( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - L: gtx.Field[Dims[E], float], - dualL: gtx.Field[Dims[E], float], - tangent_orientation: gtx.Field[Dims[E], float], - A: gtx.Field[Dims[C], float], - dualA: gtx.Field[Dims[V], float], - edge_orientation_vertex: gtx.Field[Dims[V, V2EDim], float], - edge_orientation_cell: gtx.Field[Dims[C, C2EDim], float], -) -> gtx.Field[Dims[E], float]: - # compute laplacian_fvm - uv_nabla2 = L - - return uv_nabla2 -``` - -```{code-cell} ipython3 -def test_laplacian(): - - backend = None - # backend = gtfn_cpu - # backend = gtfn_gpu - - edge_domain = gtx.domain({E: n_edges}) - vertex_domain = gtx.domain({V: n_vertices}) - cell_domain = gtx.domain({C: n_cells}) - - - u = random_field(edge_domain, allocator=backend) - v = random_field(edge_domain, allocator=backend) - nx = random_field(edge_domain, allocator=backend) - ny = random_field(edge_domain, allocator=backend) - L = random_field(edge_domain, allocator=backend) - dualL = random_field(edge_domain, allocator=backend) - tangent_orientation = random_field(edge_domain, allocator=backend) - A = random_field(cell_domain, allocator=backend) - dualA = random_field(vertex_domain, allocator=backend) - edge_orientation_vertex = random_sign( - gtx.domain({V: n_vertices, V2EDim: 6}), allocator=backend - ) - edge_orientation_cell = random_sign( - gtx.domain({C: n_cells, C2EDim: 3}), allocator=backend - ) - - laplacian_ref = laplacian_numpy( - c2e_table, - v2e_table, - e2v_table, - e2c_table, - u.asnumpy(), - v.asnumpy(), - nx.asnumpy(), - ny.asnumpy(), - L.asnumpy(), - dualL.asnumpy(), - tangent_orientation.asnumpy(), - A.asnumpy(), - dualA.asnumpy(), - edge_orientation_vertex.asnumpy(), - edge_orientation_cell.asnumpy(), - ) - - c2e_connectivity = gtx.NeighborTableOffsetProvider(c2e_table, C, E, 3, has_skip_values=False) - v2e_connectivity = gtx.NeighborTableOffsetProvider(v2e_table, V, E, 6, has_skip_values=False) - e2v_connectivity = gtx.NeighborTableOffsetProvider(e2v_table, E, V, 2, has_skip_values=False) - e2c_connectivity = gtx.NeighborTableOffsetProvider(e2c_table, E, C, 2, has_skip_values=False) - - - laplacian_gt4py = gtx.zeros(edge_domain, allocator=backend) - - laplacian_fvm( - u, - v, - nx, - ny, - L, - dualL, - tangent_orientation, - A, - dualA, - edge_orientation_vertex, - edge_orientation_cell, - out = laplacian_gt4py, - offset_provider = {C2E.value: c2e_connectivity, - V2E.value: v2e_connectivity, - E2V.value: e2v_connectivity, - E2C.value: e2c_connectivity, - }, - ) - - assert np.allclose(laplacian_gt4py.asnumpy(), laplacian_ref) -``` - -```{code-cell} ipython3 -test_laplacian() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/5_vector_laplace_exercise_solution.md b/docs/user/next/exercises/5_vector_laplace_exercise_solution.md deleted file mode 100644 index b58bb1c487..0000000000 --- a/docs/user/next/exercises/5_vector_laplace_exercise_solution.md +++ /dev/null @@ -1,286 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# 5. Vector Laplacian - -+++ - -Starting off from laplacian formula - -![](../images/laplacian.png "Laplacian") - -Look at vector Laplacian formula (relevant for Navier stokes viscous stress term) - -![](../images/vector_laplacian.png "Vector Laplacian") - -Compute normal component of vector laplacian on finite volume meshes - -![](../images/vector_laplacian_normal_component.png "Normal Component of Vector Laplacian") - -Can reuse divergence and curl as defined in previous exercises, however need two more directional gradients: - -![](../images/directional_gradient_n_picture.png "Gradient in n direction") - -![](../images/directional_gradient_n_formula.png "Gradient in n direction") - -![](../images/directional_gradient_tau_picture.png "Gradient in tau direction") - -![](../images/directional_gradient_tau_formula.png "Gradient in tau direction") - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx -``` - -```{code-cell} ipython3 -def divergence_numpy( - c2e: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - L: np.array, - A: np.array, - edge_orientation: np.array, -) -> np.array: - uv_div = ( - np.sum( - (u[c2e] * nx[c2e] + v[c2e] * ny[c2e]) * L[c2e] * edge_orientation, axis=1 - ) - / A - ) - return uv_div -``` - -```{code-cell} ipython3 -def curl_numpy( - v2e: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - dualL: np.array, - dualA: np.array, - edge_orientation: np.array, -) -> np.array: - uv_curl = ( - np.sum( - (u[v2e] * nx[v2e] + v[v2e] * ny[v2e]) * dualL[v2e] * edge_orientation, - axis=1, - ) - / dualA - ) - - return uv_curl -``` - -```{code-cell} ipython3 -def laplacian_numpy( - c2e: np.array, - v2e: np.array, - e2v: np.array, - e2c: np.array, - u: np.array, - v: np.array, - nx: np.array, - ny: np.array, - L: np.array, - dualL: np.array, - tangent_orientation: np.array, - A: np.array, - dualA: np.array, - edge_orientation_vertex: np.array, - edge_orientation_cell: np.array, -) -> np.array: - # compute curl (on vertices) - uv_curl = curl_numpy(v2e, u, v, nx, ny, dualL, dualA, edge_orientation_vertex) - - # compute divergence (on cells) - uv_div = divergence_numpy(c2e, u, v, nx, ny, L, A, edge_orientation_cell) - - # first term of of nabla2 (gradient of curl) - grad_of_curl = (uv_curl[e2v[:, 1]] - uv_curl[e2v[:, 0]])*tangent_orientation/L - - # second term of of nabla2 (gradient of divergence) - grad_of_div = (uv_div[e2c[:, 1]] - uv_div[e2c[:, 0]])/dualL - - # finalize nabla2 (difference between the two gradients) - uv_nabla2 = grad_of_div - grad_of_curl - - return uv_nabla2 -``` - -```{code-cell} ipython3 -@gtx.field_operator -def divergence( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - L: gtx.Field[Dims[E], float], - A: gtx.Field[Dims[C], float], - edge_orientation: gtx.Field[Dims[C, C2EDim], float], -) -> gtx.Field[Dims[C], float]: - uv_div = ( - neighbor_sum( - (u(C2E) * nx(C2E) + v(C2E) * ny(C2E)) * L(C2E) * edge_orientation, - axis=C2EDim, - ) - / A - ) - return uv_div -``` - -```{code-cell} ipython3 -@gtx.field_operator -def curl( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - dualL: gtx.Field[Dims[E], float], - dualA: gtx.Field[Dims[V], float], - edge_orientation: gtx.Field[Dims[V, V2EDim], float], -) -> gtx.Field[Dims[V], float]: - uv_curl = ( - neighbor_sum( - (u(V2E) * nx(V2E) + v(V2E) * ny(V2E)) * dualL(V2E) * edge_orientation, - axis=V2EDim, - ) - / dualA - ) - - return uv_curl -``` - -```{code-cell} ipython3 -@gtx.field_operator -def laplacian_fvm( - u: gtx.Field[Dims[E], float], - v: gtx.Field[Dims[E], float], - nx: gtx.Field[Dims[E], float], - ny: gtx.Field[Dims[E], float], - L: gtx.Field[Dims[E], float], - dualL: gtx.Field[Dims[E], float], - tangent_orientation: gtx.Field[Dims[E], float], - A: gtx.Field[Dims[C], float], - dualA: gtx.Field[Dims[V], float], - edge_orientation_vertex: gtx.Field[Dims[V, V2EDim], float], - edge_orientation_cell: gtx.Field[Dims[C, C2EDim], float], -) -> gtx.Field[Dims[E], float]: - - # compute curl (on vertices) - uv_curl = curl(u, v, nx, ny, dualL, dualA, edge_orientation_vertex) - - # compute divergence (on cells) - uv_div = divergence(u, v, nx, ny, L, A, edge_orientation_cell) - - # first term of of nabla2 (gradient of curl) - grad_of_curl = (uv_curl(E2V[1]) - uv_curl(E2V[0]))*tangent_orientation/L - - # second term of of nabla2 (gradient of divergence) - grad_of_div = (uv_div(E2C[1]) - uv_div(E2C[0]))/dualL - - # finalize nabla2 (difference between the two gradients) - uv_nabla2 = grad_of_div - grad_of_curl - - return uv_nabla2 -``` - -```{code-cell} ipython3 -def test_laplacian(): - - backend = None - # backend = gtfn_cpu - # backend = gtfn_gpu - - edge_domain = gtx.domain({E: n_edges}) - vertex_domain = gtx.domain({V: n_vertices}) - cell_domain = gtx.domain({C: n_cells}) - - - u = random_field(edge_domain, allocator=backend) - v = random_field(edge_domain, allocator=backend) - nx = random_field(edge_domain, allocator=backend) - ny = random_field(edge_domain, allocator=backend) - L = random_field(edge_domain, allocator=backend) - dualL = random_field(edge_domain, allocator=backend) - tangent_orientation = random_field(edge_domain, allocator=backend) - A = random_field(cell_domain, allocator=backend) - dualA = random_field(vertex_domain, allocator=backend) - edge_orientation_vertex = random_sign( - gtx.domain({V: n_vertices, V2EDim: 6}), allocator=backend - ) - edge_orientation_cell = random_sign( - gtx.domain({C: n_cells, C2EDim: 3}), allocator=backend - ) - - laplacian_ref = laplacian_numpy( - c2e_table, - v2e_table, - e2v_table, - e2c_table, - u.asnumpy(), - v.asnumpy(), - nx.asnumpy(), - ny.asnumpy(), - L.asnumpy(), - dualL.asnumpy(), - tangent_orientation.asnumpy(), - A.asnumpy(), - dualA.asnumpy(), - edge_orientation_vertex.asnumpy(), - edge_orientation_cell.asnumpy(), - ) - - c2e_connectivity = gtx.NeighborTableOffsetProvider(c2e_table, C, E, 3, has_skip_values=False) - v2e_connectivity = gtx.NeighborTableOffsetProvider(v2e_table, V, E, 6, has_skip_values=False) - e2v_connectivity = gtx.NeighborTableOffsetProvider(e2v_table, E, V, 2, has_skip_values=False) - e2c_connectivity = gtx.NeighborTableOffsetProvider(e2c_table, E, C, 2, has_skip_values=False) - - - laplacian_gt4py = gtx.zeros(edge_domain, allocator=backend) - - laplacian_fvm( - u, - v, - nx, - ny, - L, - dualL, - tangent_orientation, - A, - dualA, - edge_orientation_vertex, - edge_orientation_cell, - out = laplacian_gt4py, - offset_provider = {C2E.value: c2e_connectivity, - V2E.value: v2e_connectivity, - E2V.value: e2v_connectivity, - E2C.value: e2c_connectivity, - }, - ) - - assert np.allclose(laplacian_gt4py.asnumpy(), laplacian_ref) -``` - -```{code-cell} ipython3 -test_laplacian() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/6_where_domain.md b/docs/user/next/exercises/6_where_domain.md deleted file mode 100644 index da954864f2..0000000000 --- a/docs/user/next/exercises/6_where_domain.md +++ /dev/null @@ -1,163 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# Where, Offset, and domain - -+++ - -## Conditional: where - -+++ - -The `where` builtin works analogously to the numpy version (https://numpy.org/doc/stable/reference/generated/numpy.where.html) - -Both require the same 3 input arguments: -- mask: a field of booleans or an expression evaluating to this type -- true branch: a tuple, a field, or a scalar -- false branch: a tuple, a field, of a scalar - -+++ - -Take a simple numpy example, the `mask` here is a condition: - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx - -backend = None -# backend = gtfn_cpu -# backend = gtfn_gpu -``` - -```{code-cell} ipython3 -a_np = np.arange(10.0) -b_np = np.where(a_np < 6.0, a_np, a_np*10.0) -print("a_np array: {}".format(a_np)) -print("b_np array: {}".format(b_np)) -``` - -### **Task**: replicate this example in gt4py - -```{code-cell} ipython3 -# TODO implement the field_operator - - -@gtx.program(backend=backend) -def program_where(a: gtx.Field[Dims[K], float], b: gtx.Field[Dims[K], float]): - fieldop_where(a, out=b) -``` - -```{code-cell} ipython3 -def test_where(): - a = gtx.as_field([K], np.arange(10.0), allocator=backend) - b = gtx.as_field([K], np.zeros(shape=10), allocator=backend) - program_where(a, b, offset_provider={}) - - assert np.allclose(b_np, b.asnumpy()) -``` - -```{code-cell} ipython3 -test_where() -print("Test successful") -``` - -## Domain - -+++ - -The same operation can be performed in gt4py by including the `domain` keyowrd argument on `field_operator` call - -### **Task**: implement the same operation as above using `domain` instead of `where` - -```{code-cell} ipython3 -@gtx.field_operator -def fieldop_domain(a: gtx.Field[Dims[K], float]) -> gtx.Field[Dims[K], float]: - return a * 10.0 - - -@gtx.program(backend=backend) -def program_domain(a: gtx.Field[Dims[K], float], b: gtx.Field[Dims[K], float]): - ... # TODO write the call to fieldop_domain -``` - -```{code-cell} ipython3 -def test_domain(): - a = gtx.as_field([K], np.arange(10.0), allocator=backend) - b = gtx.as_field([K], np.arange(10.0), allocator=backend) - program_domain(a, b, offset_provider={}) - - assert np.allclose(b_np, b.asnumpy()) -``` - -```{code-cell} ipython3 -test_domain() -print("Test successful") -``` - -## where and domain - -A combination of `where` and `domain` is useful in cases when an offset is used which exceeds the field size. - -e.g. a field `a: gtx.Field[Dims[K], float]` with shape (10,) is applied an offset (`Koff`). - -+++ - -### **Task**: combine `domain` and `where` to account for extra indices - -Edit the code below such that: - 1. operations on field `a` are performed only up until the 8th index - 2. the domain is properly set accound for the offset - -+++ - -#### Python reference - -```{code-cell} ipython3 -a_np_result = np.zeros(shape=10) -for i in range(len(a_np)): - if a_np[i] < 8.0: - a_np_result[i] = a_np[i + 1] + a_np[i] - elif i < 9: - a_np_result[i] = a_np[i] -print("a_np_result array: {}".format(a_np_result)) -``` - -```{code-cell} ipython3 -@gtx.field_operator -def fieldop_domain_where(a: gtx.Field[Dims[K], float]) -> gtx.Field[Dims[K], float]: - return # TODO - -@gtx.program(backend=backend) -def program_domain_where(a: gtx.Field[Dims[K], float], b: gtx.Field[Dims[K], float]): - ... # TODO -``` - -```{code-cell} ipython3 -def test_domain_where(): - a = gtx.as_field([K], np.arange(10.0), allocator=backend) - b = gtx.as_field([K], np.zeros(shape=10), allocator=backend) - program_domain_where(a, b, offset_provider={"Koff": K}) - - assert np.allclose(a_np_result, b.asnumpy()) -``` - -```{code-cell} ipython3 -test_domain_where() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/6_where_domain_solutions.md b/docs/user/next/exercises/6_where_domain_solutions.md deleted file mode 100644 index 85577b7712..0000000000 --- a/docs/user/next/exercises/6_where_domain_solutions.md +++ /dev/null @@ -1,164 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# Where, Offset, and domain - -+++ - -## Conditional: where - -+++ - -The `where` builtin works analogously to the numpy version (https://numpy.org/doc/stable/reference/generated/numpy.where.html) - -Both require the same 3 input arguments: -- mask: a field of booleans or an expression evaluating to this type -- true branch: a tuple, a field, or a scalar -- false branch: a tuple, a field, of a scalar - -+++ - -Take a simple numpy example, the `mask` here is a condition: - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx - -backend = None -# backend = gtfn_cpu -# backend = gtfn_gpu -``` - -```{code-cell} ipython3 -a_np = np.arange(10.0) -b_np = np.where(a_np < 6.0, a_np, a_np*10.0) -print("a_np array: {}".format(a_np)) -print("b_np array: {}".format(b_np)) -``` - -### **Task**: replicate this example in gt4py - -```{code-cell} ipython3 -@gtx.field_operator -def fieldop_where(a: gtx.Field[Dims[K], float]) -> gtx.Field[Dims[K], float]: - return where(a < 6.0, a, a*10.0) - -@gtx.program(backend=backend) -def program_where(a: gtx.Field[Dims[K], float], b: gtx.Field[Dims[K], float]): - fieldop_where(a, out=b) -``` - -```{code-cell} ipython3 -def test_where(): - a = gtx.as_field([K], np.arange(10.0), allocator=backend) - b = gtx.as_field([K], np.zeros(shape=10), allocator=backend) - program_where(a, b, offset_provider={}) - - assert np.allclose(b_np, b.asnumpy()) -``` - -```{code-cell} ipython3 -test_where() -print("Test successful") -``` - -## Domain - -+++ - -The same operation can be performed in gt4py by including the `domain` keyowrd argument on `field_operator` call - -### **Task**: implement the same operation as above using `domain` instead of `where` - -```{code-cell} ipython3 -@gtx.field_operator -def fieldop_domain(a: gtx.Field[Dims[K], float]) -> gtx.Field[Dims[K], float]: - return a*10.0 - -@gtx.program(backend=backend) -def program_domain(a: gtx.Field[Dims[K], float], - b: gtx.Field[Dims[K], float]): - fieldop_domain(a, out=b, domain={K: (6, 10)}) -``` - -```{code-cell} ipython3 -def test_domain(): - a = gtx.as_field([K], np.arange(10.0), allocator=backend) - b = gtx.as_field([K], np.arange(10.0), allocator=backend) - program_domain(a, b, offset_provider={}) - - assert np.allclose(b_np, b.asnumpy()) -``` - -```{code-cell} ipython3 -test_domain() -print("Test successful") -``` - -## where and domain - -A combination of `where` and `domain` is useful in cases when an offset is used which exceeds the field size. - -e.g. a field `a: gtx.Field[Dims[K], float]` with shape (10,) is applied an offset (`Koff`). - -+++ - -### **Task**: combine `domain` and `where` to account for extra indices - -Edit the code below such that: - 1. operations on field `a` are performed only up until the 8th index - 2. the domain is properly set accound for the offset - -+++ - -#### Python reference - -```{code-cell} ipython3 -a_np_result = np.zeros(shape=10) -for i in range(len(a_np)): - if a_np[i] < 8.0: - a_np_result[i] = a_np[i + 1] + a_np[i] - elif i < 9: - a_np_result[i] = a_np[i] -print("a_np_result array: {}".format(a_np_result)) -``` - -```{code-cell} ipython3 -@gtx.field_operator -def fieldop_domain_where(a: gtx.Field[Dims[K], float]) -> gtx.Field[Dims[K], float]: - return where(a<8.0, a(Koff[1])+a, a) - -@gtx.program(backend=backend) -def program_domain_where(a: gtx.Field[Dims[K], float], b: gtx.Field[Dims[K], float]): - fieldop_domain_where(a, out=b, domain={K: (0, 9)}) -``` - -```{code-cell} ipython3 -def test_domain_where(): - a = gtx.as_field([K], np.arange(10.0), allocator=backend) - b = gtx.as_field([K], np.zeros(shape=10), allocator=backend) - program_domain_where(a, b, offset_provider={"Koff": K}) - - assert np.allclose(a_np_result, b.asnumpy()) -``` - -```{code-cell} ipython3 -test_domain_where() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/7_scan_operator.md b/docs/user/next/exercises/7_scan_operator.md deleted file mode 100644 index dc6313df5a..0000000000 --- a/docs/user/next/exercises/7_scan_operator.md +++ /dev/null @@ -1,189 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -## Scan operator - -+++ - -The unique feature of this operator is that it provides the return state of the previous iteration as its first argument (i.e., the result from the previous grid point). In other words, all the arguments of the current `return` will be available (as a tuple) in the next iteration from the first argument of the defined function. - -Example: A FORTRAN pseudocode for integrating a moisture variable (e.g., cloud water or water vapour) over a column could look as follows: - - -```FORTRAN -SUBROUTINE column_integral( var_in, rho, dz, var_out, ie, je, ke ) - ! Return the column integral of a moist species. - INTEGER, INTENT (IN) :: & - ie, je, ke ! array dimensions of the I/O-fields (horizontal, horizontal, vertical) - - REAL (KIND=wp), INTENT (OUT) :: & - q_colsum (ie,je) ! Vertically-integrated mass of water species - - REAL (KIND=wp), INTENT (IN) :: & - rho (ie,je,ke), & - dz (ie,je,ke), & ! height of model half levels - var_in (ie,je,ke) ! humidity mass concentration at time-level nnow - - !$acc parallel present( iq ) if (lzacc) - !$acc loop gang - DO j=1,je - !$acc loop vector - DO i=1,ie - q_sum(i,j) = 0.0 - END DO - END DO - !$acc end parallel - - - !$acc parallel present( iq, rho, hhl, q ) if (lzacc) - DO k = 1, ke ! Vertical loop - !$acc loop gang - DO j=1,je - !$acc loop vector - DO i=1,ie - q_colsum(i,j) = q_colsum(i,j) + var_in(i,j,k) * rho(i,j,k)* dz(i,j,k) - END DO - END DO - END DO - !$acc end parallel -END SUBROUTINE column_integral -``` - -Where: -- `var_in` is the 3D variable that will be summed up -- `q_colsum` is the resulting 2D variable -- `rho` the air density -- `dz`the thickness of the vertical layers - -In the first loop nest, `column_sum` is set to zero for all grid columns. The vertical dependency enters on the RHS of the second loop nest `q_colsum(i,j) = q_colsum(i,j) + ...` - -Using the `scan_operator` this operation would be written like this: - -```python -@scan_operator(axis=KDim, forward=True, init=0.0) -def column_integral(float: state, float: var, float: rho, float: dz) - """Return the column integral of a moist species.""" - return var * rho * dz + state -``` - -Here the vertical dependency is expressed by the first function argument (`state`). This argument carries the return from the previous k-level and does not need to be specified when the function is called (similar to the `self` argument of Python classes). The argument is intialized to `init=0.0` in the function decorator (first loop nest above) and the dimension of the integral is specified with `axis=KDim`. - - -```python -q_colsum = column_integral(qv, rho, dz) -``` - -+++ - -#### Exercise: port a toy cloud microphysics scheme from python/numpy using the template of a `scan_operator` below - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx - -backend = None -# backend = gtfn_cpu -# backend = gtfn_gpu -``` - -```{code-cell} ipython3 -def toy_microphysics_numpy(qc, qr, autoconversion_rate=0.1, sedimentaion_constant=0.05): - """A toy model of a microphysics scheme contaning autoconversion and scavenging""" - - sedimentation_flux = 0.0 - - for cell, k in np.ndindex(qc.shape): - # Autoconversion: Cloud Drops -> Rain Drops - autoconversion_tendency = qc[cell, k] * autoconversion_rate - - qc[cell, k] -= autoconversion_tendency - qr[cell, k] += autoconversion_tendency - - ## Apply sedimentation flux from level above - qr[cell, k] += sedimentation_flux - - ## Remove mass due to sedimentation flux from the current cell - qr[cell, k] -= sedimentation_flux -``` - -```{code-cell} ipython3 -@gtx.scan_operator(axis=K, forward=True, init=(0.0, 0.0, 0.0)) -def _graupel_toy_scan( - state: tuple[float, float, float], qc_in: float, qr_in: float -) -> tuple[float, float, float]: - autoconversion_rate = 0.1 - sedimentaion_constant = 0.05 - - # unpack state of previous iteration - # TODO - - # Autoconversion: Cloud Drops -> Rain Drops - # TODO - - ## Add sedimentation flux from level above - # TODO - - # Remove mass due to sedimentation flux - # TODO - - return # TODO -``` - -```{code-cell} ipython3 -@gtx.field_operator(backend=backend) -def graupel_toy_scan( - qc: gtx.Field[Dims[C, K], float], qr: gtx.Field[Dims[C, K], float] -) -> tuple[ - gtx.Field[Dims[C, K], float], - gtx.Field[Dims[C, K], float] -]: - qc, qr, _ = _graupel_toy_scan(qc, qr) - - return qc, qr -``` - -```{code-cell} ipython3 -def test_scan_operator(): - cell_k_domain = gtx.domain({C: n_cells, K: n_levels}) - - qc = random_field(cell_k_domain, allocator=backend) - qr = random_field(cell_k_domain, allocator=backend) - - qc_new = gtx.zeros(cell_k_domain, allocator=backend) - qr_new = gtx.zeros(cell_k_domain, allocator=backend) - - # Initialize Numpy fields from GT4Py fields - qc_numpy = qc.asnumpy().copy() - qr_numpy = qr.asnumpy().copy() - - # Execute the Numpy version of scheme - toy_microphysics_numpy(qc_numpy, qr_numpy) - - # Execute the GT4Py version of scheme - graupel_toy_scan(qc, qr, out=(qc_new, qr_new), offset_provider={}) - - # Compare results - assert np.allclose(qc_new.asnumpy(), qc_numpy) - assert np.allclose(qr_new.asnumpy(), qr_numpy) -``` - -```{code-cell} ipython3 -test_scan_operator() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/7_scan_operator_solutions.md b/docs/user/next/exercises/7_scan_operator_solutions.md deleted file mode 100644 index 6477692485..0000000000 --- a/docs/user/next/exercises/7_scan_operator_solutions.md +++ /dev/null @@ -1,191 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -## Scan operator - -+++ - -The unique feature of this operator is that it provides the return state of the previous iteration as its first argument (i.e., the result from the previous grid point). In other words, all the arguments of the current `return` will be available (as a tuple) in the next iteration from the first argument of the defined function. - -Example: A FORTRAN pseudocode for integrating a moisture variable (e.g., cloud water or water vapour) over a column could look as follows: - - -```FORTRAN -SUBROUTINE column_integral( var_in, rho, dz, var_out, ie, je, ke ) - ! Return the column integral of a moist species. - INTEGER, INTENT (IN) :: & - ie, je, ke ! array dimensions of the I/O-fields (horizontal, horizontal, vertical) - - REAL (KIND=wp), INTENT (OUT) :: & - q_colsum (ie,je) ! Vertically-integrated mass of water species - - REAL (KIND=wp), INTENT (IN) :: & - rho (ie,je,ke), & - dz (ie,je,ke), & ! height of model half levels - var_in (ie,je,ke) ! humidity mass concentration at time-level nnow - - !$acc parallel present( iq ) if (lzacc) - !$acc loop gang - DO j=1,je - !$acc loop vector - DO i=1,ie - q_sum(i,j) = 0.0 - END DO - END DO - !$acc end parallel - - - !$acc parallel present( iq, rho, hhl, q ) if (lzacc) - DO k = 1, ke ! Vertical loop - !$acc loop gang - DO j=1,je - !$acc loop vector - DO i=1,ie - q_colsum(i,j) = q_colsum(i,j) + var_in(i,j,k) * rho(i,j,k)* dz(i,j,k) - END DO - END DO - END DO - !$acc end parallel -END SUBROUTINE column_integral -``` - -Where: -- `var_in` is the 3D variable that will be summed up -- `q_colsum` is the resulting 2D variable -- `rho` the air density -- `dz`the thickness of the vertical layers - -In the first loop nest, `column_sum` is set to zero for all grid columns. The vertical dependency enters on the RHS of the second loop nest `q_colsum(i,j) = q_colsum(i,j) + ...` - -Using the `scan_operator` this operation would be written like this: - -```python -@scan_operator(axis=KDim, forward=True, init=0.0) -def column_integral(float: state, float: var, float: rho, float: dz) - """Return the column integral of a moist species.""" - return var * rho * dz + state -``` - -Here the vertical dependency is expressed by the first function argument (`state`). This argument carries the return from the previous k-level and does not need to be specified when the function is called (similar to the `self` argument of Python classes). The argument is intialized to `init=0.0` in the function decorator (first loop nest above) and the dimension of the integral is specified with `axis=KDim`. - - -```python -q_colsum = column_integral(qv, rho, dz) -``` - -+++ - -#### Exercise: port a toy cloud microphysics scheme from python/numpy using the template of a `scan_operator` below - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx - -backend = None -# backend = gtfn_cpu -# backend = gtfn_gpu -``` - -```{code-cell} ipython3 -def toy_microphysics_numpy(qc, qr, autoconversion_rate=0.1, sedimentaion_constant=0.05): - """A toy model of a microphysics scheme contaning autoconversion and scavenging""" - - sedimentation_flux = 0.0 - - for cell, k in np.ndindex(qc.shape): - # Autoconversion: Cloud Drops -> Rain Drops - autoconversion_tendency = qc[cell, k] * autoconversion_rate - - qc[cell, k] -= autoconversion_tendency - qr[cell, k] += autoconversion_tendency - - ## Apply sedimentation flux from level above - qr[cell, k] += sedimentation_flux - - ## Remove mass due to sedimentation flux from the current cell - qr[cell, k] -= sedimentation_flux -``` - -```{code-cell} ipython3 -@gtx.scan_operator(axis=K, forward=True, init=(0.0, 0.0, 0.0)) -def _graupel_toy_scan( - state: tuple[float, float, float], qc_in: float, qr_in: float -) -> tuple[float, float, float]: - autoconversion_rate = 0.1 - sedimentaion_constant = 0.05 - - # unpack state of previous iteration - _, _, sedimentation_flux = state - - # Autoconversion: Cloud Drops -> Rain Drops - autoconv_t = qc_in * autoconversion_rate - qc = qc_in - autoconv_t - qr = qr_in + autoconv_t - - ## Add sedimentation flux from level above - qr = qr + sedimentation_flux - - # Remove mass due to sedimentation flux - qr = qr - sedimentation_flux - - return qc, qr, sedimentation_flux -``` - -```{code-cell} ipython3 -@gtx.field_operator(backend=backend) -def graupel_toy_scan( - qc: gtx.Field[Dims[C, K], float], qr: gtx.Field[Dims[C, K], float] -) -> tuple[ - gtx.Field[Dims[C, K], float], - gtx.Field[Dims[C, K], float] -]: - qc, qr, _ = _graupel_toy_scan(qc, qr) - - return qc, qr -``` - -```{code-cell} ipython3 -def test_scan_operator(): - cell_k_domain = gtx.domain({C: n_cells, K: n_levels}) - - qc = random_field(cell_k_domain, allocator=backend) - qr = random_field(cell_k_domain, allocator=backend) - - qc_new = gtx.zeros(cell_k_domain, allocator=backend) - qr_new = gtx.zeros(cell_k_domain, allocator=backend) - - # Initialize Numpy fields from GT4Py fields - qc_numpy = qc.asnumpy().copy() - qr_numpy = qr.asnumpy().copy() - - # Execute the Numpy version of scheme - toy_microphysics_numpy(qc_numpy, qr_numpy) - - # Execute the GT4Py version of scheme - graupel_toy_scan(qc, qr, out=(qc_new, qr_new), offset_provider={}) - - # Compare results - assert np.allclose(qc_new.asnumpy(), qc_numpy) - assert np.allclose(qr_new.asnumpy(), qr_numpy) -``` - -```{code-cell} ipython3 -test_scan_operator() -print("Test successful") -``` - -```{code-cell} ipython3 - -``` diff --git a/docs/user/next/exercises/8_diffusion_exercise_solution.md b/docs/user/next/exercises/8_diffusion_exercise_solution.md deleted file mode 100644 index 664a68288b..0000000000 --- a/docs/user/next/exercises/8_diffusion_exercise_solution.md +++ /dev/null @@ -1,159 +0,0 @@ ---- -jupytext: - formats: ipynb,md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.1 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -```{code-cell} ipython3 -from helpers import * - -import gt4py.next as gtx -``` - -```{code-cell} ipython3 -def diffusion_step_numpy( - e2c2v: np.array, - v2e: np.array, - TE: np.array, - TE_t: np.array, - inv_primal_edge_length: np.array, - inv_vert_vert_length: np.array, - nnbhV: np.array, - boundary_edge: np.array, - kappa: float, - dt: float, -) -> np.array: - - # initialize - TEinit = TE - - # predict - TE = TEinit + 0.5*dt*TE_t - - # interpolate temperature from edges to vertices - TV = neighbor_sum(TE(v2e), axis=1) / nnbhV - - # compute nabla2 using the finite differences - TEnabla2 = neighbor_sum( - (TV(e2c2v[0]) + TV(e2c2v[1])) * inv_primal_edge_length ** 2 - + (TV(e2c2v[3]) + TV(e2c2v[4])) * inv_vert_vert_length ** 2, - axis=1, - ) - - TEnabla2 = TEnabla2 - ( - (2.0 * TE * inv_primal_edge_length ** 2) - + (2.0 * TE * inv_vert_vert_length ** 2) - ) - - # build ODEs - TE_t = where( - boundary_edge, - 0., - kappa*TEnabla2, - ) - - # correct - TE = TEinit + dt*TE_t -``` - -```{code-cell} ipython3 -@gtx.field_operator -def diffusion_step( - TE: gtx.Field[Dims[E], float], - TE_t: gtx.Field[Dims[E], float], - inv_primal_edge_length: gtx.Field[Dims[E], float], - inv_vert_vert_length: gtx.Field[Dims[E], float], - nnbhV: gtx.Field[Dims[V], float], - boundary_edge: gtx.Field[Dims[E], bool], - kappa: float, - dt: float, -) -> gtx.tuple[ - gtx.Field[Dims[E], float], - gtx.Field[Dims[E], float], -]: - - # initialize - TEinit = TE - - # predict - TE = TEinit + 0.5*dt*TE_t - - # interpolate temperature from edges to vertices - TV = neighbor_sum(TE(V2E), axis=V2EDim) / nnbhV - - # compute nabla2 using the finite differences - TEnabla2 = neighbor_sum( - (TV(E2C2V[0]) + TV(E2C2V[1])) * inv_primal_edge_length ** 2 - + (TV(E2C2V[3]) + TV(E2C2V[4])) * inv_vert_vert_length ** 2, - axis=E2C2VDim, - ) - - TEnabla2 = TEnabla2 - ( - (2.0 * TE * inv_primal_edge_length ** 2) - + (2.0 * TE * inv_vert_vert_length ** 2) - ) - - # build ODEs - TE_t = where( - boundary_edge, - 0., - kappa*TEnabla2, - ) - - # correct - TE = TEinit + dt*TE_t - - return TE_t, TE -``` - -```{code-cell} ipython3 -def test_diffusion_step(): - u = random_field((n_edges), E) - v = random_field((n_edges), E) - nx = random_field((n_edges), E) - ny = random_field((n_edges), E) - L = random_field((n_edges), E) - dualL = random_field((n_edges), E) - A = random_field((n_cells), C) - dualA = random_field((n_vertices), V) - edge_orientation_vertex = random_field((n_cells, 6), V, V2EDim) - edge_orientation_cell = random_field((n_cells, 3), C, C2EDim) - - laplacian_ref = laplacian_numpy( - c2e_table, - u.asnumpy(), - v.asnumpy(), - nx.asnumpy(), - ny.asnumpy(), - L.asnumpy(), - A.asnumpy(), - edge_orientation.asnumpy(), - ) - - c2e_connectivity = gtx.NeighborTableOffsetProvider(c2e_table, C, E, 3) - - laplacian_gt4py = zero_field((n_edges), E) - - laplacian_fvm( - u, v, nx, ny, L, A, edge_orientation, out = divergence_gt4py, offset_provider = {C2E.value: c2e_connectivity} - ) - - assert np.allclose(divergence_gt4py.asnumpy(), divergence_ref) -``` - -```{code-cell} ipython3 -test_diffusion_step() -print("Test successful") -``` - -```{code-cell} ipython3 - -```