Skip to content

Instantly share code, notes, and snippets.

@guyer

guyer/Constraints&BCs.ipynb

Last active January 11, 2021 16:40
Show Gist options
  • Save guyer/2bd27a25bdd40b4aa26e40222cf3d046 to your computer and use it in GitHub Desktop.
Save guyer/2bd27a25bdd40b4aa26e40222cf3d046 to your computer and use it in GitHub Desktop.
Download ZIP
Illustration of applying boundary conditions in FiPy
Display the source blob
Display the rendered blob
Raw
Constraints&BCs.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T13:52:01.407921Z",
"start_time": "2021-01-11T13:51:58.929584Z"
}
},
"outputs": [],
"source": [
"import fipy as fp"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T13:52:12.697567Z",
"start_time": "2021-01-11T13:52:12.676331Z"
}
},
"outputs": [],
"source": [
"mesh = fp.Grid1D(nx=5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## No `BoundaryCondition`s or constraints"
]
},
{
"cell_type": "code",
"execution_count": 118,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T16:27:43.909197Z",
"start_time": "2021-01-11T16:27:43.864361Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.000000 -1.000000 --- --- --- \n",
"-1.000000 3.000000 -1.000000 --- --- \n",
" --- -1.000000 3.000000 -1.000000 --- \n",
" --- --- -1.000000 3.000000 -1.000000 \n",
" --- --- --- -1.000000 2.000000 \n",
"[[0.]\n",
" [0.]\n",
" [0.]\n",
" [0.]\n",
" [0.]]\n",
"solution: [0. 0. 0. 0. 0.]\n",
"face value: [0. 0. 0. 0. 0. 0.]\n",
"face gradient: [[0. 0. 0. 0. 0. 0.]]\n"
]
}
],
"source": [
"var = fp.CellVariable(mesh=mesh)\n",
"eq = fp.TransientTerm() == fp.DiffusionTerm()\n",
"\n",
"solver = eq._prepareLinearSystem(var=var, solver=None, boundaryConditions=(), dt=1.)\n",
"\n",
"print(solver.matrix)\n",
"print(solver.RHSvector[..., fp.numerix.newaxis])\n",
"\n",
"eq.solve(var=var, dt=1.)\n",
"\n",
"print(\"solution:\", var)\n",
"print(\"face value:\", var.faceValue)\n",
"print(\"face gradient:\", var.faceGrad)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dirichlet conditions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### `FixedValue`\n",
"\n",
"Note: A lot of subtle caching goes on, so we completely redefine the problem each time."
]
},
{
"cell_type": "code",
"execution_count": 120,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T16:28:41.592240Z",
"start_time": "2021-01-11T16:28:41.499733Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 4.000000 -1.000000 --- --- --- \n",
"-1.000000 3.000000 -1.000000 --- --- \n",
" --- -1.000000 3.000000 -1.000000 --- \n",
" --- --- -1.000000 3.000000 -1.000000 \n",
" --- --- --- -1.000000 2.000000 \n",
"[[10.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [ 0.]]\n",
"solution: [2.76422764 1.05691057 0.40650407 0.16260163 0.08130081]\n",
"face value: [2.76422764 1.91056911 0.73170732 0.28455285 0.12195122 0.08130081]\n",
"face gradient: [[ 0. -1.70731707 -0.6504065 -0.24390244 -0.08130081 0. ]]\n"
]
}
],
"source": [
"var = fp.CellVariable(mesh=mesh)\n",
"eq = fp.TransientTerm() == fp.DiffusionTerm()\n",
"\n",
"BCs = (fp.FixedValue(value=5., faces=mesh.facesLeft),)\n",
"solver = eq._prepareLinearSystem(var=var, solver=None, boundaryConditions=BCs, dt=1.)\n",
"\n",
"print(solver.matrix)\n",
"print(solver.RHSvector[..., fp.numerix.newaxis])\n",
"\n",
"eq.solve(var=var, boundaryConditions=BCs, dt=1.)\n",
"\n",
"print(\"solution:\", var)\n",
"print(\"face value:\", var.faceValue)\n",
"print(\"face gradient:\", var.faceGrad)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the left-hand face value is identical to the left-hand cell value. This is a result of the zero-exterior-gradient assumption intrinsic to cell-centered finite volume. The variable and its gradient are oblivious to the boundary conditions used to solve the problem."
]
},
{
"cell_type": "markdown",
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T16:02:18.053417Z",
"start_time": "2021-01-11T16:02:18.051189Z"
}
},
"source": [
"### `.value` constraint"
]
},
{
"cell_type": "code",
"execution_count": 121,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T16:29:15.776934Z",
"start_time": "2021-01-11T16:29:15.728067Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 4.000000 -1.000000 --- --- --- \n",
"-1.000000 3.000000 -1.000000 --- --- \n",
" --- -1.000000 3.000000 -1.000000 --- \n",
" --- --- -1.000000 3.000000 -1.000000 \n",
" --- --- --- -1.000000 2.000000 \n",
"[[10.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [ 0.]]\n",
"solution: [2.76422764 1.05691057 0.40650407 0.16260163 0.08130081]\n",
"face value: [5. 1.91056911 0.73170732 0.28455285 0.12195122 0.08130081]\n",
"face gradient: [[-4.47154472 -1.70731707 -0.6504065 -0.24390244 -0.08130081 0. ]]\n"
]
}
],
"source": [
"var = fp.CellVariable(mesh=mesh)\n",
"eq = fp.TransientTerm() == fp.DiffusionTerm()\n",
"\n",
"var.constrain(5., where=mesh.facesLeft)\n",
"\n",
"solver = eq._prepareLinearSystem(var=var, solver=None, boundaryConditions=(), dt=1.)\n",
"\n",
"print(solver.matrix)\n",
"print(solver.RHSvector[..., fp.numerix.newaxis])\n",
"\n",
"eq.solve(var=var, dt=1.)\n",
"\n",
"print(\"solution:\", var)\n",
"print(\"face value:\", var.faceValue)\n",
"print(\"face gradient:\", var.faceGrad)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The left-hand face value matches the constraint, and has the correct corresponding face gradient , but the solution is otherwise identical to that found above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Neumann conditions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### `FixedFlux`"
]
},
{
"cell_type": "code",
"execution_count": 122,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T16:29:41.319784Z",
"start_time": "2021-01-11T16:29:41.275165Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.000000 -1.000000 --- --- --- \n",
"-1.000000 3.000000 -1.000000 --- --- \n",
" --- -1.000000 3.000000 -1.000000 --- \n",
" --- --- -1.000000 3.000000 -1.000000 \n",
" --- --- --- -1.000000 2.000000 \n",
"[[ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [-3.]]\n",
"solution: [-0.05454545 -0.10909091 -0.27272727 -0.70909091 -1.85454545]\n",
"face value: [-0.05454545 -0.08181818 -0.19090909 -0.49090909 -1.28181818 -1.85454545]\n",
"face gradient: [[ 0. -0.05454545 -0.16363636 -0.43636364 -1.14545455 0. ]]\n"
]
}
],
"source": [
"var = fp.CellVariable(mesh=mesh)\n",
"eq = fp.TransientTerm() == fp.DiffusionTerm()\n",
"\n",
"BCs = (fp.FixedFlux(value=3., faces=mesh.facesRight),)\n",
"solver = eq._prepareLinearSystem(var=var, solver=None, boundaryConditions=BCs, dt=1.)\n",
"\n",
"print(solver.matrix)\n",
"print(solver.RHSvector[..., fp.numerix.newaxis])\n",
"\n",
"eq.solve(var=var, boundaryConditions=BCs, dt=1.)\n",
"\n",
"print(\"solution:\", var)\n",
"print(\"face value:\", var.faceValue)\n",
"print(\"face gradient:\", var.faceGrad)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While the solution and face value are as expected, the face gradient does not show signs of the boundary condition."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Boundary flux source"
]
},
{
"cell_type": "code",
"execution_count": 123,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T16:29:46.669452Z",
"start_time": "2021-01-11T16:29:46.615066Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.000000 -1.000000 --- --- --- \n",
"-1.000000 3.000000 -1.000000 --- --- \n",
" --- -1.000000 3.000000 -1.000000 --- \n",
" --- --- -1.000000 3.000000 -1.000000 \n",
" --- --- --- -1.000000 2.000000 \n",
"[[ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [-3.]]\n",
"solution: [-0.05454545 -0.10909091 -0.27272727 -0.70909091 -1.85454545]\n",
"face value: [-0.05454545 -0.08181818 -0.19090909 -0.49090909 -1.28181818 -1.85454545]\n",
"face gradient: [[ 0. -0.05454545 -0.16363636 -0.43636364 -1.14545455 0. ]]\n"
]
}
],
"source": [
"var = fp.CellVariable(mesh=mesh)\n",
"eq = fp.TransientTerm() == fp.DiffusionTerm() + (mesh.facesRight * -3 * mesh.faceNormals).divergence\n",
"\n",
"solver = eq._prepareLinearSystem(var=var, solver=None, boundaryConditions=(), dt=1.)\n",
"\n",
"print(solver.matrix)\n",
"print(solver.RHSvector[..., fp.numerix.newaxis])\n",
"\n",
"eq.solve(var=var, dt=1.)\n",
"\n",
"print(\"solution:\", var)\n",
"print(\"face value:\", var.faceValue)\n",
"print(\"face gradient:\", var.faceGrad)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution is identical to that above and the boundary gradient is insensitive to the source placed at that boundary."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### `.faceGrad` constraint\n",
"\n",
"For this trivial problem, the relationship between the flux and the gradient is explicit,\n",
"$\\vec{J} = -D \\nabla \\phi$, so we can invert and constrain the face gradient."
]
},
{
"cell_type": "code",
"execution_count": 124,
"metadata": {
"ExecuteTime": {
"end_time": "2021-01-11T16:29:51.581209Z",
"start_time": "2021-01-11T16:29:51.537713Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.000000 -1.000000 --- --- --- \n",
"-1.000000 3.000000 -1.000000 --- --- \n",
" --- -1.000000 3.000000 -1.000000 --- \n",
" --- --- -1.000000 3.000000 -1.000000 \n",
" --- --- --- -1.000000 2.000000 \n",
"[[ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [ 0.]\n",
" [-3.]]\n",
"solution: [-0.05454545 -0.10909091 -0.27272727 -0.70909091 -1.85454545]\n",
"face value: [-0.05454545 -0.08181818 -0.19090909 -0.49090909 -1.28181818 -1.85454545]\n",
"face gradient: [[ 0. -0.05454545 -0.16363636 -0.43636364 -1.14545455 -3. ]]\n"
]
}
],
"source": [
"var = fp.CellVariable(mesh=mesh)\n",
"eq = fp.TransientTerm() == fp.DiffusionTerm()\n",
"\n",
"var.faceGrad.constrain(-3 * mesh.faceNormals, where=mesh.facesRight)\n",
"\n",
"solver = eq._prepareLinearSystem(var=var, solver=None, boundaryConditions=(), dt=1.)\n",
"\n",
"print(solver.matrix)\n",
"print(solver.RHSvector[..., fp.numerix.newaxis])\n",
"\n",
"eq.solve(var=var, dt=1.)\n",
"\n",
"print(\"solution:\", var)\n",
"print(\"face value:\", var.faceValue)\n",
"print(\"face gradient:\", var.faceGrad)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution is everywhere the same as the last two cases, with the exception of the boundary gradient, which reflects the constraint put on it."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6.7"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment