FUNC_PREFIX void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT _data_u2)\n",
"{\n",
" for (int ctr_0 = 1; ctr_0 < 59; ctr_0 += 1)\n",
" {\n",
" double * RESTRICT _data_u2_00 = _data_u2 + 70*ctr_0;\n",
" double * RESTRICT const _data_u1_00 = _data_u1 + 70*ctr_0;\n",
" double * RESTRICT const _data_u0_00 = _data_u0 + 70*ctr_0;\n",
" double * RESTRICT const _data_u1_01 = _data_u1 + 70*ctr_0 + 70;\n",
" double * RESTRICT const _data_u1_0m1 = _data_u1 + 70*ctr_0 - 70;\n",
" for (int ctr_1 = 1; ctr_1 < 69; ctr_1 += 1)\n",
" {\n",
" _data_u2_00[ctr_1] = 0.25*_data_u1_00[ctr_1 + 1] + 0.25*_data_u1_00[ctr_1 - 1] + 0.25*_data_u1_01[ctr_1] + 0.25*_data_u1_0m1[ctr_1] - 1.0*_data_u0_00[ctr_1] + 1.0*_data_u1_00[ctr_1];\n",
" }\n",
" }\n",
"}\n",
"
\n"
],
"text/plain": [
"\n",
"FUNC_PREFIX void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT _data_u2)\n",
"{\n",
" for (int ctr_0 = 1; ctr_0 < 59; ctr_0 += 1)\n",
" {\n",
" double * RESTRICT _data_u2_00 = _data_u2 + 70*ctr_0;\n",
" double * RESTRICT const _data_u1_00 = _data_u1 + 70*ctr_0;\n",
" double * RESTRICT const _data_u0_00 = _data_u0 + 70*ctr_0;\n",
" double * RESTRICT const _data_u1_01 = _data_u1 + 70*ctr_0 + 70;\n",
" double * RESTRICT const _data_u1_0m1 = _data_u1 + 70*ctr_0 - 70;\n",
" for (int ctr_1 = 1; ctr_1 < 69; ctr_1 += 1)\n",
" {\n",
" _data_u2_00[ctr_1] = 0.25*_data_u1_00[ctr_1 + 1] + 0.25*_data_u1_00[ctr_1 - 1] + 0.25*_data_u1_01[ctr_1] + 0.25*_data_u1_0m1[ctr_1] - 1.0*_data_u0_00[ctr_1] + 1.0*_data_u1_00[ctr_1];\n",
" }\n",
" }\n",
"}"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"update_rule = update_rule.subs({discretize.dx: 0.1, discretize.dt: 0.05})\n",
"ast = ps.create_kernel(update_rule)\n",
"kernel = ast.compile()\n",
"\n",
"ps.show_code(ast)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[_data_u0, _data_u1, _data_u2]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ast.get_parameters()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{u0, u1, u2}"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ast.fields_accessed"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To run simulation a suitable initial condition and boundary treatment is required. We chose an initial condition which is zero at the borders of the domain. The outermost layer is not changed by the update kernel, so we have an implicit homogenous Dirichlet boundary condition."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"X,Y = np.meshgrid( np.linspace(0, 1, size[1]), np.linspace(0,1, size[0]))\n",
"Z = np.sin(2*X*np.pi) * np.sin(2*Y*np.pi)\n",
"\n",
"# Initialize the previous and current values with the initial function\n",
"np.copyto(u_arrays[0], Z)\n",
"np.copyto(u_arrays[1], Z)\n",
"# The values for the next timesteps do not matter, since they are overwritten\n",
"u_arrays[2][:, :] = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One timestep now consists of applying the kernel once, then shifting the arrays."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def run(timesteps=1):\n",
" for t in range(timesteps):\n",
" kernel(u0=u_arrays[0], u1=u_arrays[1], u2=u_arrays[2])\n",
" u_arrays[0], u_arrays[1], u_arrays[2] = u_arrays[1], u_arrays[2], u_arrays[0]\n",
" return u_arrays[2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lets create an animation of the solution:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"