How to write values on elements in parallel?

[bzjan]

Instead of writing the value of a function on the nodes, I would like to write it on the (triangular) elements. The idea is to for each element write the desired value on its (3) constituent nodes and then later divide the value on the node by the node degree.

// write 1's to vector
for(const auto & elem : system.get_mesh().active_local_element_ptr_range()){
	system.solution->add(elem->node_ptr(0)->dof_number(system.number(),0,0),1);
	system.solution->add(elem->node_ptr(1)->dof_number(system.number(),0,0),1);
	system.solution->add(elem->node_ptr(2)->dof_number(system.number(),0,0),1);
}
system.solution->close();

// create node degrees = averaging weights
avg_weights.resize(mesh.n_nodes());
for(const auto & elem : mesh.active_local_element_ptr_range()){
	for(const auto & i : elem->node_index_range()) avg_weights[ elem->node_ptr(i)->id() ]++;
}
mesh.comm().barrier();
mesh.comm().sum(avg_weights);

// weight values
for(const auto & idx : libMesh::index_range(*system.solution)){
	libMesh::Real avg_weight = avg_weights[ idx/system.n_vars() ];
	if(avg_weight>0) system.solution->set(idx, system.solution->el(idx)/avg_weight);
}
system.solution->close();

This works very well in serial (function is = 1 everywhere):

but fails horribly in parallel (2 processors):

The miscounts occur along the boundary of the 2 processor domains.

After comparing the weight values in serial and parallel mode, it turns out that they are identical, so it seems that the problem is in the solution vector. Does close() not add up the the components of the vector across processors?

What am I missing to get this to work when running in parallel?

[roystgnr]

The solution vector is of type PARALLEL, and can only return components owned by the processor you're querying it from. Even current_local_solution is of type GHOSTED and will only return semilocal components (typically those that are either on local elements or their neighbors. Normally we'd scream and die if you tried to ask for anything else, though. I'm guessing you're running in opt mode and we end up returning 0 or returning uninitialized data in that case? The first answer to "What do I do if I see a bug in opt mode" is always "Run it in devel or dbg mode instead".

I'm not sure I understand the motivation behind what you're doing, so it feels like you're asking an XY Problem and I'm only answering Y, but your proximal bug might be fixed as simply as doing an update() to current_local_solution, reading from that, and iterating over local elements rather than over global solution indices.

[jwpeterson]

for(const auto & idx : libMesh::index_range(*system.solution)){
	libMesh::Real avg_weight = avg_weights[ idx/system.n_vars() ];
	if(avg_weight>0) system.solution->set(idx, system.solution->el(idx)/avg_weight);
}

Looping over libMesh::index_range(*system.solution) will only hit the local indices of the parallel vector, so that part looks OK to me... you won't access values that you don't have in the solution vector. I'm a bit concerned about the indexing step:

libMesh::Real avg_weight = avg_weights[ idx/system.n_vars() ];

It seems that you assume a certain ordering of the DOFs here which might not be correct.

[roystgnr]

I forgot about that bit of cleverness in index_range; please ignore my suggested fix.

John's almost certainly right about the real problem; the way we order DoFs in parallel depends on partitioning and certainly is never guaranteed to match the way we order nodes.

[bzjan]

Thanks for clearing the index issue up. I assumed that local index and node index were the same thing.

To viusalize that, I plotted them for serial and parallel code:

Using node indices explicitly seems to work well. I now have:

// create node degrees = averaging weights
for(const auto & elem : mesh.active_local_element_ptr_range()){
	for(const auto & i : elem->node_index_range()) system.solution->add(elem->node_ptr(i)->dof_number(system.number(),1,0),1.0);
}

// averaging
for(const auto & local_node_ptr : mesh.local_node_ptr_range()){
	libMesh::Real avg_weight = system.solution->el(local_node_ptr->dof_number(system.number(),1,0));
	if(avg_weight>0) system.solution->set(
	local_node_ptr->dof_number(system.number(),0,0), system.solution->el(local_node_ptr->dof_number(system.number(),0,0))/avg_weight);
}
system.solution->close();

This gives the right solution:

Is there a way of reducing/simplifying this code further?

[jwpeterson]

Is there a way of reducing/simplifying this code further?

I guess it depends a bit on the format of your input. In your example, you used the constant "1" but in practice I guess there will be a CONSTANT, MONOMIAL field variable that you want to average in this way? In that case, we do a very similar averaging process in EquationSystems::build_solution_vector() and EquationSystems::build_parallel_solution_vector() (see src/systems/equation_systems.C) so you may want to check whether that code would just work for your case as-is.