Re: [deal.II] Compilation error in v9.5.1

2023-07-21 Thread Wolfgang Bangerth

On 7/21/23 15:26, Tomáš Levý wrote:


Any advice please?


Tomas:
that's awkward -- I guess nobody tests deal.II with PETSc but without MPI :-(

Can you try if it compiles for you if you simply comment out the entire 
content of the export_to_ghosted_array_start() and 
export_to_ghosted_array_finish() functions? Or are there other functions that 
then still don't compile?


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/9ef2a7a5-ed3e-935a-a16f-7bc710e1863c%40colostate.edu.


Re: [deal.II] What do the AffineConstraints::distribute_local_to_global(...) functions do on a matrix/vector level?

2023-07-21 Thread Wolfgang Bangerth

On 7/21/23 04:48, Luca Heltai wrote:

Boundary conditions and hanging nodes have a very delicate interplay.


As just another data point, step-26 re-builds the matrix in every time step 
just so that we can apply boundary conditions and hanging node constraints 
from scratch.


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/03b57226-f66e-011a-3adb-58cdf995d661%40colostate.edu.


[deal.II] How to transpose a distributed PETScWrappers::MPI::SparseMatrix

2023-07-21 Thread Laryssa Abdala
Hello everyone,

It seems like PETScWrappers::MatrixBase::transpose()

only
works in serial. Is there a way to transpose a distributed
PETScWrappers::MPI::SparseMatrix in place? It seems like the way to do this
is to:
Mat temp;
MatTranspose( petsc_wrapper_matrix,  MAT_INITIAL_MATRIX, );

>From there, is there a way to copy a Mat into a
PETScWrappers::MPI::SparseMatrix? Please let me know if I am
missing something or if you have other suggestions.

Thank you,
Laryssa

-- 
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/CAGtr4iBT5iQ-c0cJMa6AQWFVbahTONWUjzmE1yfdT52BvRcj_w%40mail.gmail.com.


[deal.II] Compilation error in v9.5.1

2023-07-21 Thread Tomáš Levý
Dear all,

I'm trying to compile the 9.5.1 version without MPI on Ubuntu 20.04 but I'm 
getting the following error:

/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc: In member 
function ‘void 
dealii::PETScWrappers::CommunicationPattern::export_to_ghosted_array_start(const
 
dealii::ArrayView&, const 
dealii::ArrayView&) const’:
/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc:264:37: error: 
‘mpi_type_id_for_type’ is not a member of ‘dealii::Utilities::MPI’
  264 | auto datatype = Utilities::MPI::mpi_type_id_for_type;
  | ^~~~
/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc:264:64: error: 
expected primary-expression before ‘>’ token
  264 | auto datatype = Utilities::MPI::mpi_type_id_for_type;
  |^
/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc:264:65: error: 
expected primary-expression before ‘;’ token
  264 | auto datatype = Utilities::MPI::mpi_type_id_for_type;
  | ^
/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc: In member 
function ‘void 
dealii::PETScWrappers::CommunicationPattern::export_to_ghosted_array_finish(const
 
dealii::ArrayView&, const 
dealii::ArrayView&) const’:
/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc:280:37: error: 
‘mpi_type_id_for_type’ is not a member of ‘dealii::Utilities::MPI’
  280 | auto datatype = Utilities::MPI::mpi_type_id_for_type;
  | ^~~~
/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc:280:64: error: 
expected primary-expression before ‘>’ token
  280 | auto datatype = Utilities::MPI::mpi_type_id_for_type;
  |^
/home/tom/p/dealii/source/lac/petsc_communication_pattern.cc:280:65: error: 
expected primary-expression before ‘;’ token
  280 | auto datatype = Utilities::MPI::mpi_type_id_for_type;

and so on. 

mpi_type_id_for_type is defined in source/base/mpi.cc, but inside #ifdef 
... #endif:

