<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=us-ascii"><meta name=Generator content="Microsoft Word 14 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Consolas;
panose-1:2 11 6 9 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri","sans-serif";
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=EN-US link=blue vlink=purple><div class=WordSection1><p class=MsoNormal>Dear PETSc Users/Developers.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>I have been successfully using PETsc on Windows without MPI for a while now. I have now attempted to implement PETSc with MPI on Windows 10. I have built a release version of PETSc 3.18.6 with MS MPI 10.1.2, Intel MKL 3.279 (2020) and Visual Studio 2019 as a static library. I am testing a finite difference 3D Poisson solver (cylindrical coordinates) with this MPI PETSc.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>On a Windows 10 laptop with 4 cores (8 logical processors), the solver works fine with one core. However, it gives errors with 2 cores. The errors are listed below followed by a complete code. The errors are generated in KSPSolve. I perhaps miss some settings to cause this error “No method muldiagonalblock for Mat of type seqaij”. I don’t know why for a 2-core MPI program that a seqaij matrix is used in KSPSolve.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>I would be very grateful if someone can provide me with some direction to fix these issues.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Many thanks for your help,<o:p></o:p></p><p class=MsoNormal>Thuc Bui<o:p></o:p></p><p class=MsoNormal>Senior R&D Engineer<o:p></o:p></p><p class=MsoNormal>Calabazas Creek Research, Inc<o:p></o:p></p><p class=MsoNormal>(650) 948-5361<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>///////////////errors/////////////////////////<o:p></o:p></p><p class=MsoNormal>rank=0 iStart=0 iEnd=2375<o:p></o:p></p><p class=MsoNormal>rank=1 iStart=2375 iEnd=4750<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>--------------------- Error Message --------------------------------------------------------------<o:p></o:p></p><p class=MsoNormal>--------------------- Error Message --------------------------------------------------------------<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>No support for this operation for this object type<o:p></o:p></p><p class=MsoNormal>No support for this operation for this object type<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>No method multdiagonalblock for Mat of type seqaij<o:p></o:p></p><p class=MsoNormal>No method multdiagonalblock for Mat of type seqaij<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>See https://petsc.org/release/faq/ for trouble shooting.<o:p></o:p></p><p class=MsoNormal>See https://petsc.org/release/faq/ for trouble shooting.<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>Petsc Release Version 3.18.6, Mar 30, 2023 <o:p></o:p></p><p class=MsoNormal>Petsc Release Version 3.18.6, Mar 30, 2023 <o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>Poisson3D.exe on a v140-x64SReleaseSP named HERMES by bui Mon Jul 24 15:16:55 2023<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: Configure options --with-cc="win32fe cl" --with-fc=0 --with-cxx="win32fe cl" --with-debugging=0 --with-openmp --with-mpi-include=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include:/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include/x64 --with-mpi-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpi.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpifec.lib]" --with-blaslapack-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_lp64.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_core.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_thread.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/libiomp5md.lib]" --with-mpiexec=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/bin/x64/mpiexec --with-shared-libraries=0<o:p></o:p></p><p class=MsoNormal>Poisson3D.exe on a v140-x64SReleaseSP named HERMES by bui Mon Jul 24 15:16:55 2023<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#1 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<o:p></o:p></p><p class=MsoNormal>Configure options --with-cc="win32fe cl" --with-fc=0 --with-cxx="win32fe cl" --with-debugging=0 --with-openmp --with-mpi-include=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include:/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include/x64 --with-mpi-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpi.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpifec.lib]" --with-blaslapack-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_lp64.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_core.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_thread.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/libiomp5md.lib]" --with-mpiexec=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/bin/x64/mpiexec --with-shared-libraries=0<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#2 MatMultDiagonalBlock_MPIAIJ() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\impls\aij\mpi\mpiaij.c:978<o:p></o:p></p><p class=MsoNormal>#1 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#3 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<o:p></o:p></p><p class=MsoNormal>#2 MatMultDiagonalBlock_MPIAIJ() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\impls\aij\mpi\mpiaij.c:978<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#4 PCApply_Eisenstat() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\impls\eisens\eisen.c:51<o:p></o:p></p><p class=MsoNormal>#3 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#5 PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\INTERF~1\precon.c:441<o:p></o:p></p><p class=MsoNormal>#4 PCApply_Eisenstat() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\impls\eisens\eisen.c:51<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#6 KSP_PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\include\petsc/private/kspimpl.h:380<o:p></o:p></p><p class=MsoNormal>#5 PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\INTERF~1\precon.c:441<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#7 KSPInitialResidual() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itres.c:64<o:p></o:p></p><p class=MsoNormal>#6 KSP_PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\include\petsc/private/kspimpl.h:380<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#8 KSPSolve_GMRES() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\impls\gmres\gmres.c:227<o:p></o:p></p><p class=MsoNormal>#7 KSPInitialResidual() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itres.c:64<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#9 KSPSolve_Private() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:899<o:p></o:p></p><p class=MsoNormal>#8 KSPSolve_GMRES() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\impls\gmres\gmres.c:227<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#10 KSPSolve() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:1071<o:p></o:p></p><p class=MsoNormal>#9 KSPSolve_Private() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:899<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>#11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253<o:p></o:p></p><p class=MsoNormal>#10 KSPSolve() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:1071<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>No PETSc Option Table entries<o:p></o:p></p><p class=MsoNormal>#11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253<o:p></o:p></p><p class=MsoNormal>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal>----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov----------<o:p></o:p></p><p class=MsoNormal>No PETSc Option Table entries<o:p></o:p></p><p class=MsoNormal>[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to <a href="mailto:petsc-maint@mcs.anl.gov----------">petsc-maint@mcs.anl.gov----------</a><o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>//complete codes///////////////////////////////////////////////<o:p></o:p></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:gray'>#include</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"petscsys.h"</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:gray'>#include</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"petscksp.h"</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:gray'>#ifndef</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> PETSC_SUCCESS<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:gray'>#define</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:#6F008A'>PETSC_SUCCESS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> 0<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:gray'>#endif</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:green'>//rectangular matrix mxn</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>** allocateIMatrixR(</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>m</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>n</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>)<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>** v = (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>**)malloc(</span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>m</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> * </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>sizeof</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>(</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>*));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> i = 0; i < </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>m</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>; ++i) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> v[i] = (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>*)malloc(</span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>n</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> * </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>sizeof</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>(</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> i = 0; i < </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>m</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>; ++i) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> j = 0; j < </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>n</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>; ++</span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>n</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> v[i][j] = -1;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>return</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> v;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> mapIndex(</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>k</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>i</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>j</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>Nz</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>Nr</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>)<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//k for z, i for r, j for theta</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//base 0</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//3D table: theta r z</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// 0 0 0</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// 0 0 1</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// 0 0 2</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// ...</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// 0 1 0</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// 0 1 1</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// 0 1 2</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>// ...</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> n = </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>k</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> + </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>i</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> * </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>Nz</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> + </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>j</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> * </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>Nz</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> * </span><span style='font-size:10.0pt;font-family:Consolas;color:gray'>Nr</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>return</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> n;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>PetscErrorCode MPIAssembler(Mat* A, Vec* b, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Kz, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Ir, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Jt)<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscFunctionBeginUser;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> M;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecGetSize(*b, &M));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//contruct an array of Mx3, row index is the equation number, each row is (i,j,k), </span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//i is r-index, j theta-index, k z-index</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>const</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> dim = 3;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>** Nijk = allocateIMatrixR(M, dim);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> k = 0; k < Kz; ++k) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> k1 = k + 1;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> i = 0; i < Ir; ++i) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> j = 0; j < Jt; ++j) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//for each equation row m</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> m = mapIndex(k, i, j, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nijk[m][0] = i;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nijk[m][1] = j;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nijk[m][2] = k;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>/////////////////////////////////////////////////////////////////////////////////</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//Cylinder dimensions</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> R = 11.0; </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//cm, cylinder radius</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> L = 11.0; </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//cm, cylinder length in z</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Dr = R / (Ir + 1);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Dz = L / (Kz + 1);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Dr2 = Dr * Dr;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Dr_2 = 0.5 * Dr;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Dr_22 = Dr_2 * Dr_2;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Dz_2 = 1.0 / (Dz * Dz);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> DT = 2.0 * M_PI / Jt;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> DT2 = DT * DT;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Vr = 0.0; </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//Dirichlet on R</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> chargeDens = -4.0;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> rank;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//get the range of matrix rows owned by this processor</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> startingNumber, endingNumber;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatGetOwnershipRange(*A, &startingNumber, &endingNumber));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> fprintf(stderr, </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"rank=%i iStart=%i iEnd=%i\n"</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, rank, startingNumber, endingNumber);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Ir_1 = Ir - 1;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Jt_1 = Jt - 1;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> m = startingNumber; m < endingNumber; ++m) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> i = Nijk[m][0];<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> j = Nijk[m][1];<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> k = Nijk[m][2];<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> ri = Dr * (i + 1);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> ri_2 = ri - Dr_2;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> ri2 = ri + Dr_2;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> sqri = ri * ri;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Ai = (ri_2 / ri) / Dr2;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Bi = (ri2 / ri) / Dr2;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Ci = 1.0 / (sqri * DT2);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Ei = (ri2 + ri_2) / (ri * Dr2) + 2.0 / (sqri * DT2) + 2.0 * Dz_2;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//contribution of charge density to RHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetValue(*b, m, -chargeDens, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//contribution to RHS by Drichlet on R</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>if</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (i == Ir_1) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetValue(*b, i, -Bi * Vr, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>if</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (!i) { </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//from the axis</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//contributing to RHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetValue(*b, m, -chargeDens * Ai * Dr_22, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//contributing to LHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> xTheta = (Ai / Jt);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> jt = 0; jt < Jt; ++jt) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> nj = mapIndex(k, i, jt, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &nj, &xTheta, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>else</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//0 < i < Ir: contribution of Ai to LHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> i_jk = mapIndex(k, i - 1, j, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &i_jk, &Ai, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>if</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (i != Ir_1) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//0 <= i < It-1: contribution of Bi to LHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> ki1j = mapIndex(k, i + 1, j, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &ki1j, &Bi, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//apply periodic BC in theta////////////////////////////</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij_;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>if</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (!j) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij_ = mapIndex(k, i, Jt_1, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>else</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij_ = mapIndex(k, i, j - 1, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//1st contribution of Ci to LHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &kij_, &Ci, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij1;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>if</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (j == Jt_1) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij1 = mapIndex(k, i, 0, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>else</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij1 = mapIndex(k, i, j + 1, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//2nd contribution of Ci to LHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &kij1, &Ci, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>////////////////////////////////////////////////</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//apply z Neumann BC ///////////////////////////////////</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> k_ij;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>if</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (!k) { </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//z=0</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> k_ij = mapIndex(k, i, j, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>else</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> { </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//0 < z < L</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> k_ij = mapIndex(k - 1, i, j, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &k_ij, &Dz_2, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>/////////////////////////////////////////////////////////</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//contribution of Ei to LHS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij = mapIndex(k, i, j, Kz, Ir);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &kij, &Ei, ADD_VALUES));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>if</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (Nijk) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>for</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> (</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> i = 0; i < M; ++i) {<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> free(Nijk[i]);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> free(Nijk);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecAssemblyBegin(*b));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecAssemblyEnd(*b));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscFunctionReturn(</span><span style='font-size:10.0pt;font-family:Consolas;color:#6F008A'>PETSC_SUCCESS</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> main(</span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> argc, </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>char</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>** argv)<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscFunctionBeginUser;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>const</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>char</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> help[] = </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"Poisson3D solver using PETSc."</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PetscInitialize(&argc, &argv, PETSC_NULL, help));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Mat mat_;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Vec B_;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Vec X_;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nz = 26; </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//number of z-segments</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nr = 20; </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//number of r-segments</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nt = 10; </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//number of theta-segments</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PetscOptionsGetInt(NULL, NULL, </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"-z"</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, &Nz, NULL));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PetscOptionsGetInt(NULL, NULL, </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"-r"</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, &Nr, NULL));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PetscOptionsGetInt(NULL, NULL, </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"-t"</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, &Nt, NULL));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//Neumann at z=0 and z=L</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//Dirichlet at R, the cylinder radius</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//On axis, compute field</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//number of unknowns of cylinder interior nodes minus axis</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Kz = Nz - 1;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> Ir = Nr - 1;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> M = Kz * Ir * Nt;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> bw = 15; </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//for 7-point stencil</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatCreate(PETSC_COMM_WORLD, &mat_));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetSizes(mat_, PETSC_DECIDE, PETSC_DECIDE, M, M));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetFromOptions(mat_));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatMPIAIJSetPreallocation(mat_, bw, 0, bw, 0));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSeqAIJSetPreallocation(mat_, bw, 0));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//source</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecCreate(PETSC_COMM_WORLD, &B_));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetSizes(B_, PETSC_DECIDE, M));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetFromOptions(B_));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//sol</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecDuplicate(B_, &X_));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> MPIAssembler(&mat_, &B_, Kz, Ir, Nt);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PCType pcType = PCEISENSTAT;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> KSPType solverType = KSPGMRES;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> maxIters = 100;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> relTol = 1.0e-6;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> KSP ksp;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetOperators(ksp, mat_, mat_));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetType(ksp, solverType));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//preconditioner</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PC pc;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPGetPC(ksp, &pc));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PCSetType(pc, pcType));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetTolerances(ksp, relTol, PETSC_DEFAULT, PETSC_DEFAULT, maxIters));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetFromOptions(ksp));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:green'>//solving</span><span style='font-size:10.0pt;font-family:Consolas;color:black'><o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSolve(ksp, B_, X_));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>double</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> xMin, xMax;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>int</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> mnIndx, mxIndx;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecMin(X_, &mnIndx, &xMin));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecMax(X_, &mxIndx, &xMax));<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> fprintf(stderr, </span><span style='font-size:10.0pt;font-family:Consolas;color:#A31515'>"mnIndx=%i mxIndx=%i xMin=%g xMax=%g\n"</span><span style='font-size:10.0pt;font-family:Consolas;color:black'>, mnIndx, mxIndx, xMin, xMax);<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> <o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PetscFinalize());<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><span style='font-size:10.0pt;font-family:Consolas;color:blue'>return</span><span style='font-size:10.0pt;font-family:Consolas;color:black'> 0;<o:p></o:p></span></p><p class=MsoNormal style='background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}<o:p></o:p></span></p><p class=MsoNormal><o:p> </o:p></p></div></body></html>