Skip to content

Commit 6bf06c8

Browse files
committed
What is the actual biharmonic 2d green's function?
1 parent 954efe1 commit 6bf06c8

9 files changed

Lines changed: 266 additions & 1653 deletions

sumpy/recurrence.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -277,7 +277,7 @@ def recurrence_from_pde(pde: LinearPDESystemOperator) -> sp.Expr:
277277
"""
278278
ode_in_r, var, ode_order = pde_to_ode_in_r(pde)
279279
ode_in_x = ode_in_r_to_x(ode_in_r, var, ode_order).simplify()
280-
ode_in_x_cleared = (ode_in_x * var[0]**(ode_order+1)).simplify()
280+
ode_in_x_cleared = (ode_in_x * var[0]**(pde.order*2-1)).simplify()
281281
# ode_in_x_cleared shouldn't have rational function coefficients
282282
assert sp.together(ode_in_x_cleared) == ode_in_x_cleared
283283
f_x_derivs = _make_sympy_vec("f_x", ode_order+1)

test/biharmonic.ipynb

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 3,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import numpy as np\n",
10+
"import sympy as sp\n",
11+
"\n",
12+
"from sumpy.expansion.diff_op import (\n",
13+
" make_identity_diff_op,\n",
14+
")\n",
15+
"from collections import namedtuple\n",
16+
"DerivativeIdentifier = namedtuple(\"DerivativeIdentifier\", [\"mi\", \"vec_idx\"])\n",
17+
"\n",
18+
"from sumpy.recurrence import _make_sympy_vec, get_reindexed_and_center_origin_on_axis_recurrence\n",
19+
"\n",
20+
"from immutabledict import immutabledict\n",
21+
"from sumpy.expansion.diff_op import LinearPDESystemOperator"
22+
]
23+
},
24+
{
25+
"cell_type": "code",
26+
"execution_count": 5,
27+
"metadata": {},
28+
"outputs": [],
29+
"source": [
30+
"var = _make_sympy_vec(\"x\", 2)\n",
31+
"var_t = _make_sympy_vec(\"t\", 2)\n",
32+
"abs_dist = sp.sqrt((var[0]-var_t[0])**2 + (var[1]-var_t[1])**2)\n",
33+
"w = make_identity_diff_op(2)\n",
34+
"\n",
35+
"partial_1x = DerivativeIdentifier((4,0), 0)\n",
36+
"partial_1y = DerivativeIdentifier((0,4), 0)\n",
37+
"biharmonic_op = {partial_1x: 1, partial_1y: 1}\n",
38+
"list_pde = immutabledict(biharmonic_op)\n",
39+
"\n",
40+
"biharmonic_pde = LinearPDESystemOperator(2, (list_pde,))\n",
41+
"g_x_y = abs_dist**2 * (sp.log(abs_dist)-1)\n",
42+
"\n",
43+
"n_init, _, r = get_reindexed_and_center_origin_on_axis_recurrence(biharmonic_pde)\n",
44+
"\n",
45+
"derivs = [sp.diff(g_x_y,\n",
46+
" var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0)\n",
47+
" for i in range(8)]\n",
48+
"\n",
49+
"x_coord = np.random.rand() # noqa: NPY002\n",
50+
"y_coord = np.random.rand() # noqa: NPY002\n",
51+
"coord_dict = {var[0]: var[0], var[1]: var[1]}\n",
52+
"derivs = [d.subs(coord_dict) for d in derivs]\n",
53+
"\n",
54+
"n = sp.symbols(\"n\")\n",
55+
"s = sp.Function(\"s\")\n",
56+
"\n",
57+
"# pylint: disable-next=not-callable\n",
58+
"subs_dict = {s(0): derivs[0], s(1): derivs[1], s(2): derivs[1], s(3): derivs[1]}\n",
59+
"check = []\n",
60+
"\n",
61+
"assert n_init == 4\n",
62+
"max_order_check = 8\n",
63+
"for i in range(n_init, max_order_check):\n",
64+
" check.append(r.subs(n, i).subs(subs_dict) - derivs[i])\n",
65+
" # pylint: disable-next=not-callable\n",
66+
" subs_dict[s(i)] = derivs[i]"
67+
]
68+
},
69+
{
70+
"cell_type": "code",
71+
"execution_count": 7,
72+
"metadata": {},
73+
"outputs": [],
74+
"source": [
75+
"max_order_check = 8\n",
76+
"for i in range(n_init, max_order_check):\n",
77+
" check.append(r.subs(n, i).subs(subs_dict) - derivs[i])\n",
78+
" # pylint: disable-next=not-callable\n",
79+
" subs_dict[s(i)] = derivs[i]"
80+
]
81+
},
82+
{
83+
"cell_type": "code",
84+
"execution_count": 14,
85+
"metadata": {},
86+
"outputs": [
87+
{
88+
"data": {
89+
"text/latex": [
90+
"$\\displaystyle 16.9964566798618$"
91+
],
92+
"text/plain": [
93+
"16.9964566798618"
94+
]
95+
},
96+
"execution_count": 14,
97+
"metadata": {},
98+
"output_type": "execute_result"
99+
}
100+
],
101+
"source": [
102+
"r.subs(n, 4).subs(subs_dict).subs({var[0]: 1.2, var[1]: 2.3})"
103+
]
104+
},
105+
{
106+
"cell_type": "code",
107+
"execution_count": 6,
108+
"metadata": {},
109+
"outputs": [
110+
{
111+
"data": {
112+
"text/plain": [
113+
"{s(0): (x0**2 + x1**2)*(log(sqrt(x0**2 + x1**2)) - 1),\n",
114+
" s(1): -2*x0*(log(sqrt(x0**2 + x1**2)) - 1) - x0,\n",
115+
" s(2): -2*x0*(log(sqrt(x0**2 + x1**2)) - 1) - x0,\n",
116+
" s(3): -2*x0*(log(sqrt(x0**2 + x1**2)) - 1) - x0,\n",
117+
" s(4): 2*(-24*x0**4/(x0**2 + x1**2)**2 + 8*x0**2*(4*x0**2/(x0**2 + x1**2) - 3)/(x0**2 + x1**2) + 12*x0**2/(x0**2 + x1**2) + 3)/(x0**2 + x1**2),\n",
118+
" s(5): -4*x0*(-24*x0**4/(x0**2 + x1**2)**2 + 40*x0**2/(x0**2 + x1**2) - 15)/(x0**2 + x1**2)**2,\n",
119+
" s(6): 12*(-320*x0**6/(x0**2 + x1**2)**3 + 360*x0**4/(x0**2 + x1**2)**2 + 24*x0**2*(16*x0**4/(x0**2 + x1**2)**2 - 20*x0**2/(x0**2 + x1**2) + 5)/(x0**2 + x1**2) - 60*x0**2/(x0**2 + x1**2) - 5)/(x0**2 + x1**2)**2,\n",
120+
" s(7): -48*x0*(-160*x0**6/(x0**2 + x1**2)**3 + 336*x0**4/(x0**2 + x1**2)**2 - 210*x0**2/(x0**2 + x1**2) + 35)/(x0**2 + x1**2)**3}"
121+
]
122+
},
123+
"execution_count": 6,
124+
"metadata": {},
125+
"output_type": "execute_result"
126+
}
127+
],
128+
"source": [
129+
"subs_dict"
130+
]
131+
},
132+
{
133+
"cell_type": "code",
134+
"execution_count": null,
135+
"metadata": {},
136+
"outputs": [],
137+
"source": []
138+
}
139+
],
140+
"metadata": {
141+
"kernelspec": {
142+
"display_name": "inteq",
143+
"language": "python",
144+
"name": "python3"
145+
},
146+
"language_info": {
147+
"codemirror_mode": {
148+
"name": "ipython",
149+
"version": 3
150+
},
151+
"file_extension": ".py",
152+
"mimetype": "text/x-python",
153+
"name": "python",
154+
"nbconvert_exporter": "python",
155+
"pygments_lexer": "ipython3",
156+
"version": "3.11.9"
157+
}
158+
},
159+
"nbformat": 4,
160+
"nbformat_minor": 2
161+
}

test/gridfree_taylor_recurrence.ipynb

Lines changed: 0 additions & 167 deletions
This file was deleted.

0 commit comments

Comments
 (0)