#ifdef DEAL_II_WITH_MPI
// Provide definitions of template variables for all valid 
instantiations.
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type;
template const MPI_Datatype mpi_type_id_for_type>;
template const MPI_Datatype mpi_type_id_for_type>;
#endif

which should't be defined without MPI if I'm correct so the compilation 
error is about right.

Version 9.4.x compiles without problem.

Any advice please? Thanks

Tom

-- 
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/b5b1160e-96df-443a-aaf9-7384b41e659dn%40googlegroups.com.


Re: [deal.II] Q: Question about extracting part of a vector

2023-07-21 Thread Najwa Alshehri
Dear Daniel, your answer makes sense, Finally it worked.

Dear Luca, thank you for your answer. Indeed, the ExactSolution is a class 
I built myself. I followed your suggestion and it was the magic line that 
solved the issue. *However,* say that I have in another code instead of 
using a ExactSolution that I build myself, I used an FeFieldFunction for a 
reference solution, how would I make this function allocate two components?

Best,
Najwa

On Friday, July 21, 2023 at 1:54:42 PM UTC+3 luca@gmail.com wrote:

> Dear Najwa,
>
> how many components do the `ExactSolution2` function has?
>
> The error complains about `exact_solution.n_components, n_components` not 
> being equal, in particular it tells you that `ExactSolution2` only has 
> one component. 
>
> If you built the ExactSolution class yourself, make sure at construction 
> time you tell Function to allocate two components:
>
> ExactSolution2::ExactSolution2( … ) 
> : Functions::Function(2)
>
> Best,
> Luca.
>
>
> > On 20 Jul 2023, at 12:15, Najwa Alshehri  wrote:
> > 
> > Hello again,
> > 
> > I have a follow-up question. Does this ComponentSelectFunction work also 
> with vectors that are not written as blocked vectors? I have applied it 
> before when I was dealing with blocked vectors and it worked perfectly. 
> > 
> > So, I did this, 
> > const ComponentSelectFunction primal_mask(0,2);
> > Later,
> > VectorTools::integrate_difference(omega2_dh,
> > u_omega2,
> > ExactSolution2(),
> > cellwise_errors2,
> > quadrature,
> > VectorTools::L2_norm
> > ,_mask);
> > const double u2_l2_error =
> > VectorTools::compute_global_error(triangulation_omega2,
> > cellwise_errors2,
> > VectorTools::L2_norm);
> > 
> > And I got the following error!!
> > An error occurred in line <455> of file 
> <../include/deal.II/numerics/vector_tools_integrate_difference.templates.h> 
> in function
> > void dealii::VectorTools::internal::do_integrate_difference(const 
> dealii::hp::MappingCollection&, const 
> dealii::DoFHandler&, const InVector&, const 
> dealii::Function&, OutVector&, 
> const dealii::hp::QCollection&, const dealii::VectorTools::NormType&, 
> const dealii::Function*, double) [with int dim = 2; int spacedim 
> = 2; InVector = dealii::Vector; OutVector = dealii::Vector; 
> typename InVector::value_type = double]
> > The violated condition was: 
> > 
> ::dealii::deal_II_exceptions::internals::compare_for_equality(exact_solution.n_components,
>  
> n_components)
> > Additional information: 
> > Dimension 1 not equal to 2.
> > 
> > Obviously, we have dimensionality mismatching between u_omega2 and the 
> exact solution. which means that the mask is not really picking up the 
> first component of the solution u_omega2. Do you have any suggestions to 
> fix this issue?
> > 
> > Appreciate your help,
> > Najwa
> > On Thursday, July 20, 2023 at 11:48:04 AM UTC+3 Najwa Alshehri wrote:
> > Thank you Daniel for the clear quick answer. I will follow it.
> > 
> > Best,
> > Najwa
> > 
> > On Wednesday, July 19, 2023 at 5:16:14 PM UTC+3 d.arnd...@gmail.com 
> wrote:
> > Najwa,
> > 
> > The documentation of VectorTools::integrate_difference(
> https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a676190d2c897ac5da68a9c460fa95832)
>  
> says
> > 
> > 
> > The additional argument weight allows to evaluate weighted norms. The 
> weight function may be scalar, establishing a spatially variable weight in 
> the domain for all components equally. This may be used, for instance, to 
> only integrate over parts of the domain. The weight function may also be 
> vector-valued, with as many components as the finite element: Then, 
> different components get different weights. A typical application is when 
> the error with respect to only one or a subset of the solution variables is 
> to be computed, in which case the other components would have weight values 
> equal to zero. The ComponentSelectFunction class is particularly useful for 
> this purpose as it provides such a "mask" weight. The weight function is 
> expected to be positive, but negative values are not filtered. The default 
> value of this function, a null pointer, is interpreted as "no weighting 
> function", i.e., weight=1 in the whole domain for all vector components 
> uniformly.
> > 
> > Best,
> > Daniel
> > 
> > On Wed, Jul 19, 2023 at 8:01 AM Najwa Alshehri  
> wrote:
> > Dear group members,
> > 
> > I have one question, 
> > 
> > I am trying to calculate the L2 norm of the error of a solution. In 
> particular, I have a vector of the solution u_omega and a FeSystem with two 
> fes.
> > I am interested in computing the L2 norm of the error related to the 
> first fe using integrate_difference function. (Note here u_omega is not a 
> blocked vector, however, it has two solutions stacked together in one 
> vector). Can I extract the solution of the first part from u_omega?
> > 
> > 
> > Thank you in advance for your help.
> > Best,
> > Najwa
> > 
> > -- 
> > The deal.II project is located 

