Skip to content

Add plots, remove extra changes #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
anushkasinghh wants to merge 9 commits into
base: master
Choose a base branch
from
Open

Add plots, remove extra changes #20

Changes from 1 commit
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
42cff1a
Add solution plots to poisson solver notebook
anushkasinghh Mar 10, 2025
92be763
Add solution plotting function
anushkasinghh Mar 11, 2025
269fc63
Merge pull request #1 from anushkasinghh/anushka
anushkasinghh Mar 11, 2025
b5792df
Reverting incorrect changes
Mar 11, 2025
169331d
cleaned nb
Mar 12, 2025
4ccca1f
Merge branch 'pyccel:master' into anushka
anushkasinghh Mar 13, 2025
2698e29
nan
Mar 13, 2025
464d966
nan
Mar 13, 2025
958ede9
Merge branch 'master' into add-plots
yguclu Jul 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Add solution plotting function
  • Loading branch information
@anushkasinghh
anushkasinghh authored Mar 11, 2025
commit 92be763111685bf40ef6804fe6e98a157caf6e7b
127 changes: 121 additions & 6 deletions chapter1/poisson.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,12 +55,16 @@
"from sympde.topology import ScalarFunctionSpace, Square, element_of\n",
"from sympde.calculus import grad, dot\n",
"\n",
"from sympy import pi, sin\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"\n",
"domain = Square()\n",
"from sympy import pi, sin, lambdify\n",
"\n",
"V = ScalarFunctionSpace('V', domain)\n",
"domain = Square()\n",
"\n",
"V = ScalarFunctionSpace('V', domain) \n",
" \n",
"x,y = domain.coordinates\n",
"\n",
"u,v = [element_of(V, name=i) for i in ['u', 'v']]\n",
Expand Down Expand Up @@ -202,11 +206,11 @@
"outputs": [],
"source": [
"ue = sin(pi*x)*sin(pi*y)\n",
"\n",
"u = element_of(V, name='u')\n",
"\n",
"# create the formal Norm object\n",
"l2norm = Norm(u - ue, domain, kind='l2')\n",
"print(l2norm)\n",
"\n",
"# discretize the norm\n",
"l2norm_h = discretize(l2norm, domain_h, Vh)\n",
Expand Down Expand Up @@ -257,7 +261,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "d829e410",
"id": "827c3e69-77ac-4312-a4dd-a1c26a40b27c",

@yguclu yguclu Mar 17, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that this id can be changed by hand to its original value using a text editor (or from the GitHub web interface). If it causes a problem, our CI build will spot it.

"metadata": {},
"outputs": [],
"source": [
Expand All @@ -273,9 +277,120 @@
"# print the result\n",
"print(h1_error)"
]
},
{
"cell_type": "markdown",
"id": "14d92a71-8c0f-4a91-8f56-c405f3fc76d1",
"metadata": {},
"source": [
"### Visualization\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "968a2165-ab75-4bf0-b943-daab3bd9fa2d",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"\n",
"def refine_array_1d(breaks, N):\n",
" \"\"\"Refine a 1D array by inserting N points between each pair of original values.\"\"\"\n",
" refined = np.concatenate([np.linspace(breaks[i], breaks[i+1], N+1)[:-1] for i in range(len(breaks)-1)])\n",
" return np.append(refined, breaks[-1]) # Add the last point\n",
"\n",
"def plot_solutions(Vh, uh, ue, N=10):\n",
" \"\"\"\n",
" Refine the grid, evaluate solutions, compute errors, and plot results.\n",
"\n",
" Parameters:\n",
" Vh: The finite element space (must have `spaces` attribute with `breaks`).\n",
" uh: The numerical solution function.\n",
" ue: The exact solution function.\n",
" N: Number of points to insert between breaks (default: 10).\n",
" \"\"\"\n",
" # Generate refined grid for visualization\n",
" eta1 = refine_array_1d(Vh.spaces[0].breaks, N)\n",
" eta2 = refine_array_1d(Vh.spaces[1].breaks, N)\n",
"\n",
" # Evaluate numerical solution on the refined grid\n",
" num = np.array([[uh(e1, e2) for e2 in eta2] for e1 in eta1])\n",
"\n",
" # Compute exact solution\n",
" ex = np.array([[phi_exact(e1, e2) for e2 in eta2] for e1 in eta1])\n",
" err = num - ex\n",
"\n",
" # Generate physical grid coordinates\n",
" xx, yy = np.meshgrid(eta1, eta2, indexing='ij')\n",
"\n",
" # Create figure with 3 subplots\n",
" fig, axes = plt.subplots(1, 3, figsize=(12.8, 4.8))\n",
"\n",
" def add_colorbar(im, ax):\n",
" \"\"\"Add a colorbar to the given axis.\"\"\"\n",
" divider = make_axes_locatable(ax)\n",
" cax = divider.append_axes(\"right\", size=0.2, pad=0.2)\n",
" cbar = ax.get_figure().colorbar(im, cax=cax)\n",
" return cbar\n",
"\n",
" # Plot exact solution\n",
" ax = axes[0]\n",
" im = ax.contourf(xx, yy, ex, 40, cmap='jet')\n",
" add_colorbar(im, ax)\n",
" ax.set_xlabel(r'$x$', rotation='horizontal')\n",
" ax.set_ylabel(r'$y$', rotation='horizontal')\n",
" ax.set_title(r'$\\phi_{exact}(x,y)$')\n",
" ax.set_aspect('equal')\n",
"\n",
" # Plot numerical solution\n",
" ax = axes[1]\n",
" im = ax.contourf(xx, yy, num, 40, cmap='jet')\n",
" add_colorbar(im, ax)\n",
" ax.set_xlabel(r'$x$', rotation='horizontal')\n",
" ax.set_ylabel(r'$y$', rotation='horizontal')\n",
" ax.set_title(r'$\\phi(x,y)$')\n",
" ax.set_aspect('equal')\n",
"\n",
" # Plot numerical error\n",
" ax = axes[2]\n",
" im = ax.contourf(xx, yy, err, 40, cmap='jet')\n",
" add_colorbar(im, ax)\n",
" ax.set_xlabel(r'$x$', rotation='horizontal')\n",
" ax.set_ylabel(r'$y$', rotation='horizontal')\n",
" ax.set_title(r'$\\phi(x,y) - \\phi_{exact}(x,y)$')\n",
" ax.set_aspect('equal')\n",
"\n",
" # Show figure\n",
" plt.tight_layout()\n",
" fig.show()\n",
"\n",
"ue = lambdify((x, y), u, 'numpy') # convert sympy function to a callable function \n",
"plot_solutions(Vh, uh, ue, N=10)"
]
}
],
"metadata": {},
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}