|
22 | 22 | }, |
23 | 23 | "outputs": [], |
24 | 24 | "source": [ |
25 | | - "import torch\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom neuralop.layers.spectral_projection import spectral_projection_divergence_free\nfrom neuralop.losses.differentiation import FourierDiff, FiniteDiff\n \ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(f\"Using device: {device}\")" |
| 25 | + "import torch\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom neuralop.layers.spectral_projection import spectral_projection_divergence_free\nfrom neuralop.losses.differentiation import FourierDiff, FiniteDiff\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(f\"Using device: {device}\")" |
26 | 26 | ] |
27 | 27 | }, |
28 | 28 | { |
29 | 29 | "cell_type": "markdown", |
30 | 30 | "metadata": {}, |
31 | 31 | "source": [ |
32 | | - ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Divergence error computation functions\nWe define two functions to compute the divergence error \nusing spectral differentiation and finite differences.\n\n" |
| 32 | + ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Divergence error computation functions\nWe define two functions to compute the divergence error\nusing spectral differentiation and finite differences.\n\n" |
33 | 33 | ] |
34 | 34 | }, |
35 | 35 | { |
|
47 | 47 | "cell_type": "markdown", |
48 | 48 | "metadata": {}, |
49 | 49 | "source": [ |
50 | | - ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Setting considered\nWe start from a divergence-free velocity field on [0, 2*pi] x [0, 2*pi] \nconstructed from the stream function \u03c8(x,y) = sin(x) * cos(6*y) \n\nVelocity components: \n u_x = \u2202\u03c8/\u2202y = -6 * sin(x) * sin(6*y) \n\n u_y = -\u2202\u03c8/\u2202x = - cos(x) * cos(6*y) \n\n-------------------------------------------------------\n\nMathematical verification of divergence-free property: \n\n \u2207\u00b7u = \u2202u_x/\u2202x + \u2202u_y/\u2202y \n\n = \u2202/\u2202x[-6*sin(x) * sin(6*y)] + \u2202/\u2202y[-cos(x) * cos(6*y)] \n\n = -6*cos(x) * sin(6*y) + 6*cos(x) * sin(6*y) \n\n = 0 \u2713 \n\n-------------------------------------------------------\n\n\nWe then add 10% noise to break the divergence-free property. \n\nWe then apply the spectral projection to restore the divergence-free property. \n\nWe then compute the divergence error for the original, noisy, and projected fields. \n\nWe repeat this at various resolutions, [256, 512, 1024, 2048, 4096, 8192].\n\n" |
| 50 | + ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Setting considered\nWe start from a divergence-free velocity field on [0, 2*pi] x [0, 2*pi]\nconstructed from the stream function \u03c8(x,y) = sin(x) * cos(6*y)\n\nVelocity components:\n u_x = \u2202\u03c8/\u2202y = -6 * sin(x) * sin(6*y)\n\n u_y = -\u2202\u03c8/\u2202x = - cos(x) * cos(6*y)\n\n-------------------------------------------------------\n\nMathematical verification of divergence-free property:\n\n \u2207\u00b7u = \u2202u_x/\u2202x + \u2202u_y/\u2202y\n\n = \u2202/\u2202x[-6*sin(x) * sin(6*y)] + \u2202/\u2202y[-cos(x) * cos(6*y)]\n\n = -6*cos(x) * sin(6*y) + 6*cos(x) * sin(6*y)\n\n = 0 \u2713\n\n-------------------------------------------------------\n\n\nWe then add 10% noise to break the divergence-free property.\n\nWe then apply the spectral projection to restore the divergence-free property.\n\nWe then compute the divergence error for the original, noisy, and projected fields.\n\nWe repeat this at various resolutions, [256, 512, 1024, 2048, 4096, 8192].\n\n" |
51 | 51 | ] |
52 | 52 | }, |
53 | 53 | { |
|
58 | 58 | }, |
59 | 59 | "outputs": [], |
60 | 60 | "source": [ |
61 | | - "L = 2 * np.pi\nnoise_level = 0.1\nresolutions = [256, 512, 1024, 2048, 4096, 8192]\n\n\nerrors_original_spectral = []\nerrors_original_finite = []\nerrors_noisy_spectral = []\nerrors_noisy_finite = []\nerrors_prog_spectral = []\nerrors_prog_finite = []\n\nfor target_resolution in resolutions:\n \n # Create coordinate grids for this resolution\n xs = torch.arange(target_resolution, device=device, dtype=torch.float64) * (L / target_resolution)\n ys = torch.arange(target_resolution, device=device, dtype=torch.float64) * (L / target_resolution)\n X, Y = torch.meshgrid(xs, ys, indexing='ij')\n \n # Create divergence-free field using the stream function defined earlier\n u_x = -6.0 * torch.sin(X) * torch.sin(6.0 * Y)\n u_y = -torch.cos(X) * torch.cos(6.0 * Y)\n u = torch.stack([u_x, u_y], dim=0).unsqueeze(0).to(device=device, dtype=torch.float64)\n \n # Add noise to break divergence-free property\n mean_magnitude = torch.mean(torch.sqrt(u[:, 0]**2 + u[:, 1]**2))\n noise = torch.randn_like(u, dtype=torch.float64) * noise_level * mean_magnitude\n u_noisy = u + noise\n \n # Apply spectral projection to restore divergence-free property\n u_proj = spectral_projection_divergence_free(u_noisy, L, constraint_modes=(64, 64))\n \n # Compute divergence errors for all three fields\n errors_original_spectral.append(div_error_fourier(u, L))\n errors_original_finite.append(div_error_finite_diff(u, L))\n errors_noisy_spectral.append(div_error_fourier(u_noisy, L))\n errors_noisy_finite.append(div_error_finite_diff(u_noisy, L))\n errors_prog_spectral.append(div_error_fourier(u_proj, L))\n errors_prog_finite.append(div_error_finite_diff(u_proj, L))" |
| 61 | + "L = 2 * np.pi\nnoise_level = 0.1\nresolutions = [256, 512, 1024, 2048, 4096, 8192]\n\n\nerrors_original_spectral = []\nerrors_original_finite = []\nerrors_noisy_spectral = []\nerrors_noisy_finite = []\nerrors_prog_spectral = []\nerrors_prog_finite = []\n\nfor target_resolution in resolutions:\n # Create coordinate grids for this resolution\n xs = torch.arange(target_resolution, device=device, dtype=torch.float64) * (L / target_resolution)\n ys = torch.arange(target_resolution, device=device, dtype=torch.float64) * (L / target_resolution)\n X, Y = torch.meshgrid(xs, ys, indexing=\"ij\")\n \n # Create divergence-free field using the stream function defined earlier\n u_x = -6.0 * torch.sin(X) * torch.sin(6.0 * Y)\n u_y = -torch.cos(X) * torch.cos(6.0 * Y)\n u = torch.stack([u_x, u_y], dim=0).unsqueeze(0).to(device=device, dtype=torch.float64)\n \n # Add noise to break divergence-free property\n mean_magnitude = torch.mean(torch.sqrt(u[:, 0] ** 2 + u[:, 1] ** 2))\n noise = torch.randn_like(u, dtype=torch.float64) * noise_level * mean_magnitude\n u_noisy = u + noise\n\n # Apply spectral projection to restore divergence-free property\n u_proj = spectral_projection_divergence_free(u_noisy, L, constraint_modes=(64, 64))\n\n # Compute divergence errors for all three fields\n errors_original_spectral.append(div_error_fourier(u, L))\n errors_original_finite.append(div_error_finite_diff(u, L))\n errors_noisy_spectral.append(div_error_fourier(u_noisy, L))\n errors_noisy_finite.append(div_error_finite_diff(u_noisy, L))\n errors_prog_spectral.append(div_error_fourier(u_proj, L))\n errors_prog_finite.append(div_error_finite_diff(u_proj, L))" |
62 | 62 | ] |
63 | 63 | }, |
64 | 64 | { |
65 | 65 | "cell_type": "markdown", |
66 | 66 | "metadata": {}, |
67 | 67 | "source": [ |
68 | | - ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Divergence Errors using Spectral Differentiation \nThe Fourier differentiation method computes derivatives in the spectral domain \nby transforming the field to Fourier space, applying the appropriate wavenumber \noperators, and transforming back. \n\nWe display the divergence error for the original, noisy, and projected fields \nat the different resolutions. \n\nNote that at lower resolutions, finite differences are not accurate enough \nto properly compute the divergence error, which is why the errors appear higher \ninitially but improve at higher resolutions. This is a limitation of \nfinite differences for computing derivatives, not an issue with the \nspectral projection itself. Spectral differentiation provides more accurate \nderivative calculations at lower resolutions.\n\n" |
| 68 | + ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Divergence Errors using Spectral Differentiation\nThe Fourier differentiation method computes derivatives in the spectral domain\nby transforming the field to Fourier space, applying the appropriate wavenumber\noperators, and transforming back.\n\nWe display the divergence error for the original, noisy, and projected fields\nat the different resolutions.\n\nNote that at lower resolutions, finite differences are not accurate enough\nto properly compute the divergence error, which is why the errors appear higher\ninitially but improve at higher resolutions. This is a limitation of\nfinite differences for computing derivatives, not an issue with the\nspectral projection itself. Spectral differentiation provides more accurate\nderivative calculations at lower resolutions.\n\n" |
69 | 69 | ] |
70 | 70 | }, |
71 | 71 | { |
|
76 | 76 | }, |
77 | 77 | "outputs": [], |
78 | 78 | "source": [ |
79 | | - "# Spectral Differentiation table\nprint(\"-\" * 55)\nprint(f\"{'Resolution':<12} {'Original':<15} {'Noisy':<15} {'Projected':<15}\")\nprint(\"-\" * 55)\n\nfor i, res in enumerate(resolutions):\n print(f\"{res:<12} {errors_original_spectral[i]:<15.2e} {errors_noisy_spectral[i]:<15.2e} {errors_prog_spectral[i]:<15.2e}\")\n\n\n# Spectral Differentiation plot\nplt.figure(figsize=(10, 6))\nplt.semilogy(resolutions, errors_original_spectral, 'o-', label='Original', color='black', linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_noisy_spectral, 'o-', label='Noisy', color='green', linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_prog_spectral, 'o-', label='Projected', color='blue', linewidth=2.5, markersize=6)\nplt.xlabel('Resolution', fontsize=16)\nplt.ylabel('Divergence Error (L2 norm)', fontsize=16)\nplt.title('Spectral Divergence Errors', fontsize=18)\nplt.legend(fontsize=14)\nplt.grid(True, alpha=0.3)\nplt.xticks(fontsize=14)\nplt.yticks(fontsize=14)\nplt.tight_layout()\nplt.show()" |
| 79 | + "# Spectral Differentiation table\nprint(\"-\" * 55)\nprint(f\"{'Resolution':<12} {'Original':<15} {'Noisy':<15} {'Projected':<15}\")\nprint(\"-\" * 55)\n\nfor i, res in enumerate(resolutions):\n print(\n f\"{res:<12} {errors_original_spectral[i]:<15.2e} {errors_noisy_spectral[i]:<15.2e} {errors_prog_spectral[i]:<15.2e}\"\n )\n\n\n# Spectral Differentiation plot\nplt.figure(figsize=(10, 6))\nplt.semilogy(resolutions, errors_original_spectral, \"o-\", label=\"Original\", color=\"black\", linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_noisy_spectral, \"o-\", label=\"Noisy\", color=\"green\", linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_prog_spectral, \"o-\", label=\"Projected\", color=\"blue\", linewidth=2.5, markersize=6)\nplt.xlabel(\"Resolution\", fontsize=16)\nplt.ylabel(\"Divergence Error (L2 norm)\", fontsize=16)\nplt.title(\"Spectral Divergence Errors\", fontsize=18)\nplt.legend(fontsize=14)\nplt.grid(True, alpha=0.3)\nplt.xticks(fontsize=14)\nplt.yticks(fontsize=14)\nplt.tight_layout()\nplt.show()" |
80 | 80 | ] |
81 | 81 | }, |
82 | 82 | { |
83 | 83 | "cell_type": "markdown", |
84 | 84 | "metadata": {}, |
85 | 85 | "source": [ |
86 | | - ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Divergence Errors using Finite Differences\nThe finite difference method approximates derivatives using central differences. \n\nWe display the divergence error for the original, noisy, and projected fields \nat the different resolutions.\n\n" |
| 86 | + ".. raw:: html\n\n <div style=\"margin-top: 3em;\"></div>\n\n## Divergence Errors using Finite Differences\nThe finite difference method approximates derivatives using central differences.\n\nWe display the divergence error for the original, noisy, and projected fields\nat the different resolutions.\n\n" |
87 | 87 | ] |
88 | 88 | }, |
89 | 89 | { |
|
94 | 94 | }, |
95 | 95 | "outputs": [], |
96 | 96 | "source": [ |
97 | | - "# Finite differences table\nprint(\"-\" * 55)\nprint(f\"{'Resolution':<12} {'Original':<15} {'Noisy':<15} {'Projected':<15}\")\nprint(\"-\" * 55)\n\nfor i, res in enumerate(resolutions):\n print(f\"{res:<12} {errors_original_finite[i]:<15.2e} {errors_noisy_finite[i]:<15.2e} {errors_prog_finite[i]:<15.2e}\")\n\n\n# Finite differences plot\nplt.figure(figsize=(10, 6))\nplt.semilogy(resolutions, errors_original_finite, 'o-', label='Original', color='black', linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_noisy_finite, 'o-', label='Noisy', color='green', linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_prog_finite, 'o-', label='Projected', color='blue', linewidth=2.5, markersize=6)\nplt.xlabel('Resolution', fontsize=16)\nplt.ylabel('Divergence Error (L2 norm)', fontsize=16)\nplt.title('Finite Difference Divergence Errors', fontsize=18)\nplt.legend(fontsize=14)\nplt.grid(True, alpha=0.3)\nplt.xticks(fontsize=14)\nplt.yticks(fontsize=14)\nplt.tight_layout()\nplt.show()" |
| 97 | + "# Finite differences table\nprint(\"-\" * 55)\nprint(f\"{'Resolution':<12} {'Original':<15} {'Noisy':<15} {'Projected':<15}\")\nprint(\"-\" * 55)\n\nfor i, res in enumerate(resolutions):\n print(\n f\"{res:<12} {errors_original_finite[i]:<15.2e} {errors_noisy_finite[i]:<15.2e} {errors_prog_finite[i]:<15.2e}\"\n )\n\n\n# Finite differences plot\nplt.figure(figsize=(10, 6))\nplt.semilogy(resolutions, errors_original_finite, \"o-\", label=\"Original\", color=\"black\", linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_noisy_finite, \"o-\", label=\"Noisy\", color=\"green\", linewidth=2.5, markersize=6)\nplt.semilogy(resolutions, errors_prog_finite, \"o-\", label=\"Projected\", color=\"blue\", linewidth=2.5, markersize=6)\nplt.xlabel(\"Resolution\", fontsize=16)\nplt.ylabel(\"Divergence Error (L2 norm)\", fontsize=16)\nplt.title(\"Finite Difference Divergence Errors\", fontsize=18)\nplt.legend(fontsize=14)\nplt.grid(True, alpha=0.3)\nplt.xticks(fontsize=14)\nplt.yticks(fontsize=14)\nplt.tight_layout()\nplt.show()" |
98 | 98 | ] |
99 | 99 | } |
100 | 100 | ], |
|
0 commit comments