{
"cells": [
{
"cell_type": "markdown",
"id": "b40cde8e-a20e-4904-bfec-753a7d82a300",
"metadata": {},
"source": [
"# How to discretize equations \n",
"\n",
"Now we demonstrate how to discretize a PDE object.\n",
"\n",
"Pre-coded sample objects are stored in `sample` attribute of *phyem*. Invoke these samples by"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d603a38d-c76a-4cf6-b5fc-bbd27c54db4f",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:57.387182Z",
"iopub.status.busy": "2024-11-28T03:41:57.387182Z",
"iopub.status.idle": "2024-11-28T03:41:58.651258Z",
"shell.execute_reply": "2024-11-28T03:41:58.651258Z"
},
"is_executing": true,
"tags": []
},
"outputs": [],
"source": [
"import sys\n",
"ph_dir = '../../../../../' # the path to dir that containing the phyem package.\n",
"sys.path.append(ph_dir)\n",
"import phyem as ph # import the phyem package\n",
"ph.config._set_matplot_block(False)\n",
"samples = ph.samples"
]
},
{
"cell_type": "markdown",
"id": "498e0ee3-b788-4265-81a9-ce957c466895",
"metadata": {},
"source": [
"The partial differential equations of the canocical linear por-Hamiltonian are pre-coded. Call it through"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "a343b362-5677-42d7-8dc2-cdb3cca95e45",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:58.767724Z",
"iopub.status.busy": "2024-11-28T03:41:58.767724Z",
"iopub.status.idle": "2024-11-28T03:41:58.770716Z",
"shell.execute_reply": "2024-11-28T03:41:58.770716Z"
},
"is_executing": true,
"tags": []
},
"outputs": [],
"source": [
"oph = samples.pde_canonical_pH(n=3, p=3, periodic=True)[0] # where o on `oph` means outer"
]
},
{
"cell_type": "markdown",
"id": "702b2712-bd83-468e-9478-3e9c608d1a55",
"metadata": {},
"source": [
"Check `oph`, *outer oriented port-Hamiltonian*, with"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "8fd0eb4d-6350-4268-a3a3-c905926b46a4",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:58.800932Z",
"iopub.status.busy": "2024-11-28T03:41:58.800932Z",
"iopub.status.idle": "2024-11-28T03:41:59.238138Z",
"shell.execute_reply": "2024-11-28T03:41:59.238138Z"
},
"is_executing": true,
"tags": []
},
"outputs": [
{
"data": {
"text/plain": [
"<Figure size 800x600 with 1 Axes>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"oph.pr()"
]
},
{
"cell_type": "markdown",
"id": "95eaceec-a95e-4da5-82c7-7e7b1a8770a5",
"metadata": {},
"source": [
"We can take out the knowns and label them by"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "6bb4c1e9-e284-4de6-bf5a-746a707fb08f",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:59.315376Z",
"iopub.status.busy": "2024-11-28T03:41:59.314309Z",
"iopub.status.idle": "2024-11-28T03:41:59.317308Z",
"shell.execute_reply": "2024-11-28T03:41:59.317308Z"
},
"is_executing": true,
"tags": []
},
"outputs": [],
"source": [
"a3, b2 = oph.unknowns"
]
},
{
"cell_type": "markdown",
"id": "4c0f26a4-49b6-4d3c-aca3-e355f35449c9",
"metadata": {},
"source": [
"We now test `oph` with test functions from the spaces where `a3` and `b2` come from, and label the test functions by $v^3$ and $u^2$,"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "01e41ea7-fbd8-4e0a-80e0-771bd1498ded",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:59.345538Z",
"iopub.status.busy": "2024-11-28T03:41:59.345538Z",
"iopub.status.idle": "2024-11-28T03:41:59.348954Z",
"shell.execute_reply": "2024-11-28T03:41:59.348954Z"
},
"is_executing": true,
"tags": []
},
"outputs": [],
"source": [
"wf = oph.test_with(oph.unknowns, sym_repr=[r'v^3', r'u^2'])"
]
},
{
"cell_type": "markdown",
"id": "a627673d-fbaa-4b35-9f4b-7184f265d382",
"metadata": {
"tags": []
},
"source": [
"Now, we apply integration by parts to the term indexed by `'1-1'`."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "fec1529c-8a7a-4b5a-9e4c-0c4580f06729",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:59.377659Z",
"iopub.status.busy": "2024-11-28T03:41:59.377659Z",
"iopub.status.idle": "2024-11-28T03:41:59.380746Z",
"shell.execute_reply": "2024-11-28T03:41:59.380746Z"
},
"is_executing": true,
"tags": []
},
"outputs": [],
"source": [
"wf = wf.derive.integration_by_parts('1-1')"
]
},
{
"cell_type": "markdown",
"id": "65b936bb-c33f-4a14-a0cc-c67ae5d1c00f",
"metadata": {
"tags": []
},
"source": [
"We now apply a particular discretization to this weak formulation,"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "99d09aa1-ef0f-48ae-b104-bb128fd380bb",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:59.409431Z",
"iopub.status.busy": "2024-11-28T03:41:59.409431Z",
"iopub.status.idle": "2024-11-28T03:41:59.416556Z",
"shell.execute_reply": "2024-11-28T03:41:59.416556Z"
},
"is_executing": true,
"tags": []
},
"outputs": [],
"source": [
"td = wf.td\n",
"td.set_time_sequence() # initialize a time sequence\n",
"td.define_abstract_time_instants('k-1', 'k-1/2', 'k')\n",
"td.differentiate('0-0', 'k-1', 'k')\n",
"td.average('0-1', b2, ['k-1', 'k'])\n",
"\n",
"td.differentiate('1-0', 'k-1', 'k')\n",
"td.average('1-1', a3, ['k-1', 'k'])\n",
"dt = td.time_sequence.make_time_interval('k-1', 'k')\n",
"\n",
"wf = td()\n",
"\n",
"wf.unknowns = [\n",
" a3 @ td.time_sequence['k'],\n",
" b2 @ td.time_sequence['k'],\n",
"]\n",
"\n",
"wf = wf.derive.split(\n",
" '0-0', 'f0',\n",
" [a3 @ td.ts['k'], a3 @ td.ts['k-1']],\n",
" ['+', '-'],\n",
" factors=[1/dt, 1/dt],\n",
")\n",
"\n",
"wf = wf.derive.split(\n",
" '0-2', 'f0',\n",
" [ph.d(b2 @ td.ts['k-1']), ph.d(b2 @ td.ts['k'])],\n",
" ['+', '+'],\n",
" factors=[1/2, 1/2],\n",
")\n",
"\n",
"wf = wf.derive.split(\n",
" '1-0', 'f0',\n",
" [b2 @ td.ts['k'], b2 @ td.ts['k-1']],\n",
" ['+', '-'],\n",
" factors=[1/dt, 1/dt]\n",
")\n",
"\n",
"wf = wf.derive.split(\n",
" '1-2', 'f0',\n",
" [a3 @ td.ts['k-1'], a3 @ td.ts['k']],\n",
" ['+', '+'],\n",
" factors=[1/2, 1/2],\n",
")\n",
"\n",
"wf = wf.derive.rearrange(\n",
" {\n",
" 0: '0, 3 = 2, 1',\n",
" 1: '3, 0 = 2, 1',\n",
" }\n",
")"
]
},
{
"cell_type": "markdown",
"id": "a0de472f-32b1-4c27-93cc-0abe5cbab7b9",
"metadata": {},
"source": [
"We now can write the weak formulation with matrix proxies. Before doing that, we need to "
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "5e30efae-221f-4969-b285-243686f4c1f9",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:59.455344Z",
"iopub.status.busy": "2024-11-28T03:41:59.455344Z",
"iopub.status.idle": "2024-11-28T03:41:59.460313Z",
"shell.execute_reply": "2024-11-28T03:41:59.460313Z"
},
"is_executing": true,
"tags": []
},
"outputs": [],
"source": [
"ph.space.finite(3)\n",
"mp = wf.mp()"
]
},
{
"cell_type": "markdown",
"id": "2bf04ccc-e555-4eed-9eae-3dc55f1ac324",
"metadata": {},
"source": [
"The matrix format of the weak formulation leads to a linear system,"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "7bffe47c-a9bd-4bcf-95e8-f8ccf14a0ada",
"metadata": {
"execution": {
"iopub.execute_input": "2024-11-28T03:41:59.486018Z",
"iopub.status.busy": "2024-11-28T03:41:59.486018Z",
"iopub.status.idle": "2024-11-28T03:41:59.573818Z",
"shell.execute_reply": "2024-11-28T03:41:59.573818Z"
},
"is_executing": true,
"tags": []
},
"outputs": [
{
"data": {
"text/plain": [
"<Figure size 1200x600 with 1 Axes>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ls = mp.ls()\n",
"ls.pr()"
]
},
{
"cell_type": "markdown",
"id": "ba265f8a-8b4f-4704-af51-2e9cc45a222c",
"metadata": {},
"source": [
"Note that, till this moment, everything is still abstract. To do the numerical simulation, we need to bring them to a particular implementation, for example, the `msepy`, standing for *mimetic spectral elements python*, by calling `ph.fem.apply` function, which will be shown in other notebooks."
]
}
],
"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.12.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}