Skip to content

Added solution plot function #16

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

Closed
wants to merge 5 commits into from
Closed

Added solution plot function #16

Changes from all commits
Commits
Show all changes
5 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
b7ef4ac
Merge branch 'master' into master
yguclu Mar 12, 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
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",
Copy link

Choose a reason for hiding this comment

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

git sees some changes because this import now comes after the line break. If you keep the same line breaks than in the original code, git will only see the relevant changes and the commit will be a bit cleaner.

"\n",
"V = ScalarFunctionSpace('V', domain)\n",
"domain = Square()\n",
Copy link

Choose a reason for hiding this comment

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

same here: the only changes here are the position of the line breaks "\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",
"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",
Copy link

Choose a reason for hiding this comment

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

it seems to me that this line triggers an error when the notebook is executed. @anushkasinghh did you verify that you could run the notebooks ?

When the CI tests are run with the PR (probably after each new commit) you can also check that the "book" (which will be deployed after this PR is merged) looks as it should, by downloading the artifacts as described in the readme.

But you should first check that you can run the notebook correctly before committing, that will be simpler :)

" 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"
}
},
Comment on lines +375 to +393
Copy link

Choose a reason for hiding this comment

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

It seems to me that this metadata can/should be cleaned before a commit. @yguclu do you agree?

Copy link
Member

Choose a reason for hiding this comment

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

I am no Jupyter Notebook expert, but I think so, yes.

Copy link
Member

Choose a reason for hiding this comment

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

Not Yaman but I agree @campospinto. The metadata here is hardcoding Python version info which we're trying to avoid. @anushkasinghh please clean the notebook as stated in the guide:

# Clean notebook
nb-clean clean --remove-empty-cells --remove-all-notebook-metadata chapter1/poisson.ipynb

# Check if notebook is clean. No output means the notebook is already clean (otherwise you'll get an error)
nb-clean check chapter1/poisson.ipynb

Copy link
Member

Choose a reason for hiding this comment

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

@kvrigor I have noticed the cell IDs being changed in this PR. Is this something unavoidable, or we could keep the original values somehow?

Copy link
Member

Choose a reason for hiding this comment

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

My surface reading of this doc suggests that cell IDs should change whenever cells change, since the primary purpose of a cell ID is to uniquely identify cell contents. Unfortunately we can't avoid mutating cell IDs, and maybe that's one of the reasons why a Jupyter notebook isn't version-control friendly.

Copy link
Member

Choose a reason for hiding this comment

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

Is a cell ID something like a hash code? It seems to me that a new one could be generated when opening the notebook, rather than hardcoded in the file...

Copy link
Member

Choose a reason for hiding this comment

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

Is a cell ID something like a hash code?

It's possible; hardcoding it might be necessary for efficiency and also for scripting access by IDE tools.

"nbformat": 4,
"nbformat_minor": 5
}
Loading