Re: [deal.II] How to set refinement flag on the neighbouring cell when it is owned by another processor?

2023-07-21 Thread Abbas Ballout
I ended up passing the assemble_own_interior_faces_both flag to mesh-worker 
mesh loop. 
I will probably have to use these Consensus algorithms sooner or later 
anyway. 

Thank you. 
Abbas 

On Thursday, July 13, 2023 at 3:17:59 AM UTC+2 Wolfgang Bangerth wrote:

> On 7/11/23 08:23, Abbas Ballout wrote:
> > I don't want to just set a refinement flag for my cell but I also want 
> to set 
> > the refinement
> > flag for the neighbouring cell too and that is what 
> ncellset_refine_flag(); 
> > doesThis doesn't work if the neighbouring cell is owned by another core
>
> Correct. That's because you can't set information on cells you don't own, 
> in 
> particular because it isn't clear what should happen if different 
> non-owning 
> processes (for which this cell is a ghost cell) disagree in their 
> decisions on 
> what they want to set.
>
> If you want to do this kind of thing, you probably want to collect a list 
> of 
> (cell_id, value) that you want to send to each of the neighboring 
> processes, 
> then send it via Utilities::MPI::isend(). Then you'd receive this kind of 
> list 
> on each process via Utilities::MPI::irecv(), unpack it, and then set the 
> value 
> on each of the cells in question.
>
> If you want a higher-level interface, you can use the functions in 
> Utilities::MPI::ConsensusAlgorithms.
>
> Best
> Wolfgang
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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/b7efbdde-5de5-4b93-a4ef-91152ba63a99n%40googlegroups.com.


Re: [deal.II] What do the AffineConstraints::distribute_local_to_global(...) functions do on a matrix/vector level?

2023-07-21 Thread Luca Heltai
Kyle, 

