Issue with PETSc vector

[Giselle Sosa Jones]

Hello,

I have an issue that seems very simple but I cannot figure out. I have a piecewise constant parallel vector that I defined using a DG-0 space. I have an auxiliary function that calculates the value according to the cell, and then in my main file, I would like to assign the value to the vector. The values are calculated correctly, but it seems like assigning to the vector is not. The values should either be 1.e-8 or 1.e-11, and the vector is 0 sometimes. I noticed this because I am printing to a vtu file and plotting. A sample code is below. Instead of using the index of the cell, I also tried using the dof index, but I am getting similar results.

Thanks in advance for the help,

Giselle

PETScWrappers::MPI::Vector vec;

//

// Initialize fe space

//

// Initialize vector

dof_handler_dg0.distribute_dofs(fe_dg0);
const std::vector<IndexSet> locally_owned_dofs_per_proc_dg0 =
DoFTools::locally_owned_dofs_per_subdomain(dof_handler_dg0);
locally_owned_dofs_dg0 = locally_owned_dofs_per_proc_dg0[this_mpi_process];

vec.reinit(locally_owned_dofs_dg0, mpi_communicator);

// Calculate values

for (const auto &cell : dof_handler_dg0.active_cell_iterators())
{
double val;
if(cell->subdomain_id() == this_mpi_process)
{
val = compute_value<dim>(cell);

                vec[cell->active_cell_index()] = val;
}

       vec.compress(VectorOperation::insert);

}

[Wolfgang Bangerth]

On 12/18/23 13:23, Giselle Sosa Jones wrote:
> // Initialize vector
> dof_handler_dg0.distribute_dofs(fe_dg0);
> const std::vector<IndexSet> locally_owned_dofs_per_proc_dg0 =
> DoFTools::locally_owned_dofs_per_subdomain(dof_handler_dg0);
> locally_owned_dofs_dg0 = locally_owned_dofs_per_proc_dg0[this_mpi_process];
> vec.reinit(locally_owned_dofs_dg0, mpi_communicator);
>
> // Calculate values
> for (const auto &cell : dof_handler_dg0.active_cell_iterators())
> {
> double val;
> if(cell->subdomain_id() == this_mpi_process)
> {
> val = compute_value<dim>(cell);
>                 vec[cell->active_cell_index()] = val;
> }
>        vec.compress(VectorOperation::insert);
> }

Gisella:
There are potentially three problems with this code, though I don't know
how many of these are actual ones :-)

1/ You treat the 'vec' vector as one associated with degrees of freedom,
but you index into it with an active_cell_index. This isn't right. You
need to index into it with the DoF index on this cell. I'm pretty sure
this is already going to fix your issue.

2/ You want to do the vec.compress() call only once after assembly, not
once per cell. In fact, I'm surprised this works the way you have it:
compress() is a *collective* cell in which all processes must
participate (or you will get a deadlock), but different processes will
have different numbers of active cells, so they will call this function
a different number of times.

3/ DataOut::add_data_vector() has a discussion that explains the 'type'
argument, see
https://dealii.org/developer/doxygen/deal.II/classDataOut__DoFData.html#ac43d3b1f6e67424f36e474627fe8e401
You may want to read through that and perhaps provide an explicit value
for the third argument (if you aren't already).

Best
W.

[Giselle Sosa Jones]

Thank you so much! I did all of this and it works now. I had tried the dof index before, but I think the output to a vtu file was wrong, so now I have fixed that too.

Much appreciated,

Giselle