<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)">
<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.apple-converted-space
{mso-style-name:apple-converted-space;}
span.EmailStyle19
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:70.85pt 2.0cm 2.0cm 2.0cm;}
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="IT" link="blue" vlink="purple" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal"><span style="mso-fareast-language:EN-US">Thanks Barry.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US">I already did it, I mean a constant vec per domain with 1 on rows of that domain and zero on rows of the other one, and exactly I did this:<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> </span>
<span lang="FR" style="mso-fareast-language:EN-US">Vec *constants = nullptr;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="FR" style="mso-fareast-language:EN-US"> PetscInt nConstants = nActiveDomains;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="FR" style="mso-fareast-language:EN-US"> VecDuplicateVecs(m_solution, nConstants, &constants);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="FR" style="mso-fareast-language:EN-US"> </span>
<span lang="EN-GB" style="mso-fareast-language:EN-US">for (PetscInt i = 0; i < nConstants; ++i) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> VecSet(constants[i],0.0);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> std::vector<PetscScalar *> rawConstants(nConstants);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> for (PetscInt i = 0; i < nConstants; ++i) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> VecGetArray(constants[i], &rawConstants[i]);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> long nRows = getRowCount();<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> const CBaseLinearSolverMapping *cellRowMapping = getMapping();<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> for (long cellRow = 0; cellRow < nRows; ++cellRow) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> std::size_t cellRawId = cellRowMapping->getRowRawCell(cellRow);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> int cellDomain = m_grid->getCellStorage().getDomains().rawAt(cellRawId);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> if (m_grid->isCellDomainActive(cellDomain)) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> std::size_t activeDomainId = activeDomainIds[cellDomain];<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> rawConstants[activeDomainId][cellRow] = 1.;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> for (PetscInt i = 0; i < nConstants; ++i) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> VecRestoreArray(constants[i], &rawConstants[i]);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> VecAssemblyBegin(constants[i]);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> VecAssemblyEnd(constants[i]);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> VecNormalize(constants[i], nullptr);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"> MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants, constants, &nullspace);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US">But this does not pass the test if I divide by cell volume every row of the linear system.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US">It works if I do not perform the division<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US">Thanks<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<div>
<p class="MsoNormal">Marco Cisternino, PhD<br>
<a href="mailto:marco.cisternino@optimad.it">marco.cisternino@optimad.it</a><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-GB">______________________<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB">Optimad Engineering Srl<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-GB">Via Bligny 5, </span><span lang="EN-GB" style="font-family:"Segoe UI",sans-serif;color:#201F1E;background:white">Torino, Italia.</span><span lang="EN-GB" 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="EN-GB" style="font-family:"Segoe UI",sans-serif;border:none windowtext 1.0pt;padding:0cm;background:white">www.optimad.it</span></a><span lang="EN-GB"><o:p></o:p></span></p>
</div>
<p class="MsoNormal"><span lang="EN-GB" style="mso-fareast-language:EN-US"><o:p> </o:p></span></p>
<div>
<div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0cm 0cm 0cm">
<p class="MsoNormal"><b><span lang="EN-US">From:</span></b><span lang="EN-US"> Barry Smith <bsmith@petsc.dev>
<br>
<b>Sent:</b> martedì 7 dicembre 2021 17:53<br>
<b>To:</b> Marco Cisternino <marco.cisternino@optimad.it><br>
<b>Cc:</b> petsc-users <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Nullspaces<o:p></o:p></span></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<p class="MsoNormal"> A side note: The MatNullSpaceTest tells you that the null space you provided is in the null space of the operator, it does not say if you have found the entire null space. In your case with two subdomains the null space is actually two
dimensional; all constant on one domain (0 on the other) and 0 on the first domain and all constant on the second. So you need to pass two vectors into the MatNullSpaceCreate(). But regardless the test should work for all constant on both domains.<o:p></o:p></p>
<div>
<p class="MsoNormal"><br>
<br>
<o:p></o:p></p>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<p class="MsoNormal">On Dec 7, 2021, at 11:19 AM, Marco Cisternino <<a href="mailto:marco.cisternino@optimad.it">marco.cisternino@optimad.it</a>> wrote:<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal">Good morning,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">I’m still struggling with the Poisson equation with Neumann BCs.</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><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>
</div>
<div>
<p class="MsoNormal"><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>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">Moreover, in my computational domain there are 2 separated sub-domains.</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">I build the null space and then I use MatNullSpaceTest to check it.</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">If I do this:</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, &nullspace);</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">It works</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">If I do this:</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="FR">Vec nsp;</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="FR">VecDuplicate(m_rhs, &nsp);</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="FR">VecSet(nsp,1.0);</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="FR">VecNormalize(nsp, nullptr);</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, 1, &nsp, &nullspace);</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">It does not work</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">Probably, I have wrong expectations, but should not it be the same?</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">Thanks</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">Marco Cisternino, PhD<br>
<a href="mailto:marco.cisternino@optimad.it">marco.cisternino@optimad.it</a><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">______________________</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">Optimad Engineering Srl</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span lang="EN-GB">Via Bligny 5,<span class="apple-converted-space"> </span></span><span lang="EN-GB" style="font-family:"Segoe UI",sans-serif;color:#201F1E;background:white">Torino, Italia.</span><span lang="EN-GB" 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="EN-GB" 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>
</div>
</blockquote>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</body>
</html>