<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=utf-8"><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:Tahoma;
panose-1:2 11 6 4 3 5 4 4 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:12.0pt;
font-family:"Times New Roman","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.gmailsignatureprefix
{mso-style-name:gmail_signature_prefix;}
span.EmailStyle18
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.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><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Thank you so much Matt, for getting back to me so quickly. Yes, using another PC fixes the issues.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Many thanks again for your help,<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Thuc<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'>From:</span></b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> Matthew Knepley [mailto:knepley@gmail.com] <br><b>Sent:</b> Monday, July 24, 2023 6:50 PM<br><b>To:</b> Thuc Bui<br><b>Cc:</b> petsc-users<br><b>Subject:</b> Re: [petsc-users] 3D Poisson solver failed in KSPSolve when number of cores is larger than one<o:p></o:p></span></p><p class=MsoNormal><o:p> </o:p></p><div><div><p class=MsoNormal>On Mon, Jul 24, 2023 at 8:16 PM Thuc Bui <<a href="mailto:bui@calcreek.com">bui@calcreek.com</a>> wrote:<o:p></o:p></p></div><div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in'><div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Dear PETSc Users/Developers.<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>I would be very grateful if someone can provide me with some direction to fix these issues.<o:p></o:p></p></div></div></div></blockquote><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>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.<o:p></o:p></p></div><div><p class=MsoNormal>It should work with a different PC. We will fix this as soon as we can.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> Thanks,<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> Matt<o:p></o:p></p></div><div><p class=MsoNormal> <o:p></o:p></p></div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in'><div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Many thanks for your help,<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Thuc Bui<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Senior R&D Engineer<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Calabazas Creek Research, Inc<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>(650) 948-5361<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>///////////////errors/////////////////////////<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>rank=0 iStart=0 iEnd=2375<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>rank=1 iStart=2375 iEnd=4750<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>--------------------- Error Message --------------------------------------------------------------<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>--------------------- Error Message --------------------------------------------------------------<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>No support for this operation for this object type<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>No support for this operation for this object type<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>No method multdiagonalblock for Mat of type seqaij<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>No method multdiagonalblock for Mat of type seqaij<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>See <a href="https://petsc.org/release/faq/" target="_blank">https://petsc.org/release/faq/</a> for trouble shooting.<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>See <a href="https://petsc.org/release/faq/" target="_blank">https://petsc.org/release/faq/</a> for trouble shooting.<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Petsc Release Version 3.18.6, Mar 30, 2023 <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Petsc Release Version 3.18.6, Mar 30, 2023 <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>No PETSc Option Table entries<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>#11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[1]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[0]PETSC ERROR: <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov----------<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>No PETSc Option Table entries<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>[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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>//complete codes///////////////////////////////////////////////<o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:gray'>#endif</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:green'>//rectangular matrix mxn</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>)</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>*));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> v[i][j] = -1;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>)</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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)</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscFunctionBeginUser;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecGetSize(*b, &M));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nijk[m][0] = i;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nijk[m][1] = j;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Nijk[m][2] = k;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatGetOwnershipRange(*A, &startingNumber, &endingNumber));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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];</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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];</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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];</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetValue(*b, m, -chargeDens, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetValue(*b, i, -Bi * Vr, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetValue(*b, m, -chargeDens * Ai * Dr_22, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &nj, &xTheta, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &i_jk, &Ai, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &ki1j, &Bi, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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_;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij_ = mapIndex(k, i, Jt_1, Kz, Ir);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij_ = mapIndex(k, i, j - 1, Kz, Ir);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &kij_, &Ci, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij1 = mapIndex(k, i, 0, Kz, Ir);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> kij1 = mapIndex(k, i, j + 1, Kz, Ir);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &kij1, &Ci, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> k_ij = mapIndex(k, i, j, Kz, Ir);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> k_ij = mapIndex(k - 1, i, j, Kz, Ir);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetValues(*A, 1, &m, 1, &kij, &Ei, ADD_VALUES));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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) {</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> free(Nijk[i]);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> free(Nijk);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> }</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecAssemblyBegin(*b));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecAssemblyEnd(*b));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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)</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>{</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscFunctionBeginUser;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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'>;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PetscInitialize(&argc, &argv, PETSC_NULL, help));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Mat mat_;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Vec B_;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> Vec X_;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatCreate(PETSC_COMM_WORLD, &mat_));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetSizes(mat_, PETSC_DECIDE, PETSC_DECIDE, M, M));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSetFromOptions(mat_));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatMPIAIJSetPreallocation(mat_, bw, 0, bw, 0));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(MatSeqAIJSetPreallocation(mat_, bw, 0));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecCreate(PETSC_COMM_WORLD, &B_));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetSizes(B_, PETSC_DECIDE, M));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecSetFromOptions(B_));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecDuplicate(B_, &X_));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> MPIAssembler(&mat_, &B_, Kz, Ir, Nt);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PCType pcType = PCEISENSTAT;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> KSPType solverType = KSPGMRES;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> KSP ksp;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetOperators(ksp, mat_, mat_));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetType(ksp, solverType));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PC pc;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPGetPC(ksp, &pc));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PCSetType(pc, pcType));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetTolerances(ksp, relTol, PETSC_DEFAULT, PETSC_DEFAULT, maxIters));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSetFromOptions(ksp));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(KSPSolve(ksp, B_, X_));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecMin(X_, &mnIndx, &xMin));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(VecMax(X_, &mxIndx, &xMax));</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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);</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> </span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'> PetscCall(PetscFinalize());</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;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;</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white'><span style='font-size:10.0pt;font-family:Consolas;color:black'>}</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div></div></div></blockquote></div><p class=MsoNormal><br clear=all><o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p></div><p class=MsoNormal><span class=gmailsignatureprefix>-- </span><o:p></o:p></p><div><div><div><div><div><div><div><p class=MsoNormal>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<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><o:p></o:p></p></div></div></div></div></div></div></div></div></div></body></html>