<div dir="ltr"><div dir="ltr">On Sat, May 2, 2020 at 3:36 AM Shashwat Tiwari <<a href="mailto:shaswat121994@gmail.com">shaswat121994@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi,<br>I am writing a simple code to solve linear advection equation using a first order cell centered finite volume scheme on 2D unstructured girds using DMPlex. In order to set values to a Global vector (for example, to set the initial value of the solution vector), I am looping over the cells owned by a process (including partition ghost cells) and checking the "vtk" label for each cell, assigned by the DMPlexConstructGhostCells() function, to prevent it from writing to ghost cells. If I don't do this check, the code gives a segmentation fault, which, as far as I understand, is caused by trying to write into the ghost points which do not exist on the Global Vector.  Following is the function that I have written to set the initial condition:<br></div><div><br></div><div><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap"><div><span style="color:rgb(32,32,32)">PetscErrorCode</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">SetIC</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">DM</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">dm</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">Vec</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">U</span><span style="color:rgb(0,0,0)">)</span></div><div><span style="color:rgb(0,0,0)">{</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">PetscErrorCode</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">PetscScalar</span><span style="color:rgb(0,0,0)"> *</span><span style="color:rgb(32,32,32)">u</span><span style="color:rgb(0,0,0)">;</span></div><br><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">PetscFunctionBegin</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">PetscInt</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">c</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">cStart</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">cEnd</span><span style="color:rgb(0,0,0)">; </span><span style="color:rgb(160,160,160);font-style:italic">// cells</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">PetscReal</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">area</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">centroid</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(101,84,192)">3</span><span style="color:rgb(0,0,0)">], </span><span style="color:rgb(32,32,32)">normal</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(101,84,192)">3</span><span style="color:rgb(0,0,0)">]; </span><span style="color:rgb(160,160,160);font-style:italic">// geometric data</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(160,160,160);font-style:italic">// get cell stratum owned by processor</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">DMPlexGetHeightStratum</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">dm</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(101,84,192)">0</span><span style="color:rgb(0,0,0)">, &</span><span style="color:rgb(32,32,32)">cStart</span><span style="color:rgb(0,0,0)">, &</span><span style="color:rgb(32,32,32)">cEnd</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(160,160,160);font-style:italic">// get array for U</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">VecGetArray</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">U</span><span style="color:rgb(0,0,0)">, &</span><span style="color:rgb(32,32,32)">u</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(160,160,160);font-style:italic">// loop over cells and assign values</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(9,30,66);font-weight:bold">for</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">c</span><span style="color:rgb(0,0,0)">=</span><span style="color:rgb(32,32,32)">cStart</span><span style="color:rgb(0,0,0)">; </span><span style="color:rgb(32,32,32)">c</span><span style="color:rgb(0,0,0)"><</span><span style="color:rgb(32,32,32)">cEnd</span><span style="color:rgb(0,0,0)">; ++</span><span style="color:rgb(32,32,32)">c</span><span style="color:rgb(0,0,0)">)</span></div><div><span style="color:rgb(0,0,0)">   {</span></div><div><span style="color:rgb(0,0,0)">      </span><span style="color:rgb(32,32,32)">PetscInt</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">label</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)">      </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">DMGetLabelValue</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">dm</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(191,38,0)">"vtk"</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">c</span><span style="color:rgb(0,0,0)">, &</span><span style="color:rgb(32,32,32)">label</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">      </span><span style="color:rgb(160,160,160);font-style:italic">// write into Global vector if the cell is a real cell</span></div><div><span style="color:rgb(0,0,0)">      </span><span style="color:rgb(9,30,66);font-weight:bold">if</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">label</span><span style="color:rgb(0,0,0)"> == </span><span style="color:rgb(101,84,192)">1</span><span style="color:rgb(0,0,0)">)</span></div><div><span style="color:rgb(0,0,0)">      {</span></div><div><span style="color:rgb(0,0,0)">         </span><span style="color:rgb(32,32,32)">PetscReal</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">X</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(101,84,192)">2</span><span style="color:rgb(0,0,0)">]; </span><span style="color:rgb(160,160,160);font-style:italic">// cell centroid</span></div><div><span style="color:rgb(0,0,0)">         </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">DMPlexComputeCellGeometryFVM</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">dm</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">c</span><span style="color:rgb(0,0,0)">, &</span><span style="color:rgb(32,32,32)">area</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">centroid</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(32,32,32)">normal</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">         </span><span style="color:rgb(32,32,32)">X</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(101,84,192)">0</span><span style="color:rgb(0,0,0)">] = </span><span style="color:rgb(32,32,32)">centroid</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(101,84,192)">0</span><span style="color:rgb(0,0,0)">]; </span><span style="color:rgb(32,32,32)">X</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(101,84,192)">1</span><span style="color:rgb(0,0,0)">] = </span><span style="color:rgb(32,32,32)">centroid</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(101,84,192)">1</span><span style="color:rgb(0,0,0)">];</span></div><div><span style="color:rgb(0,0,0)">         </span><span style="color:rgb(32,32,32)">u</span><span style="color:rgb(0,0,0)">[</span><span style="color:rgb(32,32,32)">c</span><span style="color:rgb(0,0,0)">] = </span><span style="color:rgb(32,32,32)">initial_condition</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">X</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">      }</span></div><div><span style="color:rgb(0,0,0)">   }</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">VecRestoreArray</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">U</span><span style="color:rgb(0,0,0)">, &</span><span style="color:rgb(32,32,32)">u</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">   </span><span style="color:rgb(32,32,32)">PetscFunctionReturn</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(101,84,192)">0</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">}</span></div></div></div><div><br></div><div>This gives me the desired output, but, I wanted to ask if there is a better and more efficient way to achieve this, i.e. to write to a global vector, than checking the label for each cell. I am also attaching a sample code which sets the initial condition and writes a corresponding vtk file. Kindly correct me if I am wrong in my understanding and give your suggestions to improve upon this.</div></div></blockquote><div><br></div><div>1) This will work, but I think the intent is to use the "ghost" label to identify ghost cells, It has the same information as "vtk" since you do not want to output those</div><div>     cells for visualization either, but it might be more clear in the code.</div><div><br></div><div>2) If I were doing this many times, I would create a new IS with the cells I was looping over, namely [cStart, cEnd) - {ghost cells}. That way you have only a lookup rather</div><div>    than a search. However, if you do it only once or twice, I don't think there is much advantage.</div><div><br></div><div>  Thanks,</div><div><br></div><div>    Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Regards, <br></div><div>Shashwat</div></div></blockquote></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>