Skip to content

Instantly share code, notes, and snippets.

@guyer

guyer/BlockNonDiagonal.ipynb

Last active February 16, 2022 23:22
Show Gist options
  • Save guyer/c71281d62998ae28e9f3ba5050d86699 to your computer and use it in GitHub Desktop.
Save guyer/c71281d62998ae28e9f3ba5050d86699 to your computer and use it in GitHub Desktop.
Download ZIP
Demonstration of FiPy's difficulty with coupled equations when using PETSc
Display the source blob
Display the rendered blob
Raw
BlockNonDiagonal.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "9eded64f-c34e-4622-86ef-2edb252adae0",
"metadata": {},
"source": [
"# Minimal demonstration of FiPy's difficulty with block-non-diagonal equations"
]
},
{
"cell_type": "markdown",
"id": "f7dfbac2-8bd1-4518-b808-d2af4606b366",
"metadata": {},
"source": [
"## Imports"
]
},
{
"cell_type": "markdown",
"id": "e6475ef1-424b-4f86-ba60-b7ac37b9d34e",
"metadata": {},
"source": [
"**note:** The PETSc solvers are unable to solve matrix generated by this problem, as it's block structure is not block-diagonal (see [#840][#840]).\n",
"\n",
"[#840]: https://github.com/usnistgov/fipy/issues/840"
]
},
{
"cell_type": "markdown",
"id": "40352656-3b1a-4c04-a783-99e53535b351",
"metadata": {},
"source": [
"import os\n",
"os.environ['FIPY_SOLVERS'] = 'petsc'"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "429d3cbb-f6fd-4c18-a008-80883f7a4dda",
"metadata": {},
"outputs": [],
"source": [
"import fipy as fp"
]
},
{
"cell_type": "markdown",
"id": "2b257203-5af2-433c-859d-6ea5d2008839",
"metadata": {},
"source": [
"## Domain"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "16867d32-ce36-4734-9391-f7b850087964",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"dx = 1\n",
"nx = 3\n",
"Lx = nx * dx\n",
"mesh = fp.Grid1D(nx=nx, dx=dx)"
]
},
{
"cell_type": "markdown",
"id": "a3d966ad-cb6c-482b-ab0b-d69f239a6817",
"metadata": {},
"source": [
"## Equations"
]
},
{
"cell_type": "markdown",
"id": "3071ae9f-756d-4ad8-b4b6-a880249a1024",
"metadata": {},
"source": [
"The 6th order PDE\n",
"$$\n",
"\\frac{\\partial n}{\\partial t} = \\nabla^2\\left[\n",
" \\left(\n",
" 1 + \\nabla^2 \\left(2 + \\nabla^2\\right)\n",
" \\right)n + n^3\n",
"\\right]\n",
"$$\n",
"can be represented as three 2nd order PDEs"
]
},
{
"cell_type": "markdown",
"id": "f2eab18f-04fa-4642-82b7-dd484cc29815",
"metadata": {},
"source": [
"$$\\begin{align*}\n",
"\\frac{\\partial n}{\\partial t}\n",
"&= \\Gamma\\nabla^2\\mu\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 53,
"id": "49b80234-5387-441a-88b8-79cf0479992e",
"metadata": {},
"outputs": [],
"source": [
"n = fp.CellVariable(mesh=mesh, name=r\"$n$\", hasOld=True)"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "6d11a9e9-6024-489d-af12-c7f4274b1f5b",
"metadata": {},
"outputs": [],
"source": [
"mu = fp.CellVariable(mesh=mesh, name=r\"$\\mu$\")"
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "00d44f92-4093-45a4-bde9-067459a30c3b",
"metadata": {},
"outputs": [],
"source": [
"n_eq = (fp.TransientTerm(coeff=1, var=n) == fp.DiffusionTerm(coeff=1, var=mu))"
]
},
{
"cell_type": "markdown",
"id": "1f23dadd-755d-4577-aaf3-5fd5753ea5d6",
"metadata": {},
"source": [
"where\n",
"\n",
"$$\\mu \\equiv\n",
"\\left(n + \\nabla^2 \\xi\\right)\n",
"+ n^3\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "30a0906d-09a2-4c3e-922d-1635a4382a07",
"metadata": {},
"outputs": [],
"source": [
"xi = fp.CellVariable(mesh=mesh, name=r\"$\\xi$\")"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "ad398c19-283f-4655-b730-5c7d6c495b41",
"metadata": {},
"outputs": [],
"source": [
"mu_eq = (fp.ImplicitSourceTerm(coeff=1, var=mu)\n",
" == fp.ImplicitSourceTerm(coeff=1 + n**2, var=n)\n",
" + fp.DiffusionTerm(coeff=1, var=xi))"
]
},
{
"cell_type": "markdown",
"id": "73b02fcb-bca8-4f2e-ac4c-f4898524de67",
"metadata": {},
"source": [
"and\n",
"\n",
"$$\\xi \\equiv \\left(2 + \\nabla^2\\right)n$$"
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "2c96ec97-049d-4bcc-b979-8ab12bb1c250",
"metadata": {},
"outputs": [],
"source": [
"xi_eq = (fp.ImplicitSourceTerm(coeff=1, var=xi)\n",
" == fp.ImplicitSourceTerm(coeff=2, var=n)\n",
" + fp.DiffusionTerm(coeff=1, var=n))"
]
},
{
"cell_type": "markdown",
"id": "3b056cdb-a0ea-41c5-875b-48e70b7032df",
"metadata": {},
"source": [
"We solve $n$, $\\mu$, and $\\xi$ together"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "9fa7f93b-df72-4cf6-9ef9-b4c2808390f4",
"metadata": {},
"outputs": [],
"source": [
"eq = n_eq & mu_eq & xi_eq"
]
},
{
"cell_type": "markdown",
"id": "6195a945-7913-473b-bb0a-09bca44c036f",
"metadata": {},
"source": [
"## Solution"
]
},
{
"cell_type": "code",
"execution_count": 60,
"id": "b9f94ad4-7719-47d2-80ba-f41a9ecddc95",
"metadata": {},
"outputs": [],
"source": [
"solver = eq._prepareLinearSystem(var=None, solver=None, boundaryConditions=(), dt=1)"
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "580e0d74-3d46-476c-89a1-d19c390bdc12",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.000000 --- --- --- --- --- 1.000000 -1.000000 --- \n",
" --- 1.000000 --- --- --- --- -1.000000 2.000000 -1.000000 \n",
" --- --- 1.000000 --- --- --- --- -1.000000 1.000000 \n",
" --- --- --- 1.000000 -1.000000 --- 1.000000 --- --- \n",
" --- --- --- -1.000000 2.000000 -1.000000 --- 1.000000 --- \n",
" --- --- --- --- -1.000000 1.000000 --- --- 1.000000 \n",
" 1.000000 -1.000000 --- 1.000000 --- --- --- --- --- \n",
"-1.000000 2.000000 -1.000000 --- 1.000000 --- --- --- --- \n",
" --- -1.000000 1.000000 --- --- 1.000000 --- --- --- \n"
]
}
],
"source": [
"print(solver.matrix)"
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "a2a06172-8b22-453f-bd02-2344efce927e",
"metadata": {},
"outputs": [
{
"ename": "Error",
"evalue": "error code 73\n[0] KSPSolve() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/ksp/interface/itfunc.c:1086\n[0] KSPSolve_Private() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/ksp/interface/itfunc.c:852\n[0] KSPSetUp() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/ksp/interface/itfunc.c:408\n[0] PCSetUp() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/pc/interface/precon.c:1016\n[0] PCSetUp_ILU() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/pc/impls/factor/ilu/ilu.c:144\n[0] MatILUFactorSymbolic() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/mat/interface/matrix.c:6943\n[0] MatILUFactorSymbolic_SeqAIJ() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/mat/impls/aij/seq/aijfact.c:1697\n[0] Object is in wrong state\n[0] Matrix is missing diagonal entry 6",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mError\u001b[0m Traceback (most recent call last)",
"Input \u001b[0;32mIn [62]\u001b[0m, in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43meq\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msolve\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdt\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\n",
"File \u001b[0;32m~/Documents/research/FiPy/fipy/fipy/terms/term.py:178\u001b[0m, in \u001b[0;36mTerm.solve\u001b[0;34m(self, var, solver, boundaryConditions, dt)\u001b[0m\n\u001b[1;32m 157\u001b[0m \u001b[38;5;124mr\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 158\u001b[0m \u001b[38;5;124;03mBuilds and solves the `Term`'s linear system once. This method\u001b[39;00m\n\u001b[1;32m 159\u001b[0m \u001b[38;5;124;03mdoes not return the residual. It should be used when the\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 173\u001b[0m \u001b[38;5;124;03m Timestep size.\u001b[39;00m\n\u001b[1;32m 174\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 176\u001b[0m solver \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_prepareLinearSystem(var, solver, boundaryConditions, dt)\n\u001b[0;32m--> 178\u001b[0m \u001b[43msolver\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_solve\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
"File \u001b[0;32m~/Documents/research/FiPy/fipy/fipy/solvers/petsc/petscSolver.py:59\u001b[0m, in \u001b[0;36mPETScSolver._solve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 53\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ((\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmatrix \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m)\n\u001b[1;32m 54\u001b[0m \u001b[38;5;129;01mor\u001b[39;00m (\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmatrix\u001b[38;5;241m.\u001b[39mmatrix\u001b[38;5;241m.\u001b[39msizes[\u001b[38;5;241m0\u001b[39m][\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmatrix\u001b[38;5;241m.\u001b[39mmatrix\u001b[38;5;241m.\u001b[39msizes[\u001b[38;5;241m1\u001b[39m][\u001b[38;5;241m1\u001b[39m])\n\u001b[1;32m 55\u001b[0m \u001b[38;5;129;01mor\u001b[39;00m (\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmatrix\u001b[38;5;241m.\u001b[39mmatrix\u001b[38;5;241m.\u001b[39msizes[\u001b[38;5;241m0\u001b[39m][\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m!=\u001b[39m overlappingVector\u001b[38;5;241m.\u001b[39msize)):\n\u001b[1;32m 57\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m SolutionVariableNumberError\n\u001b[0;32m---> 59\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_solve_\u001b[49m\u001b[43m(\u001b[49m\u001b[43mglobalMatrix\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmatrix\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 60\u001b[0m \u001b[43m \u001b[49m\u001b[43moverlappingVector\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 61\u001b[0m \u001b[43m \u001b[49m\u001b[43moverlappingRHSvector\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 63\u001b[0m value \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmatrix\u001b[38;5;241m.\u001b[39m_petsc2fipyGhost(vec\u001b[38;5;241m=\u001b[39moverlappingVector)\n\u001b[1;32m 64\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mvar\u001b[38;5;241m.\u001b[39mvalue \u001b[38;5;241m=\u001b[39m numerix\u001b[38;5;241m.\u001b[39mreshape(value, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mvar\u001b[38;5;241m.\u001b[39mshape)\n",
"File \u001b[0;32m~/Documents/research/FiPy/fipy/fipy/solvers/petsc/petscKrylovSolver.py:64\u001b[0m, in \u001b[0;36mPETScKrylovSolver._solve_\u001b[0;34m(self, L, x, b)\u001b[0m\n\u001b[1;32m 62\u001b[0m ksp\u001b[38;5;241m.\u001b[39msetOperators(L)\n\u001b[1;32m 63\u001b[0m ksp\u001b[38;5;241m.\u001b[39msetFromOptions()\n\u001b[0;32m---> 64\u001b[0m \u001b[43mksp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msolve\u001b[49m\u001b[43m(\u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 66\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mFIPY_VERBOSE_SOLVER\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;129;01min\u001b[39;00m os\u001b[38;5;241m.\u001b[39menviron:\n\u001b[1;32m 67\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mfipy\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtools\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mdebug\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m PRINT\n",
"File \u001b[0;32mPETSc/KSP.pyx:397\u001b[0m, in \u001b[0;36mpetsc4py.PETSc.KSP.solve\u001b[0;34m()\u001b[0m\n",
"\u001b[0;31mError\u001b[0m: error code 73\n[0] KSPSolve() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/ksp/interface/itfunc.c:1086\n[0] KSPSolve_Private() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/ksp/interface/itfunc.c:852\n[0] KSPSetUp() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/ksp/interface/itfunc.c:408\n[0] PCSetUp() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/pc/interface/precon.c:1016\n[0] PCSetUp_ILU() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/ksp/pc/impls/factor/ilu/ilu.c:144\n[0] MatILUFactorSymbolic() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/mat/interface/matrix.c:6943\n[0] MatILUFactorSymbolic_SeqAIJ() at /Users/runner/miniforge3/conda-bld/petsc_1637321992139/work/src/mat/impls/aij/seq/aijfact.c:1697\n[0] Object is in wrong state\n[0] Matrix is missing diagonal entry 6"
]
}
],
"source": [
"eq.solve(dt=1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f3e48d5d-4eb2-4aa6-b70f-a4601efab6e9",
"metadata": {},
"outputs": [],
"source": []
}
],
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment