<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"><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><br></div><div>Regards, <br></div><div>Shashwat<br></div><div><br><br></div></div>