<div dir="ltr"><div dir="ltr">On Mon, Jul 24, 2023 at 8:16 PM Thuc Bui <<a href="mailto:bui@calcreek.com">bui@calcreek.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg-7688111216358332865"><div lang="EN-US"><div class="m_-7688111216358332865WordSection1"><p class="MsoNormal">Dear PETSc Users/Developers.<u></u><u></u></p><p class="MsoNormal"><u></u> <u></u></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.<u></u><u></u></p><p class="MsoNormal"><u></u> <u></u></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.<u></u><u></u></p><p class="MsoNormal"><u></u> <u></u></p><p class="MsoNormal">I would be very grateful if someone can provide me with some direction to fix these issues.</p></div></div></div></blockquote><div><br></div><div>That is a bug with PCEISENSTAT. It is not frequently used. Since your matrix does not use inodes, a default implementation of MatMultDIagonalBlock is not found. We just need to add this for MATSEQAIJ.</div><div>It should work with a different PC. We will fix this as soon as we can.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg-7688111216358332865"><div lang="EN-US"><div class="m_-7688111216358332865WordSection1"><p class="MsoNormal">Many thanks for your help,<u></u><u></u></p><p class="MsoNormal">Thuc Bui<u></u><u></u></p><p class="MsoNormal">Senior R&D Engineer<u></u><u></u></p><p class="MsoNormal">Calabazas Creek Research, Inc<u></u><u></u></p><p class="MsoNormal">(650) 948-5361<u></u><u></u></p><p class="MsoNormal"><u></u> <u></u></p><p class="MsoNormal"><u></u> <u></u></p><p class="MsoNormal">///////////////errors/////////////////////////<u></u><u></u></p><p class="MsoNormal">rank=0 iStart=0 iEnd=2375<u></u><u></u></p><p class="MsoNormal">rank=1 iStart=2375 iEnd=4750<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">--------------------- Error Message --------------------------------------------------------------<u></u><u></u></p><p class="MsoNormal">--------------------- Error Message --------------------------------------------------------------<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">No support for this operation for this object type<u></u><u></u></p><p class="MsoNormal">No support for this operation for this object type<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">No method multdiagonalblock for Mat of type seqaij<u></u><u></u></p><p class="MsoNormal">No method multdiagonalblock for Mat of type seqaij<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">See <a href="https://petsc.org/release/faq/" target="_blank">https://petsc.org/release/faq/</a> for trouble shooting.<u></u><u></u></p><p class="MsoNormal">See <a href="https://petsc.org/release/faq/" target="_blank">https://petsc.org/release/faq/</a> for trouble shooting.<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">Petsc Release Version 3.18.6, Mar 30, 2023 <u></u><u></u></p><p class="MsoNormal">Petsc Release Version 3.18.6, Mar 30, 2023 <u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">Poisson3D.exe on a v140-x64SReleaseSP named HERMES by bui Mon Jul 24 15:16:55 2023<u></u><u></u></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<u></u><u></u></p><p class="MsoNormal">Poisson3D.exe on a v140-x64SReleaseSP named HERMES by bui Mon Jul 24 15:16:55 2023<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#1 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<u></u><u></u></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<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#2 MatMultDiagonalBlock_MPIAIJ() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\impls\aij\mpi\mpiaij.c:978<u></u><u></u></p><p class="MsoNormal">#1 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#3 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<u></u><u></u></p><p class="MsoNormal">#2 MatMultDiagonalBlock_MPIAIJ() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\impls\aij\mpi\mpiaij.c:978<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#4 PCApply_Eisenstat() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\impls\eisens\eisen.c:51<u></u><u></u></p><p class="MsoNormal">#3 MatMultDiagonalBlock() at C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#5 PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\INTERF~1\precon.c:441<u></u><u></u></p><p class="MsoNormal">#4 PCApply_Eisenstat() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\impls\eisens\eisen.c:51<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#6 KSP_PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\include\petsc/private/kspimpl.h:380<u></u><u></u></p><p class="MsoNormal">#5 PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\INTERF~1\precon.c:441<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#7 KSPInitialResidual() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itres.c:64<u></u><u></u></p><p class="MsoNormal">#6 KSP_PCApply() at C:\Users\bui\DOCUME~1\Petsc\latest\include\petsc/private/kspimpl.h:380<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#8 KSPSolve_GMRES() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\impls\gmres\gmres.c:227<u></u><u></u></p><p class="MsoNormal">#7 KSPInitialResidual() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itres.c:64<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#9 KSPSolve_Private() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:899<u></u><u></u></p><p class="MsoNormal">#8 KSPSolve_GMRES() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\impls\gmres\gmres.c:227<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#10 KSPSolve() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:1071<u></u><u></u></p><p class="MsoNormal">#9 KSPSolve_Private() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:899<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">#11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253<u></u><u></u></p><p class="MsoNormal">#10 KSPSolve() at C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:1071<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">No PETSc Option Table entries<u></u><u></u></p><p class="MsoNormal">#11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253<u></u><u></u></p><p class="MsoNormal">[1]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: <u></u><u></u></p><p class="MsoNormal">----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov----------<u></u><u></u></p><p class="MsoNormal">No PETSc Option Table entries<u></u><u></u></p><p class="MsoNormal">[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to <a href="mailto:petsc-maint@mcs.anl.gov----------" target="_blank">petsc-maint@mcs.anl.gov----------</a><u></u><u></u></p><p class="MsoNormal"><u></u> <u></u></p><p class="MsoNormal">//complete codes///////////////////////////////////////////////<u></u><u></u></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:gray">#include</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"petscsys.h"</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:gray">#include</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"petscksp.h"</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:gray">#ifndef</span><span style="font-size:10pt;font-family:Consolas;color:black"> PETSC_SUCCESS<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:gray">#define</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:rgb(111,0,138)">PETSC_SUCCESS</span><span style="font-size:10pt;font-family:Consolas;color:black"> 0<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:gray">#endif</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:green">//rectangular matrix mxn</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black">** allocateIMatrixR(</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:gray">m</span><span style="font-size:10pt;font-family:Consolas;color:black">, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:gray">n</span><span style="font-size:10pt;font-family:Consolas;color:black">)<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">{<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black">** v = (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black">**)malloc(</span><span style="font-size:10pt;font-family:Consolas;color:gray">m</span><span style="font-size:10pt;font-family:Consolas;color:black"> * </span><span style="font-size:10pt;font-family:Consolas;color:blue">sizeof</span><span style="font-size:10pt;font-family:Consolas;color:black">(</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black">*));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> i = 0; i < </span><span style="font-size:10pt;font-family:Consolas;color:gray">m</span><span style="font-size:10pt;font-family:Consolas;color:black">; ++i) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 v[i] = (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black">*)malloc(</span><span style="font-size:10pt;font-family:Consolas;color:gray">n</span><span style="font-size:10pt;font-family:Consolas;color:black"> * </span><span style="font-size:10pt;font-family:Consolas;color:blue">sizeof</span><span style="font-size:10pt;font-family:Consolas;color:black">(</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black">));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> i = 0; i < </span><span style="font-size:10pt;font-family:Consolas;color:gray">m</span><span style="font-size:10pt;font-family:Consolas;color:black">; ++i) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> j = 0; j < </span><span style="font-size:10pt;font-family:Consolas;color:gray">n</span><span style="font-size:10pt;font-family:Consolas;color:black">; ++</span><span style="font-size:10pt;font-family:Consolas;color:gray">n</span><span style="font-size:10pt;font-family:Consolas;color:black">) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          v[i][j] = -1;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">return</span><span style="font-size:10pt;font-family:Consolas;color:black"> v;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">}<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> mapIndex(</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:gray">k</span><span style="font-size:10pt;font-family:Consolas;color:black">, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:gray">i</span><span style="font-size:10pt;font-family:Consolas;color:black">, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:gray">j</span><span style="font-size:10pt;font-family:Consolas;color:black">, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:gray">Nz</span><span style="font-size:10pt;font-family:Consolas;color:black">, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:gray">Nr</span><span style="font-size:10pt;font-family:Consolas;color:black">)<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">{<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//k for z, i for r, j for theta</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//base 0</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//3D table:      theta    r                z</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        0                0                 0</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        0                0                 1</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        0                0                 2</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        ...</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        0                1                 0</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        0                1                 1</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        0                1                 2</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//                        ...</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> n = </span><span style="font-size:10pt;font-family:Consolas;color:gray">k</span><span style="font-size:10pt;font-family:Consolas;color:black"> + </span><span style="font-size:10pt;font-family:Consolas;color:gray">i</span><span style="font-size:10pt;font-family:Consolas;color:black"> * </span><span style="font-size:10pt;font-family:Consolas;color:gray">Nz</span><span style="font-size:10pt;font-family:Consolas;color:black"> + </span><span style="font-size:10pt;font-family:Consolas;color:gray">j</span><span style="font-size:10pt;font-family:Consolas;color:black"> * </span><span style="font-size:10pt;font-family:Consolas;color:gray">Nz</span><span style="font-size:10pt;font-family:Consolas;color:black"> * </span><span style="font-size:10pt;font-family:Consolas;color:gray">Nr</span><span style="font-size:10pt;font-family:Consolas;color:black">;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">return</span><span style="font-size:10pt;font-family:Consolas;color:black"> n;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">}<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">PetscErrorCode MPIAssembler(Mat* A, Vec* b, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Kz, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Ir, </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Jt)<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">{<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscFunctionBeginUser;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> M;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecGetSize(*b, &M));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;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:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//i is r-index, j theta-index, k z-index</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">const</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> dim = 3;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black">** Nijk = allocateIMatrixR(M, dim);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> k = 0; k < Kz; ++k) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> k1 = k + 1;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> i = 0; i < Ir; ++i) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> j = 0; j < Jt; ++j) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                                   </span><span style="font-size:10pt;font-family:Consolas;color:green">//for each equation row m</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                                   </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> m = mapIndex(k, i, j, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                                   Nijk[m][0] = i;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                                   Nijk[m][1] = j;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                                   Nijk[m][2] = k;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">/////////////////////////////////////////////////////////////////////////////////</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//Cylinder dimensions</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> R = 11.0; </span><span style="font-size:10pt;font-family:Consolas;color:green">//cm, cylinder radius</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> L = 11.0; </span><span style="font-size:10pt;font-family:Consolas;color:green">//cm, cylinder length in z</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Dr = R / (Ir + 1);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Dz = L / (Kz + 1);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Dr2 = Dr * Dr;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Dr_2 = 0.5 * Dr;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Dr_22 = Dr_2 * Dr_2;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Dz_2 = 1.0 / (Dz * Dz);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> DT = 2.0 * M_PI / Jt;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> DT2 = DT * DT;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Vr = 0.0; </span><span style="font-size:10pt;font-family:Consolas;color:green">//Dirichlet on R</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> chargeDens = -4.0;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> rank;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//get the range of matrix rows owned by this processor</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> startingNumber, endingNumber;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatGetOwnershipRange(*A, &startingNumber, &endingNumber));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         fprintf(stderr, </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"rank=%i iStart=%i iEnd=%i\n"</span><span style="font-size:10pt;font-family:Consolas;color:black">, rank, startingNumber, endingNumber);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Ir_1 = Ir - 1;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Jt_1 = Jt - 1;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> m = startingNumber; m < endingNumber; ++m) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> i = Nijk[m][0];<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> j = Nijk[m][1];<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> k = Nijk[m][2];<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> ri = Dr * (i + 1);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> ri_2 = ri - Dr_2;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> ri2 = ri + Dr_2;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> sqri = ri * ri;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Ai = (ri_2 / ri) / Dr2;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Bi = (ri2 / ri) / Dr2;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Ci = 1.0 / (sqri * DT2);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> Ei = (ri2 + ri_2) / (ri * Dr2) + 2.0 / (sqri * DT2) + 2.0 * Dz_2;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">//contribution of charge density to RHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 PetscCall(VecSetValue(*b, m, -chargeDens, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">//contribution to RHS by Drichlet on R</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">if</span><span style="font-size:10pt;font-family:Consolas;color:black"> (i == Ir_1) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          PetscCall(VecSetValue(*b, i, -Bi * Vr, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">if</span><span style="font-size:10pt;font-family:Consolas;color:black"> (!i) {        </span><span style="font-size:10pt;font-family:Consolas;color:green">//from the axis</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:green">//contributing to RHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          PetscCall(VecSetValue(*b, m, -chargeDens * Ai * Dr_22, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:green">//contributing to LHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> xTheta = (Ai / Jt);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> jt = 0; jt < Jt; ++jt) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                                   </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> nj = mapIndex(k, i, jt, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                                   PetscCall(MatSetValues(*A, 1, &m, 1, &nj, &xTheta, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">else</span><span style="font-size:10pt;font-family:Consolas;color:black"> {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:green">//0 < i < Ir: contribution of Ai to LHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> i_jk = mapIndex(k, i - 1, j, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          PetscCall(MatSetValues(*A, 1, &m, 1, &i_jk, &Ai, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">if</span><span style="font-size:10pt;font-family:Consolas;color:black"> (i != Ir_1) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:green">//0 <= i < It-1: contribution of Bi to LHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> ki1j = mapIndex(k, i + 1, j, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          PetscCall(MatSetValues(*A, 1, &m, 1, &ki1j, &Bi, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">//apply periodic BC in theta////////////////////////////</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> kij_;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">if</span><span style="font-size:10pt;font-family:Consolas;color:black"> (!j) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          kij_ = mapIndex(k, i, Jt_1, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">else</span><span style="font-size:10pt;font-family:Consolas;color:black"> {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          kij_ = mapIndex(k, i, j - 1, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">//1st contribution of Ci to LHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                  PetscCall(MatSetValues(*A, 1, &m, 1, &kij_, &Ci, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> kij1;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">if</span><span style="font-size:10pt;font-family:Consolas;color:black"> (j == Jt_1) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          kij1 = mapIndex(k, i, 0, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">else</span><span style="font-size:10pt;font-family:Consolas;color:black"> {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          kij1 = mapIndex(k, i, j + 1, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">//2nd contribution of Ci to LHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 PetscCall(MatSetValues(*A, 1, &m, 1, &kij1, &Ci, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">////////////////////////////////////////////////</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">//apply z Neumann BC ///////////////////////////////////</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> k_ij;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">if</span><span style="font-size:10pt;font-family:Consolas;color:black"> (!k) {        </span><span style="font-size:10pt;font-family:Consolas;color:green">//z=0</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          k_ij = mapIndex(k, i, j, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">else</span><span style="font-size:10pt;font-family:Consolas;color:black"> {           </span><span style="font-size:10pt;font-family:Consolas;color:green">//0 < z < L</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          k_ij = mapIndex(k - 1, i, j, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 PetscCall(MatSetValues(*A, 1, &m, 1, &k_ij, &Dz_2, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">/////////////////////////////////////////////////////////</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:green">//contribution of Ei to LHS</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> kij = mapIndex(k, i, j, Kz, Ir);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 PetscCall(MatSetValues(*A, 1, &m, 1, &kij, &Ei, ADD_VALUES));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">if</span><span style="font-size:10pt;font-family:Consolas;color:black"> (Nijk) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 </span><span style="font-size:10pt;font-family:Consolas;color:blue">for</span><span style="font-size:10pt;font-family:Consolas;color:black"> (</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> i = 0; i < M; ++i) {<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                          free(Nijk[i]);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 free(Nijk);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         }<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecAssemblyBegin(*b));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecAssemblyEnd(*b));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscFunctionReturn(</span><span style="font-size:10pt;font-family:Consolas;color:rgb(111,0,138)">PETSC_SUCCESS</span><span style="font-size:10pt;font-family:Consolas;color:black">);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">}<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> main(</span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> argc, </span><span style="font-size:10pt;font-family:Consolas;color:blue">char</span><span style="font-size:10pt;font-family:Consolas;color:black">** argv)<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">{<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscFunctionBeginUser;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">const</span><span style="font-size:10pt;font-family:Consolas;color:black"> </span><span style="font-size:10pt;font-family:Consolas;color:blue">char</span><span style="font-size:10pt;font-family:Consolas;color:black"> help[] = </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"Poisson3D solver using PETSc."</span><span style="font-size:10pt;font-family:Consolas;color:black">;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(PetscInitialize(&argc, &argv, PETSC_NULL, help));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         Mat mat_;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         Vec B_;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         Vec X_;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Nz = 26;     </span><span style="font-size:10pt;font-family:Consolas;color:green">//number of z-segments</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Nr = 20;     </span><span style="font-size:10pt;font-family:Consolas;color:green">//number of r-segments</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Nt = 10;     </span><span style="font-size:10pt;font-family:Consolas;color:green">//number of theta-segments</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(PetscOptionsGetInt(NULL, NULL, </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"-z"</span><span style="font-size:10pt;font-family:Consolas;color:black">, &Nz, NULL));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(PetscOptionsGetInt(NULL, NULL, </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"-r"</span><span style="font-size:10pt;font-family:Consolas;color:black">, &Nr, NULL));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(PetscOptionsGetInt(NULL, NULL, </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"-t"</span><span style="font-size:10pt;font-family:Consolas;color:black">, &Nt, NULL));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//Neumann at z=0 and z=L</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//Dirichlet at R, the cylinder radius</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//On axis, compute field</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//number of unknowns of cylinder interior nodes minus axis</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Kz = Nz - 1;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> Ir = Nr - 1;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> M = Kz * Ir * Nt;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> bw = 15;     </span><span style="font-size:10pt;font-family:Consolas;color:green">//for 7-point stencil</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatCreate(PETSC_COMM_WORLD, &mat_));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatSetSizes(mat_, PETSC_DECIDE, PETSC_DECIDE, M, M));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatSetFromOptions(mat_));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatMPIAIJSetPreallocation(mat_, bw, 0, bw, 0));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(MatSeqAIJSetPreallocation(mat_, bw, 0));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//source</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecCreate(PETSC_COMM_WORLD, &B_));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecSetSizes(B_, PETSC_DECIDE, M));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecSetFromOptions(B_));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//sol</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecDuplicate(B_, &X_));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         MPIAssembler(&mat_, &B_, Kz, Ir, Nt);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PCType pcType = PCEISENSTAT;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         KSPType solverType = KSPGMRES;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> maxIters = 100;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> relTol = 1.0e-6;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         KSP ksp;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(KSPSetOperators(ksp, mat_, mat_));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(KSPSetType(ksp, solverType));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//preconditioner</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PC pc;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(KSPGetPC(ksp, &pc));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(PCSetType(pc, pcType));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(KSPSetTolerances(ksp, relTol, PETSC_DEFAULT, PETSC_DEFAULT, maxIters));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(KSPSetFromOptions(ksp));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:green">//solving</span><span style="font-size:10pt;font-family:Consolas;color:black"><u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(KSPSolve(ksp, B_, X_));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black"> <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">double</span><span style="font-size:10pt;font-family:Consolas;color:black"> xMin, xMax;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">int</span><span style="font-size:10pt;font-family:Consolas;color:black"> mnIndx, mxIndx;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecMin(X_, &mnIndx, &xMin));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(VecMax(X_, &mxIndx, &xMax));<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         fprintf(stderr, </span><span style="font-size:10pt;font-family:Consolas;color:rgb(163,21,21)">"mnIndx=%i mxIndx=%i xMin=%g xMax=%g\n"</span><span style="font-size:10pt;font-family:Consolas;color:black">, mnIndx, mxIndx, xMin, xMax);<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">                 <u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         PetscCall(PetscFinalize());<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">         </span><span style="font-size:10pt;font-family:Consolas;color:blue">return</span><span style="font-size:10pt;font-family:Consolas;color:black"> 0;<u></u><u></u></span></p><p class="MsoNormal" style="background:white"><span style="font-size:10pt;font-family:Consolas;color:black">}<u></u><u></u></span></p><p class="MsoNormal"><u></u> <u></u></p></div></div></div></blockquote></div><br clear="all"><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>