what you are asking is a very difficult question. If you take a look at the 
(new) step-86 branch (https://github.com/dealii/dealii/pull/15540), you’ll find 
a similar approach to what you describe. 

Everything you said is correct, with possibly one problem (that I have 
encountered very recently, and may be related to what you are observing).

Boundary conditions and hanging nodes have a very delicate interplay. In 
particular, when you *build* the constraints, if one of the hanging nodes 
(dependent) has one of its dependencies on the boundary, then when you add 
non-homogeneous boundary conditions to the constraints, the hanging node 
constraints are also modified to take into account the values of the boundary 
conditions. This computation is done *when you compute the constraints for BC*, 
and it holds the same values for the rest of the life of the affine 
constraints. When you `distribute`, the constraints are computed with the BC 
that you had at the beginning. 

This information is retained even if you change later your (time dependent) 
boundary conditions, and may result in the hanging nodes being compute in a 
wrong manner (in your solution, you may see discontinuous solutions instead of 
continuous ones where hanging nodes interact with boundary nodes).

This affects the final call to “distribute”, not the 
“distribute_local_to_global” (with false as last argument), since that call 
condensates away the constraints. So your process is correct, but `distribute` 
does not do the right job, because the constraints on the hanging nodes have 
not been updated to take into account the new boundary conditions.

To verify that this is the case, try doing the following when you update the 
boundary conditions: clear the constraints, and recreate them from scratch, and 
see if this solves your problem. The `distribute_local_to_global` will behave 
exactly in the same way, but `distribute` should now do the right thing.

L.

> On 20 Jul 2023, at 22:28, Kyle Schwiebert  wrote:
> 
> Hello all,
> 
> Thank you for taking the time to look at my question. First, I'll ask a 
> couple of basic questions about the built-in functions, and then I'll give a 
> few details of why I ask.
> 
> Does the inhomogeneity-handling call distribute_local_to_global(...,false) do 
> the following:
> 
> Say that we have a constraint x_i = a. I would guess that the function takes 
> the following steps: (matrix is M, RHS is right hand side)
> 1) Set M(i,i) = L and set RHS(i) = 0, where L is somehow "nice" for the 
> matrix's spectrum.
> 2) Set the rest of the entries of the ith row and column of M = 0. Perhaps 
> technically we don't bother to zero the column since we plan to call 
> AffineConstraints::set_zero.
> 3) For each DoF j which is coupled to DoF i and is not itself constrained, we 
> do RHS(j) -= a*M(j,i).
> 
> Now considering the same constraint with instead a call to 
> distribute_local_to_global(local_rhs,dof_inds,global_rhs) I would guess that 
> we just set RHS(i) = 0 and do nothing else special.
> 
> Most importantly, is the above correct?
> 
> If so, I am having trouble understanding where I might be going wrong. I'm 
> solving a time dependent coupled problem with a decoupled solver. By 
> decoupled solver I mean that information from the second domain only appears 
> in the right hand side. To make matters one step more complex, one term in my 
> matrix is time-dependent and so I do not recompute the whole matrix. My 
> procedure is this:
> 
> 1) Assemble the static part of the matrix with a zero right hand side. This 
> matrix will never be touched again. The vector, RHS_bc, which will also never 
> be touched again contains no nonzero entries except for those created by the 
> call to distribute_local_to_global(...,false).
> 
> 2) I then assemble the time-dependent terms into a separate matrix and 
> vector, RHS_td, using distribute_local_to_global(...,false) and add the 
> RHS_bc to RHS_td.
> 
> 3) I then assemble the coupled terms from the other solver into RHS_c with 
> distribute_local_to_global(local_rhs,dof_inds,RHS_c).
> 
> 4) I then add RHS_c to RHS_td and solve (being careful to call 
> constraints.set_zero(solution) before the solve and 
> constraints.distribute(solution) after the solve). I am of course using a 
> custom class similar to the linear operator classes to bundle the vmult() 
> outputs of the static and time-dependent matrices.
> 
> However, I am seeing blowup in my solution and I'm seeing some indication 
> that the inhomogeneous boundary conditions are to blame. Note that I also 
> tried calling constraints.set_zero(RHS_c) in step 4 above and as I expected 
> that had no change (since presumably constrained entries are already zero).
> 
> Thank you again for taking the time to engage with my question and in advance 
> for any tips you are able to offer.
> 
> -Kyle
> 
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
>