Skip to content

Instantly share code, notes, and snippets.

@guyer

guyer/BlockNonDiagonal.ipynb

Last active February 16, 2022 23:22
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save guyer/c71281d62998ae28e9f3ba5050d86699 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
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