Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-04 Thread 'develo...@googlemail.com' via deal.II User Group 
Using the merged commit my runtime decreased now from 120 sec. per call to 
~1-2 sec. per call, thanks for the fast help!

develo...@googlemail.com schrieb am Donnerstag, 3. September 2020 um 
14:13:37 UTC+2:

> I'm following the discussion, and will retest my program as soon as the 
> commit has been approved and merged. Thanks!
>
> Martin Kronbichler schrieb am Donnerstag, 3. September 2020 um 13:55:40 
> UTC+2:
>
>> Indeed the other function was fast. I looked into the function already 
>> and have a fix posted:
>>
>> https://github.com/dealii/dealii/pull/10883
>>
>> Best,
>> Martin
>>
>>

-- 
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/6cf635b0-0384-4fa4-97a7-c6935721f5bfn%40googlegroups.com.

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-03 Thread 'develo...@googlemail.com' via deal.II User Group 
I'm following the discussion, and will retest my program as soon as the 
commit has been approved and merged. Thanks!

Martin Kronbichler schrieb am Donnerstag, 3. September 2020 um 13:55:40 
UTC+2:

> Indeed the other function was fast. I looked into the function already 
> and have a fix posted:
>
> https://github.com/dealii/dealii/pull/10883
>
> Best,
> Martin
>
>

-- 
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/2d5142e6-ede2-4b40-b3b9-34123ff3e635n%40googlegroups.com.

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-03 Thread Martin Kronbichler 
Indeed the other function was fast. I looked into the function already 
and have a fix posted:


https://github.com/dealii/dealii/pull/10883

Best,
Martin

--
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/377979f8-d10e-065e-f157-abecadebbe92%40gmail.com.

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-03 Thread Wolfgang Bangerth 

On 9/3/20 1:16 AM, 'Maxi Miller' via deal.II User Group wrote:


Again, extract_dofs takes most of the time here. Is there an alternative 
approach for that, too?


There is a variation of extract_dofs() that returns its results in the form of 
a std::vector (as one of the arguments) instead of an IndexSet. Could 
you try and see whether that is substantially faster?


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/e4b9665c-7423-ff8c-10ec-951c18085c32%40colostate.edu.

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-03 Thread 'Maxi Miller' via deal.II User Group 
Another issue I have is for the following code piece, where I basically
have to rescale every part of the solution by a different factor (due to
using unitless variables in the calculation):
const ComponentMask surface_A_mask = fe.component_mask(surface_A);
const ComponentMask surface_B_mask = fe.component_mask(surface_B);
const ComponentMask surface_C_mask = fe.component_mask(surface_C);
const ComponentMask surface_D_mask = fe.component_mask(surface_D);

IndexSet A_value_mask = DoFTools::extract_dofs(dof_handler, surface_A_mask);
IndexSet B_value_mask = DoFTools::extract_dofs(dof_handler, surface_B_mask);
IndexSet C_value_mask = DoFTools::extract_dofs(dof_handler, surface_C_mask);
IndexSet D_value_mask = DoFTools::extract_dofs(dof_handler, surface_D_mask);

IndexSet locally_owned_elements = scaled_solution.locally_owned_elements();

A_value_mask.compress();
B_value_mask.compress();
C_value_mask.compress();
D_value_mask.compress();
locally_owned_elements.compress();

vector_t unified_test_vector(dof_handler.locally_owned_dofs(),
mpi_communicator);

for(auto index : scaled_solution.locally_owned_elements()) {

if(A_value_mask.is_element(locally_owned_elements.index_within_set(index)))
  {
  unified_test_vector(index) =
unitless_recalculator.scale_A_solution(scaled_solution(index));
 }

if(B_value_mask.is_element(locally_owned_elements.index_within_set(index)))
 {
  unified_test_vector(index) =
unitless_recalculator.scale_B_solution(scaled_solution(index));
 }

if(C_value_mask.is_element(locally_owned_elements.index_within_set(index)))
 {
  unified_test_vector(index) =
unitless_recalculator.scale_C_solution(scaled_solution(index));
 }

if(D_value_mask.is_element(locally_owned_elements.index_within_set(index)))
 {
  unified_test_vector(index) =
unitless_recalculator.scale_D_solution(scaled_solution(index));
 }
}

Again, extract_dofs takes most of the time here. Is there an alternative
approach for that, too?

Am Do., 3. Sept. 2020 um 08:56 Uhr schrieb 'develo...@googlemail.com' via
deal.II User Group :

> Yes, that worked, thanks!
> Now, when comparing, it takes 132 seconds for my initial version, and 0.41
> seconds for the version you proposed. Why is my approach so much more
> expensive?
> Thanks!
>
> Jean-Paul Pelteret schrieb am Mittwoch, 2. September 2020 um 20:34:56
> UTC+2:
>
>> I was close — it’s
>> fe.get_unit_support_points() //
>> https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#a5b35a290aa7dd7562911a92a13b11fee
>>
>> Moreover, defining fe_values as const disables the reinit()-call, and I
>> have to call cell->is_locally_owned(), instead of cell->locally_owned().
>>
>>
>> Coding in an email client is unforgiving, so I’m glad that you were able
>> to work past these errors.
>>
>> Best,
>> Jean-Paul
>>
>> On 02 Sep 2020, at 14:11, 'develo...@googlemail.com' via deal.II User
>> Group  wrote:
>>
>> I'm using FESystem for defining fe, and therefore I get
>> "unit_support_points is protected within this context". Is there another
>> way to get access to those points?
>> Moreover, defining fe_values as const disables the reinit()-call, and I
>> have to call cell->is_locally_owned(), instead of cell->locally_owned().
>> Thanks!
>>
>> Jean-Paul Pelteret schrieb am Dienstag, 1. September 2020 um 20:28:18
>> UTC+2:
>>
>>> Hi,
>>>
>>> 30-40% of the runtime for this operation seems a bit excessive to me.
>>> Maybe you should try using the FEValuesExtractors in conjunction with
>>> FEValues in order to get the solution at the quadrature points. If you want
>>> the solution at the support points, then you can just create a quadrature
>>> rule that has the quadrature points at the support points and initialise
>>> your FEValues object with it.
>>>
>>> So this bit of pseudo-code illustrates roughly how it could be done:
>>>
>>> //
>>> https://dealii.org/current/doxygen/deal.II/classQuadrature.html#ac002803e6ec722da5a0de102574110c9
>>> //
>>> https://dealii.org/current/doxygen/deal.II/classFiniteElement.html#ab4f6e0c83686b918fbb92716ead92313
>>> const Quadrature q_cell_support_points (fe.unit_support_points());
>>> // Assumes that your FE has unit support points; could also use the unit
>>> cell vertex positions, something else...
>>> const FEValues fe_values(fe, q_cell_support_points, update_values);
>>> const FEValuesExtractors::Scalar surface_A(0);
>>>
>>> Vector solution_vector = …; // If this is an MPI vector, then
>>> still use the same basic approach
>>>
>>> double max_value = std::numeric_limits::min();
>>> for (auto cell : dof_handler.active_cell_iterators())
>>> {
>>>   if (!cell->locally_owned()) continue;
>>>
>>>   fe_values.reinit(cell);
>>>
>>>   //
>>> https://dealii.org/current/doxygen/deal.II/classFEValuesViews_1_1Scalar.html#a7aa6b0275facea985f23a64ad690d9b1
>>>   std::vector q_point_solution_values
>>> (fe_values.n_quadrature_points);
>>>   fe_values[surface_A].get_function_values(solution_vector ,
>>> q_point_

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-02 Thread 'develo...@googlemail.com' via deal.II User Group 
Yes, that worked, thanks!
Now, when comparing, it takes 132 seconds for my initial version, and 0.41 
seconds for the version you proposed. Why is my approach so much more 
expensive?
Thanks!

Jean-Paul Pelteret schrieb am Mittwoch, 2. September 2020 um 20:34:56 UTC+2:

> I was close — it’s 
> fe.get_unit_support_points() // 
> https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#a5b35a290aa7dd7562911a92a13b11fee
>
> Moreover, defining fe_values as const disables the reinit()-call, and I 
> have to call cell->is_locally_owned(), instead of cell->locally_owned().
>
>
> Coding in an email client is unforgiving, so I’m glad that you were able 
> to work past these errors.
>
> Best,
> Jean-Paul
>
> On 02 Sep 2020, at 14:11, 'develo...@googlemail.com' via deal.II User 
> Group  wrote:
>
> I'm using FESystem for defining fe, and therefore I get 
> "unit_support_points is protected within this context". Is there another 
> way to get access to those points?
> Moreover, defining fe_values as const disables the reinit()-call, and I 
> have to call cell->is_locally_owned(), instead of cell->locally_owned().
> Thanks!
>
> Jean-Paul Pelteret schrieb am Dienstag, 1. September 2020 um 20:28:18 
> UTC+2:
>
>> Hi,
>>
>> 30-40% of the runtime for this operation seems a bit excessive to me. 
>> Maybe you should try using the FEValuesExtractors in conjunction with 
>> FEValues in order to get the solution at the quadrature points. If you want 
>> the solution at the support points, then you can just create a quadrature 
>> rule that has the quadrature points at the support points and initialise 
>> your FEValues object with it.
>>
>> So this bit of pseudo-code illustrates roughly how it could be done:
>>
>> // 
>> https://dealii.org/current/doxygen/deal.II/classQuadrature.html#ac002803e6ec722da5a0de102574110c9
>> // 
>> https://dealii.org/current/doxygen/deal.II/classFiniteElement.html#ab4f6e0c83686b918fbb92716ead92313
>> const Quadrature q_cell_support_points (fe.unit_support_points()); 
>> // Assumes that your FE has unit support points; could also use the unit 
>> cell vertex positions, something else...
>> const FEValues fe_values(fe, q_cell_support_points, update_values);
>> const FEValuesExtractors::Scalar surface_A(0);
>>
>> Vector solution_vector = …; // If this is an MPI vector, then 
>> still use the same basic approach
>>
>> double max_value = std::numeric_limits::min();
>> for (auto cell : dof_handler.active_cell_iterators())
>> {
>>   if (!cell->locally_owned()) continue;
>>
>>   fe_values.reinit(cell);
>>
>>   // 
>> https://dealii.org/current/doxygen/deal.II/classFEValuesViews_1_1Scalar.html#a7aa6b0275facea985f23a64ad690d9b1
>>   std::vector q_point_solution_values 
>> (fe_values.n_quadrature_points);
>>   fe_values[surface_A].get_function_values(solution_vector , 
>> q_point_solution_values);
>>
>>   // This will get you the maximum of the solution at the queried points
>>   // for the local process.
>>   for (auto entry : q_point_solution_values)
>> max_value = std::max(max_value, entry);
>>
>>   // Since you’re doing the work with a cell loop, if appropriate then 
>> you could also
>>   // retrieve the max for the other fields at the same time.
>> }
>>
>> // If the solution is distributed, then you still need the max over all 
>> processes:
>> // 
>> https://dealii.org/current/doxygen/deal.II/namespaceUtilities_1_1MPI.html#ad2f716b789abe53715d6659f38aa7815
>> max_value = Utilities::MPI::max(max_value, mpi_communicator);
>>
>> I think that’s more or less how you could do it using the established 
>> classes and methodology. Hopefully it would yield you better performance,
>> since you don’t have to fetch a huge index set mask. The other benefit is 
>> that what I’ve written here makes no assumptions about how the coefficients 
>> of the solution vector relate to the physical solution.
>>
>> I hope that this helps you out. It would be great to hear whether or not 
>> this sort of approach yields a performance benefit.
>>
>> Best,
>> Jean-Paul
>>
>>
>> On 31 Aug 2020, at 10:58, 'develo...@googlemail.com' via deal.II User 
>> Group  wrote:
>>
>>
>> I'd like to simulate four coupled PDEs (with variables A, B, C and D) 
>> using deal.II. At each step in time I'm interested in the maximum and 
>> minimum value of each variable. Therefore, I have to loop over my solution, 
>> and gather the maximum and minimum value of each variable by using
>> vector_t local_solution;
>> present_solution.update_ghost_values();
>> local_solution = present_solution;
>>
>> IndexSet locally_owned_elements = local_solution.locally_owned_elements();
>> locally_owned_elements.compress();
>>
>> for(auto index: present_solution.locally_owned_elements()) {
>>  
>> if(A_value_mask.is_element(locally_owned_elements.index_within_set(index))) 
>> {
>>  max_A_value = std::max(max_A_value, 
>> static_cast(local_solution(index)));
>>  min_A_value = std::min(min_A_value, 
>> static_cast(local_solution(ind

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-02 Thread Jean-Paul Pelteret 
I was close — it’s 
fe.get_unit_support_points() // 
https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#a5b35a290aa7dd7562911a92a13b11fee
 


> Moreover, defining fe_values as const disables the reinit()-call, and I have 
> to call cell->is_locally_owned(), instead of cell->locally_owned().

Coding in an email client is unforgiving, so I’m glad that you were able to 
work past these errors.

Best,
Jean-Paul

> On 02 Sep 2020, at 14:11, 'develo...@googlemail.com' via deal.II User Group 
>  wrote:
> 
> I'm using FESystem for defining fe, and therefore I get "unit_support_points 
> is protected within this context". Is there another way to get access to 
> those points?
> Moreover, defining fe_values as const disables the reinit()-call, and I have 
> to call cell->is_locally_owned(), instead of cell->locally_owned().
> Thanks!
> 
> Jean-Paul Pelteret schrieb am Dienstag, 1. September 2020 um 20:28:18 UTC+2:
> Hi,
> 
> 30-40% of the runtime for this operation seems a bit excessive to me. Maybe 
> you should try using the FEValuesExtractors in conjunction with FEValues in 
> order to get the solution at the quadrature points. If you want the solution 
> at the support points, then you can just create a quadrature rule that has 
> the quadrature points at the support points and initialise your FEValues 
> object with it.
> 
> So this bit of pseudo-code illustrates roughly how it could be done:
> 
> // 
> https://dealii.org/current/doxygen/deal.II/classQuadrature.html#ac002803e6ec722da5a0de102574110c9
>  
> 
> // 
> https://dealii.org/current/doxygen/deal.II/classFiniteElement.html#ab4f6e0c83686b918fbb92716ead92313
>  
> 
> const Quadrature q_cell_support_points (fe.unit_support_points()); // 
> Assumes that your FE has unit support points; could also use the unit cell 
> vertex positions, something else...
> const FEValues fe_values(fe, q_cell_support_points, update_values);
> const FEValuesExtractors::Scalar surface_A(0);
> 
> Vector solution_vector = …; // If this is an MPI vector, then still 
> use the same basic approach
> 
> double max_value = std::numeric_limits::min();
> for (auto cell : dof_handler.active_cell_iterators())
> {
>   if (!cell->locally_owned()) continue;
> 
>   fe_values.reinit(cell);
> 
>   // 
> https://dealii.org/current/doxygen/deal.II/classFEValuesViews_1_1Scalar.html#a7aa6b0275facea985f23a64ad690d9b1
>  
> 
>   std::vector q_point_solution_values (fe_values.n_quadrature_points);
>   fe_values[surface_A].get_function_values(solution_vector , 
> q_point_solution_values);
> 
>   // This will get you the maximum of the solution at the queried points
>   // for the local process.
>   for (auto entry : q_point_solution_values)
> max_value = std::max(max_value, entry);
> 
>   // Since you’re doing the work with a cell loop, if appropriate then you 
> could also
>   // retrieve the max for the other fields at the same time.
> }
> 
> // If the solution is distributed, then you still need the max over all 
> processes:
> // 
> https://dealii.org/current/doxygen/deal.II/namespaceUtilities_1_1MPI.html#ad2f716b789abe53715d6659f38aa7815
>  
> 
> max_value = Utilities::MPI::max(max_value, mpi_communicator);
> 
> I think that’s more or less how you could do it using the established classes 
> and methodology. Hopefully it would yield you better performance,
> since you don’t have to fetch a huge index set mask. The other benefit is 
> that what I’ve written here makes no assumptions about how the coefficients 
> of the solution vector relate to the physical solution.
> 
> I hope that this helps you out. It would be great to hear whether or not this 
> sort of approach yields a performance benefit.
> 
> Best,
> Jean-Paul
> 
> 
> 
>> On 31 Aug 2020, at 10:58, 'develo...@googlemail.com 
>> ' via deal.II User Group > > wrote:
>> 
> 
>> 
>> I'd like to simulate four coupled PDEs (with variables A, B, C and D) using 
>> deal.II. At each step in time I'm interested in the maximum and minimum 
>> value of each variable. Therefore, I have to loop over my solution, and 
>> gather the maximum and minimum value of each variable by using
>> vector_t local_solution;
>> present_solution.update_ghost_values();
>> local_solution = present_solution;
>> 
>> IndexSet locally_owned_elements = local_solution.locally_owned_elements();
>> locally_owned_elements.compress();
>> 
>> for(auto index: present_solution.locally_owned_elements()) {
>>  
>> if(

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-02 Thread 'develo...@googlemail.com' via deal.II User Group 
I'm using FESystem for defining fe, and therefore I get 
"unit_support_points is protected within this context". Is there another 
way to get access to those points?
Moreover, defining fe_values as const disables the reinit()-call, and I 
have to call cell->is_locally_owned(), instead of cell->locally_owned().
Thanks!

Jean-Paul Pelteret schrieb am Dienstag, 1. September 2020 um 20:28:18 UTC+2:

> Hi,
>
> 30-40% of the runtime for this operation seems a bit excessive to me. 
> Maybe you should try using the FEValuesExtractors in conjunction with 
> FEValues in order to get the solution at the quadrature points. If you want 
> the solution at the support points, then you can just create a quadrature 
> rule that has the quadrature points at the support points and initialise 
> your FEValues object with it.
>
> So this bit of pseudo-code illustrates roughly how it could be done:
>
> // 
> https://dealii.org/current/doxygen/deal.II/classQuadrature.html#ac002803e6ec722da5a0de102574110c9
> // 
> https://dealii.org/current/doxygen/deal.II/classFiniteElement.html#ab4f6e0c83686b918fbb92716ead92313
> const Quadrature q_cell_support_points (fe.unit_support_points()); // 
> Assumes that your FE has unit support points; could also use the unit cell 
> vertex positions, something else...
> const FEValues fe_values(fe, q_cell_support_points, update_values);
> const FEValuesExtractors::Scalar surface_A(0);
>
> Vector solution_vector = …; // If this is an MPI vector, then 
> still use the same basic approach
>
> double max_value = std::numeric_limits::min();
> for (auto cell : dof_handler.active_cell_iterators())
> {
>   if (!cell->locally_owned()) continue;
>
>   fe_values.reinit(cell);
>
>   // 
> https://dealii.org/current/doxygen/deal.II/classFEValuesViews_1_1Scalar.html#a7aa6b0275facea985f23a64ad690d9b1
>   std::vector q_point_solution_values 
> (fe_values.n_quadrature_points);
>   fe_values[surface_A].get_function_values(solution_vector , 
> q_point_solution_values);
>
>   // This will get you the maximum of the solution at the queried points
>   // for the local process.
>   for (auto entry : q_point_solution_values)
> max_value = std::max(max_value, entry);
>
>   // Since you’re doing the work with a cell loop, if appropriate then you 
> could also
>   // retrieve the max for the other fields at the same time.
> }
>
> // If the solution is distributed, then you still need the max over all 
> processes:
> // 
> https://dealii.org/current/doxygen/deal.II/namespaceUtilities_1_1MPI.html#ad2f716b789abe53715d6659f38aa7815
> max_value = Utilities::MPI::max(max_value, mpi_communicator);
>
> I think that’s more or less how you could do it using the established 
> classes and methodology. Hopefully it would yield you better performance,
> since you don’t have to fetch a huge index set mask. The other benefit is 
> that what I’ve written here makes no assumptions about how the coefficients 
> of the solution vector relate to the physical solution.
>
> I hope that this helps you out. It would be great to hear whether or not 
> this sort of approach yields a performance benefit.
>
> Best,
> Jean-Paul
>
>
> On 31 Aug 2020, at 10:58, 'develo...@googlemail.com' via deal.II User 
> Group  wrote:
>
>
> I'd like to simulate four coupled PDEs (with variables A, B, C and D) 
> using deal.II. At each step in time I'm interested in the maximum and 
> minimum value of each variable. Therefore, I have to loop over my solution, 
> and gather the maximum and minimum value of each variable by using
> vector_t local_solution;
> present_solution.update_ghost_values();
> local_solution = present_solution;
>
> IndexSet locally_owned_elements = local_solution.locally_owned_elements();
> locally_owned_elements.compress();
>
> for(auto index: present_solution.locally_owned_elements()) {
>  
> if(A_value_mask.is_element(locally_owned_elements.index_within_set(index))) 
> {
>  max_A_value = std::max(max_A_value, 
> static_cast(local_solution(index)));
>  min_A_value = std::min(min_A_value, 
> static_cast(local_solution(index)));
> }
> }
>
> To create A_value_mask, I use
> const FEValuesExtractors::Scalar surface_A(0);
> const ComponentMask surface_A_mask = fe.component_mask(surface_A);
> IndexSet A_value_mask = DoFTools::extract_dofs(dof_handler, 
> surface_A_mask);
>
> For 2.4 million elements and four variables (i.e. four calls), this second 
> part (more specifically the call extract_dofs()) takes ~120 sec., which 
> takes up ~30% of the total run time. I have to use such a construct three 
> times in my code, resulting in 360 sec. spend in extracting dofs, and ~40 
> sec. in the other parts of the program. Is that by design, or is there 
> something wrong in my code?
> Thanks!
>
> -- 
> 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

Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

2020-09-01 Thread Jean-Paul Pelteret 
Hi,

30-40% of the runtime for this operation seems a bit excessive to me. Maybe you 
should try using the FEValuesExtractors in conjunction with FEValues in order 
to get the solution at the quadrature points. If you want the solution at the 
support points, then you can just create a quadrature rule that has the 
quadrature points at the support points and initialise your FEValues object 
with it.

So this bit of pseudo-code illustrates roughly how it could be done:

// 
https://dealii.org/current/doxygen/deal.II/classQuadrature.html#ac002803e6ec722da5a0de102574110c9
 

// 
https://dealii.org/current/doxygen/deal.II/classFiniteElement.html#ab4f6e0c83686b918fbb92716ead92313
 

const Quadrature q_cell_support_points (fe.unit_support_points()); // 
Assumes that your FE has unit support points; could also use the unit cell 
vertex positions, something else...
const FEValues fe_values(fe, q_cell_support_points, update_values);
const FEValuesExtractors::Scalar surface_A(0);

Vector solution_vector = …; // If this is an MPI vector, then still use 
the same basic approach

double max_value = std::numeric_limits::min();
for (auto cell : dof_handler.active_cell_iterators())
{
  if (!cell->locally_owned()) continue;

  fe_values.reinit(cell);

  // 
https://dealii.org/current/doxygen/deal.II/classFEValuesViews_1_1Scalar.html#a7aa6b0275facea985f23a64ad690d9b1
 

  std::vector q_point_solution_values (fe_values.n_quadrature_points);
  fe_values[surface_A].get_function_values(solution_vector , 
q_point_solution_values);

  // This will get you the maximum of the solution at the queried points
  // for the local process.
  for (auto entry : q_point_solution_values)
max_value = std::max(max_value, entry);

  // Since you’re doing the work with a cell loop, if appropriate then you 
could also
  // retrieve the max for the other fields at the same time.
}

// If the solution is distributed, then you still need the max over all 
processes:
// 
https://dealii.org/current/doxygen/deal.II/namespaceUtilities_1_1MPI.html#ad2f716b789abe53715d6659f38aa7815
 

max_value = Utilities::MPI::max(max_value, mpi_communicator);

I think that’s more or less how you could do it using the established classes 
and methodology. Hopefully it would yield you better performance,
since you don’t have to fetch a huge index set mask. The other benefit is that 
what I’ve written here makes no assumptions about how the coefficients of the 
solution vector relate to the physical solution.

I hope that this helps you out. It would be great to hear whether or not this 
sort of approach yields a performance benefit.

Best,
Jean-Paul


> On 31 Aug 2020, at 10:58, 'develo...@googlemail.com' via deal.II User Group 
>  wrote:
> 
> 
> I'd like to simulate four coupled PDEs (with variables A, B, C and D) using 
> deal.II. At each step in time I'm interested in the maximum and minimum value 
> of each variable. Therefore, I have to loop over my solution, and gather the 
> maximum and minimum value of each variable by using
> vector_t local_solution;
> present_solution.update_ghost_values();
> local_solution = present_solution;
> 
> IndexSet locally_owned_elements = local_solution.locally_owned_elements();
> locally_owned_elements.compress();
> 
> for(auto index: present_solution.locally_owned_elements()) {
>  
> if(A_value_mask.is_element(locally_owned_elements.index_within_set(index))) {
>  max_A_value = std::max(max_A_value, 
> static_cast(local_solution(index)));
>  min_A_value = std::min(min_A_value, 
> static_cast(local_solution(index)));
> }
> }
> 
> To create A_value_mask, I use
> const FEValuesExtractors::Scalar surface_A(0);
> const ComponentMask surface_A_mask = fe.component_mask(surface_A);
> IndexSet A_value_mask = DoFTools::extract_dofs(dof_handler, surface_A_mask);
> 
> For 2.4 million elements and four variables (i.e. four calls), this second 
> part (more specifically the call extract_dofs()) takes ~120 sec., which takes 
> up ~30% of the total run time. I have to use such a construct three times in 
> my code, resulting in 360 sec. spend in extracting dofs, and ~40 sec. in the 
> other parts of the program. Is that by design, or is there something wrong in 
> my code?
> Thanks!
> 
> -- 
> 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

[deal.II] extract_dofs() extremely slow in for single-thread run?

2020-08-31 Thread 'develo...@googlemail.com' via deal.II User Group 

I'd like to simulate four coupled PDEs (with variables A, B, C and D) using 
deal.II. At each step in time I'm interested in the maximum and minimum 
value of each variable. Therefore, I have to loop over my solution, and 
gather the maximum and minimum value of each variable by using
vector_t local_solution;
present_solution.update_ghost_values();
local_solution = present_solution;

IndexSet locally_owned_elements = local_solution.locally_owned_elements();
locally_owned_elements.compress();

for(auto index: present_solution.locally_owned_elements()) {
 
if(A_value_mask.is_element(locally_owned_elements.index_within_set(index))) 
{
 max_A_value = std::max(max_A_value, 
static_cast(local_solution(index)));
 min_A_value = std::min(min_A_value, 
static_cast(local_solution(index)));
}
}

To create A_value_mask, I use
const FEValuesExtractors::Scalar surface_A(0);
const ComponentMask surface_A_mask = fe.component_mask(surface_A);
IndexSet A_value_mask = DoFTools::extract_dofs(dof_handler, surface_A_mask);

For 2.4 million elements and four variables (i.e. four calls), this second 
part (more specifically the call extract_dofs()) takes ~120 sec., which 
takes up ~30% of the total run time. I have to use such a construct three 
times in my code, resulting in 360 sec. spend in extracting dofs, and ~40 
sec. in the other parts of the program. Is that by design, or is there 
something wrong in my code?
Thanks!

-- 
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/80a29cd6-d2e8-447d-a0a9-89b2e0fc4b32n%40googlegroups.com.