Hi all,
I have been tracking an error in my code in applying Dirichlet
boundary conditions. I use the ConstraintMatrix method of applying the
boundary conditions as in step-22. The fluid domain equation that I am
solving is the Laplace equation, and the boundary conditions are time
dependent (a mixture of Robin and Dirichlet conditions).
My approach to the time stepping is to build the bulk of the system
matrix (including fluid equation and the unchanging parts of the Robin
boundary conditions) into initial_matrix in a method
setup_initial_system. I copy this unchanging part into the actual
system_matrix at each time step in the assemble_system method and then
proceed with building the remainder of the system_matrix.
I have tracked the problem with my Dirichlet boundary conditions to
the distribute_local_to_global method, which does not seem to be
usable if the system_matrix already contains data before assembly
begins.
I verified that two different ways of distributing the cell_matrix
data and constraints give different results by using the following
snippets (one alternative commented out). The commented out option
produces the expected results and the other does not.
cell->get_dof_indices(local_dof_indices);
/* for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}*/
all_constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix, system_rhs);
}
/* all_constraints.condense (system_matrix, system_rhs);*/
It seems apparent that distribute_local_to_global is not valid if the
matrix and rhs already contain data. I believe that the actual method
run is found in constant_matrix.templates.h and I made some notes as
below:
Line 1173: start of relevant
Line 1205: Finds average of local_matrix diagonal, but does not take
into account entries already in matrix.
Line 1311: Uses averaged diagonal
Line 1370: matrix_ptr[constraint_lines[i].first] possibly zero because
majority of the fluid already built in global_matrix.
Line 1412: Similar line to above, but for the indirect references case.
Line 1479: Another similar line for indirect references case.
The distribute_local_to_global method seems complex already, without
taking into account the possibility that system_matrix already
contains data. Is it possible to extend it to cover this case? I would
prefer to use this method to the commented out alternative since it
extends to MPI Vectors and Matrices easily.
Are there any suggestions of alternative approaches? I cannot use the
FilteredMatrix class since I add entries during my assemble_system
method. I currently use system_matrix.copy_from(initial_matrix); to
copy the matrix, and have thought of copying from initial_matrix into
each cell_matrix and then using distribute_local_to_global but this
seems inefficient.
Any thoughts or suggestions are appreciated!
Michael
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
