Re: [deal.II] Getting RHS values at nodes with DBC

[Wolfgang Bangerth]

> In my case I solve a simple laplace problem, that represent groundwater flow 
> with non zero dirichlet boundary conditions
> My goal is to estimate how much water comes in and out from the dirichlet BC 
> nodes.

You're asking the wrong question. What you are looking for is a *flux*, which 
is an integral over part of the boundary. That's how you should compute it -- 
which implies you'll want to use FEFaceValues as always when face integrals 
are involved. Specifically, the flux through a single point is of course zero.


> The way I could do this, is after I solve the reduced version of the full 
> linear system A*U = F
> to calculate
> DBC = A*U-F

I'm not clear what you're trying to do here, but if A*U=F, then DBC is of 
course zero. I suspect that you mean for either A, U, or F to be something 
slightly different than what you had in the linear system above, but it's not 
clear to me -- though it may also not matter in view of what I mentioned above 
:-)

Best
  W.

-- 

Wolfgang Bangerth  email: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1abaa9aa-1464-28e0-22fc-56d5e379fa4e%40colostate.edu.

Re: [deal.II] VectorTools::Compress() fails after mesh_loop and distribute_local_to_global

['Maxi Miller' via deal.II User Group]

Will test that, just need to recompile my debug-version of deal.II

Am Samstag, 31. August 2019 18:17:34 UTC+2 schrieb Timo Heister:
>
> This is an interesting problem that might be a deal.II bug. I can't 
> see how it has anything to do with mesh_loop, but it looks like 
> operator=(0) is not setting last_action. 
>
> That said, I am confused because .reinit() should set last_action to 
> "Zero". 
>
> Can you go into the line of the assert in trilinos_vector.cc and print 
> the value of "last_action" on all processors? 
>
> If you can not figure it out, can you make a minimal complete example I 
> can try? 
>
>
> On Tue, Aug 27, 2019 at 6:59 AM 'Maxi Miller' via deal.II User Group 
> > wrote: 
> > 
> > I am trying to assemble my residual/RHS by using mesh_loop: 
> > evaluation_point = present_solution; 
> > extension_vector = newton_update; 
> > extension_vector *= alpha; 
> > 
> > evaluation_point += extension_vector; 
> > 
> > residual.reinit(locally_owned_dofs, mpi_communicator); 
> > residual_result.reinit(locally_owned_dofs, mpi_communicator); 
> > residual = 0; 
> > residual_result = 0; 
> > 
> > const QGauss  quadrature_formula(fe.degree + gauss_size); 
> > 
> > 
> > using CellFilter = FilteredIterator DoFHandler<2>::active_cell_iterator>; 
> > std::cout << "Looping over all cells\n"; 
> > 
> > MeshWorker::mesh_loop(dof_handler.begin_active(), 
> >   dof_handler.end(), 
> >   std::bind( ad_type_code>::local_assemble_residual, 
> > this, 
> > std::placeholders::_1, 
> > std::placeholders::_2, 
> > std::placeholders::_3), 
> >   std::bind( ad_type_code>::copy_local_to_global_residual, 
> > this, 
> > std::placeholders::_1), 
> >   ResidualScratchData(fe, 
> quadrature_formula, update_values | update_gradients | 
> update_quadrature_points | update_JxW_values, I_val, I_val_old), 
> >   ResidualCopyData()); 
> > 
> > std::cout << "Compressing residual\n"; 
> > residual.compress(VectorOperation::add); 
> > std::cout << "Residual compressed\n"; 
> > boundary_constraints.set_zero(residual); 
> > std::cout << "Boundary constraints set to zero\n"; 
> > residual_result = residual; 
> > 
> > 
> > 
> > with 
> > copy_local_to_global_residual 
> > 
> > defined as 
> > void MinimalSurfaceProblem ad_type_code>::copy_local_to_global_residual(const ResidualCopyData ) 
> > { 
> > 
> boundary_constraints.distribute_local_to_global(data.local_residual, 
> > 
> data.local_dof_indices, 
> > residual); 
> > } 
> > 
> > Unfortunately, my program fails after some iterations with the error 
> >  
> > An error occurred in line <624> of file 
> <~/Downloads/git-files/dealii/source/lac/trilinos_vector.cc> in function 
> > void 
> dealii::TrilinosWrappers::MPI::Vector::compress(dealii::VectorOperation::values)
>  
>
> > The violated condition was: 
> > result.max - result.min < 1e-5 
> > Additional information: 
> > Not all processors agree whether the last operation on this vector 
> was an addition or a set operation. This will prevent the compress() 
> operation from succeeding. 
> > 
> > Stacktrace: 
> > --- 
> > #0  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::TrilinosWrappers::MPI::Vector::compress(dealii::VectorOperation::values)
>  
>
> > #1  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
> Silicon_Verburg::physics_equations, 
> (dealii::Differentiation::AD::NumberTypes)3>::compute_residual(double, 
> dealii::TrilinosWrappers::MPI::Vector&) 
> > #2  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
> Silicon_Verburg::physics_equations, 
> (dealii::Differentiation::AD::NumberTypes)3>::recalculate_step_length() 
> > #3  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
> Silicon_Verburg::physics_equations, 
> (dealii::Differentiation::AD::NumberTypes)3>::update_newton_solution(double 
> const&, double const&) 
> > #4  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
> Silicon_Verburg::physics_equations, 
> (dealii::Differentiation::AD::NumberTypes)3>::run() 
> > #5  bin/TTM-equation: main 
> >  
> > 
> > Based on the output I assume that the 
> distribute_local_to_global-function is the problem (when running with 8 
> threads, I get 8 lines with "Looping over all cells", but only 7 with 
>