<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 15 (filtered medium)">
<!--[if !mso]><style>v\:* {behavior:url(#default#VML);}
o\:* {behavior:url(#default#VML);}
w\:* {behavior:url(#default#VML);}
.shape {behavior:url(#default#VML);}
</style><![endif]--><style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:"Segoe UI";
panose-1:2 11 5 2 4 2 4 2 2 3;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
span.EmailStyle19
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri",sans-serif;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:70.85pt 2.0cm 2.0cm 2.0cm;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:488833329;
mso-list-template-ids:1773971626;}
ol
{margin-bottom:0cm;}
ul
{margin-bottom:0cm;}
--></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" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal">Thank you, Matthew.<o:p></o:p></p>
<p class="MsoNormal">I’m going to pay attention to our non-dimensionalization, avoiding division by cell volume helps a lot.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Sorry, Mark, I cannot get your point: which 1D problem are you referring to? The case I’m talking about is based on a 3D octree mesh.
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thank you all for your support.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Bests,<o:p></o:p></p>
<p class="MsoNormal"><span lang="IT">Marco Cisternino</span><span lang="IT"> </span>
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0cm 0cm 0cm">
<p class="MsoNormal"><b>From:</b> Mark Adams <mfadams@lbl.gov> <br>
<b>Sent:</b> mercoledì 19 gennaio 2022 15:16<br>
<b>To:</b> Matthew Knepley <knepley@gmail.com><br>
<b>Cc:</b> Marco Cisternino <marco.cisternino@optimad.it>; petsc-users <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Nullspaces<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal">ILU is LU for this 1D problem and its singular.<o:p></o:p></p>
<div>
<p class="MsoNormal">ILU might have some logic to deal with a singular system. Not sure. LU should fail.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">You might try -ksp_type preonly and -pc_type lu <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">And -pc_type jacobi should not have any numerical problems. Try that.<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal">On Wed, Jan 19, 2022 at 7:19 AM Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm">
<div>
<div>
<p class="MsoNormal">On Wed, Jan 19, 2022 at 4:52 AM Marco Cisternino <<a href="mailto:marco.cisternino@optimad.it" target="_blank">marco.cisternino@optimad.it</a>> wrote:<o:p></o:p></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm">
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Thank you Matthew.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">But I cannot get the point. I got the point about the test but to try to explain my doubt I’m going to prepare another toy code.<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">By words…<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">I usually have a finite volume discretization of the Laplace operator with homogeneous Neumann BC on an octree mesh and it reads<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Aij * xj = bi,<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">being the discretization of<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="IT">Int|Vi(nabla^2 pi dV) = Int|Vi(nabla dot ui)</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">(I omit constants), where Int|Vi(…dV) is a volume integral on the I cell, pi is cell pressure, ui is the cell velocity.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">The computational domain contains 2 separated sub-domains.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Let’s consider 4 cells into the whole domain and 2 cells for each sub-domain.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">I would expect a null space made of 2 vectors and from your patch they should look like [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], i.e. norm2 = 1 for both.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">With MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, &nullspace) I’m saying that [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)] is the null space, which is not but it is
in the null space.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">But this is not the case I sent you, not exactly.<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">The case I sent is 1/Vi * Aij * xj = 1/Vi bi, where Vi is the volume of the cell i.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Let’s admit that yj is in the null space of of Aij, it should be in the null space of 1/Vi * Aij, therefore Aij*yj = 0 and 1/Vi * Aij*yj = 0 too.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">But in the framework of the test, this is true with infinite precision.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">What happens if norm2(Aij*yj) = 10^-15 and Vi = 10^-5?
<span lang="IT">Norm2(1/Vi * Aij * yj) = 10^-10!!! Is yi still in the null space numerically?</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Let’s say yi is the constant vector over the whole domain, i.e. [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)]. Should this be in the null space of 1/Vi * Aij, shouldn’t it?<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">An analogous argument should be for the compatibility condition that concerns bi.<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">My current problem is that properly creating the null space for Aij, i.e. [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], allows me to solve and find xi, but multiplying
by 1/Vi, I cannot get any solution using both FGMRES+ILU and FGMRE+GAMG.<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">The tiny problem will load Aij, Vi and bi and show the problem by testing the proper null space and trying to solve. I will include the patch to my PETSc version. I hope to come
back to you very soon.<o:p></o:p></p>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">This sounds like a dimensionalization problem to me. It is best to choose length (and other) units that make the matrix entries about 1. It seems<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">like you are far from this, and it is causing a loss of accuracy (your null vector has a residual of about 1e-7).<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:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm">
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Thank you very much for your support!<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">Marco Cisternino<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 style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0cm 0cm 0cm">
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><b>From:</b> Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>
<br>
<b>Sent:</b> martedì 18 gennaio 2022 21:25<br>
<b>To:</b> Marco Cisternino <<a href="mailto:marco.cisternino@optimad.it" target="_blank">marco.cisternino@optimad.it</a>><br>
<b>Cc:</b> petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b> Re: [petsc-users] Nullspaces<o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">I made a fix for this:<o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <a href="https://gitlab.com/petsc/petsc/-/merge_requests/4729" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/4729</a><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> Thanks,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> Matt<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">On Tue, Jan 18, 2022 at 3:20 PM Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0cm;margin-bottom:5.0pt">
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">On Thu, Dec 16, 2021 at 11:09 AM Marco Cisternino <<a href="mailto:marco.cisternino@optimad.it" target="_blank">marco.cisternino@optimad.it</a>> wrote:<o:p></o:p></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0cm;margin-bottom:5.0pt">
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Hello Matthew,<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">as promised I prepared a minimal (112960 rows. I’m not able to produce anything smaller than this and triggering the issue) example of the behavior I was talking about some days
ago.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">What I did is to produce matrix, right hand side and initial solution of the linear system.<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">As I told you before, this linear system is the discretization of the pressure equation of a predictor-corrector method for NS equations in the framework of finite volume method.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">This case has homogeneous Neumann boundary conditions. Computational domain has two independent and separated sub-domains.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">I discretize the weak formulation and I divide every row of the linear system by the volume of the relative cell.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">The underlying mesh is not uniform, therefore cells have different volumes.
<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">The issue I’m going to explain does not show up if the mesh is uniform, same volume for all the cells.<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 usually build the null space sub-domain by sub-domain with<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-indent:36.0pt">
MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants, constants, &nullspace);<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Where nConstants = 2 and constants contains two normalized arrays with constant values on degrees of freedom relative to the associated sub-domain and zeros elsewhere.<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">However, as a test I tried the constant over the whole domain using 2 alternatives that should produce the same null space:<o:p></o:p></p>
<ol start="1" type="1">
<li class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;mso-list:l0 level1 lfo1">
MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, &nullspace);<o:p></o:p></li><li class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;mso-list:l0 level1 lfo1">
Vec* nsp;<o:p></o:p></li></ol>
<p>VecDuplicateVecs(solution, 1, &nsp);<o:p></o:p></p>
<p>VecSet(nsp[0],1.0);<o:p></o:p></p>
<p>VecNormalize(nsp[0], nullptr);<o:p></o:p></p>
<p>MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, 1, nsp, &nullspace);<o:p></o:p></p>
<p> <o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
Once I created the null space I test it using:<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify;text-indent:36.0pt">
<span lang="IT">MatNullSpaceTest(nullspace, m_A, &isNullSpaceValid);</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
<span lang="IT"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
The case 1 pass the test while case 2 don’t.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
I have a small code for matrix loading, null spaces creation and testing.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
Unfortunately I cannot implement a small code able to produce that linear system.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
As attachment you can find an archive containing the matrix, the initial solution (used to manually build the null space) and the rhs (not used in the test code) in binary format.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
You can also find the testing code in the same archive.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
I used petsc 3.12(gcc+openMPI) and petsc 3.15.2(intelOneAPI) same results.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
If the attachment is not delivered, I can share a link to it.<o:p></o:p></p>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Marco, please forgive me for taking so long to get to your issue. It has been crazy.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">You are correct, we had a bug. it is in MatNullSpaceTest. The normalization for the constant vector was wrong:<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">diff --git a/src/mat/interface/matnull.c b/src/mat/interface/matnull.c<br>
index f8ab2925988..0c4c3855be0 100644<br>
--- a/src/mat/interface/matnull.c<br>
+++ b/src/mat/interface/matnull.c<br>
@@ -429,7 +429,7 @@ PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)<br>
if (sp->has_cnst) {<br>
ierr = VecDuplicate(l,&r);CHKERRQ(ierr);<br>
ierr = VecGetSize(l,&N);CHKERRQ(ierr);<br>
- sum = 1.0/N;<br>
<br>
+ sum = 1.0/PetscSqrtReal(N);<br>
ierr = VecSet(l,sum);CHKERRQ(ierr);<br>
ierr = MatMult(mat,l,r);CHKERRQ(ierr);<br>
ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr);<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">With this fix, your two cases give the same answer, namely that the constant vector is not a null vector of your<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">operator, but it is close, as your can see using -mat_null_space_test_view.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> Thanks,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> Matt<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0cm;margin-bottom:5.0pt">
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
Thanks for any help.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;text-align:justify">
<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="IT">Marco Cisternino
</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="IT"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="IT"> </span><o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="IT">Marco Cisternino, PhD<br>
</span><a href="mailto:marco.cisternino@optimad.it" target="_blank"><span lang="IT">marco.cisternino@optimad.it</span></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">Optimad Engineering Srl<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Via Bligny 5,
<span style="font-family:"Segoe UI",sans-serif;color:#201F1E;background:white">Torino, Italia.</span><span style="font-family:"Segoe UI",sans-serif;color:#201F1E"><br>
<span style="background:white">+3901119719782</span><br>
</span><a href="http://www.optimad.it/" target="_blank"><span style="font-family:"Segoe UI",sans-serif;border:none windowtext 1.0pt;padding:0cm;background:white">www.optimad.it</span></a><o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<div>
<div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0cm 0cm 0cm">
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><b>From:</b> Marco Cisternino <<a href="mailto:marco.cisternino@optimad.it" target="_blank">marco.cisternino@optimad.it</a>>
<br>
<b>Sent:</b> martedì 7 dicembre 2021 19:36<br>
<b>To:</b> Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
<b>Cc:</b> petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b> Re: [petsc-users] Nullspaces<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto;background:white">
<span style="color:#212121">I will, as soon as possible...</span><o:p></o:p></p>
</div>
<div id="gmail-m_2525094775288695030gmail-m_-945672450800582229gmail-m_774479102612054469gmail-m_-3783333818967129010ms-outlook-mobile-signature">
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Scarica
<a href="https://aka.ms/AAb9ysg" target="_blank">Outlook per Android</a><o:p></o:p></p>
</div>
<div class="MsoNormal" align="center" style="text-align:center">
<hr size="2" width="98%" align="center">
</div>
<div id="gmail-m_2525094775288695030gmail-m_-945672450800582229gmail-m_774479102612054469gmail-m_-3783333818967129010divRplyFwdMsg">
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><b><span style="color:black">From:</span></b><span style="color:black"> Matthew Knepley <</span><a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a><span style="color:black">><br>
<b>Sent:</b> Tuesday, December 7, 2021 7:25:43 PM<br>
<b>To:</b> Marco Cisternino <</span><a href="mailto:marco.cisternino@optimad.it" target="_blank">marco.cisternino@optimad.it</a><span style="color:black">><br>
<b>Cc:</b> petsc-users <</span><a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a><span style="color:black">><br>
<b>Subject:</b> Re: [petsc-users] Nullspaces</span> <o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">On Tue, Dec 7, 2021 at 11:19 AM Marco Cisternino <<a href="mailto:marco.cisternino@optimad.it" target="_blank">marco.cisternino@optimad.it</a>> wrote:<o:p></o:p></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0cm;margin-bottom:5.0pt">
<div>
<div>
<p>Good morning,<o:p></o:p></p>
<p><span lang="EN-GB">I’m still struggling with the Poisson equation with Neumann BCs.</span><o:p></o:p></p>
<p><span lang="EN-GB">I discretize the equation by finite volume method and I divide every line of the linear system by the volume of the cell. I could avoid this division, but I’m trying to understand.</span><o:p></o:p></p>
<p><span lang="EN-GB">My mesh is not uniform, i.e. cells have different volumes (it is an octree mesh).</span><o:p></o:p></p>
<p><span lang="EN-GB">Moreover, in my computational domain there are 2 separated sub-domains.</span><o:p></o:p></p>
<p><span lang="EN-GB">I build the null space and then I use MatNullSpaceTest to check it.</span><o:p></o:p></p>
<p><span lang="EN-GB"> </span><o:p></o:p></p>
<p><span lang="EN-GB">If I do this:</span><o:p></o:p></p>
<p><span lang="EN-GB">MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, &nullspace);</span><o:p></o:p></p>
<p><span lang="EN-GB">It works</span><o:p></o:p></p>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">This produces the normalized constant vector.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0cm;margin-bottom:5.0pt">
<div>
<div>
<p><span lang="EN-GB">If I do this:</span><o:p></o:p></p>
<p>Vec nsp;<o:p></o:p></p>
<p>VecDuplicate(m_rhs, &nsp);<o:p></o:p></p>
<p>VecSet(nsp,1.0);<o:p></o:p></p>
<p>VecNormalize(nsp, nullptr);<o:p></o:p></p>
<p><span lang="EN-GB">MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, 1, &nsp, &nullspace);</span><o:p></o:p></p>
<p><span lang="EN-GB">It does not work</span><o:p></o:p></p>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">This is also the normalized constant vector.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">So you are saying that these two vectors give different results with MatNullSpaceTest()?<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Something must be wrong in the code. Can you send a minimal example of this? I will go<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">through and debug it.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> Thanks,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> Matt<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0cm;margin-bottom:5.0pt">
<div>
<div>
<p><span lang="EN-GB">Probably, I have wrong expectations, but should not it be the same?</span><o:p></o:p></p>
<p><span lang="EN-GB"> </span><o:p></o:p></p>
<p><span lang="EN-GB">Thanks</span><o:p></o:p></p>
<p><span lang="EN-GB"> </span><o:p></o:p></p>
<p>Marco Cisternino, PhD<br>
<a href="mailto:marco.cisternino@optimad.it" target="_blank">marco.cisternino@optimad.it</a><o:p></o:p></p>
<p><span lang="EN-GB">______________________</span><o:p></o:p></p>
<p><span lang="EN-GB">Optimad Engineering Srl</span><o:p></o:p></p>
<p><span lang="IT">Via Bligny 5, </span><span lang="IT" style="font-family:"Segoe UI",sans-serif;color:#201F1E;background:white">Torino, Italia.</span><span lang="IT" style="font-family:"Segoe UI",sans-serif;color:#201F1E"><br>
<span style="background:white">+3901119719782</span><br>
</span><a href="http://www.optimad.it/" target="_blank"><span lang="IT" style="font-family:"Segoe UI",sans-serif;border:none windowtext 1.0pt;padding:0cm;background:white">www.optimad.it</span></a><o:p></o:p></p>
<p><span lang="IT"> </span><o:p></o:p></p>
</div>
</div>
</blockquote>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="IT"><br clear="all">
</span><o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="IT"> </span><o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">--
<o:p></o:p></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">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" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><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>
</div>
</div>
</blockquote>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><br clear="all">
<o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">--
<o:p></o:p></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">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" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><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>
</blockquote>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><br clear="all">
<o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">--
<o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">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" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><o:p></o:p></p>
</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">-- <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>
</blockquote>
</div>
</div>
</body>
</html>