Hi Wasim,
It looks like you're passing your solution vector by copy for each cell
that you're assembling on. You probably want to pass it by reference.
Best,
Jean-Paul
Sent from my mobile device. Please excuse my brevity and any typos.
On Sun, 08 Jan 2023, 14:38 Wasim Niyaz Munshi ce21d400, <
ce21d...@smail.iitm.ac.in> wrote:
> Thank you Prof. Munch.
>
> I am now passing FEValues object as a parameter to the H_plus function. My
> assemble time has decreased from (4-5 seconds) to around 2 seconds, but it
> is still nowhere close to the time required if there is no call to h_plus
> function(say eg. step-3 tutorial problem).
> I don't expect to match that time but is there still anything that I can
> do to reduce the assembly time further?
>
> Thanks and regards
> Wasim
>
> On Sunday, January 8, 2023 at 4:11:27 PM UTC+5:30 Wasim Niyaz Munshi
> ce21d400 wrote:
>
>> Do I need to make this change in both declaration and definition?
>> Also, do I need a* const *because then* fe_values.reinit(cell) *will be
>> discarded?
>>
>> Thanks
>> Wasim
>>
>> On Sun, Jan 8, 2023 at 3:44 PM Peter Munch <peterr...@gmail.com> wrote:
>>
>>> Change:
>>>
>>> > float H_plus(Vector<double> solution_elastic, const auto cell,const
>>> unsigned int q_point,
>>> * FEValues<2> fe_values_damage*);
>>>
>>> to:
>>>
>>> > float H_plus(Vector<double> solution_elastic, const auto cell,const
>>> unsigned int q_point,
>>> * const FEValues<2> & fe_values_damage*);
>>>
>>> Peter
>>> On Sunday, 8 January 2023 at 11:11:34 UTC+1 ce21...@smail.iitm.ac.in
>>> wrote:
>>>
>>>> Thank you, Prof. Munch.
>>>>
>>>> I tried to pass the FEValues object as a parameter to the H_plus
>>>> function as follows:
>>>>
>>>> I changed the function declaration from:
>>>> float H_plus(Vector<double> solution_elastic, const auto cell,const
>>>> unsigned int q_point);
>>>> to this:
>>>> float H_plus(Vector<double> solution_elastic, const auto cell,const
>>>> unsigned int q_point,
>>>> * FEValues<2> fe_values_damage*);
>>>>
>>>> I made the same change in the function definition also.
>>>>
>>>> I also changed the call to my function from:
>>>> H_plus(solution_elastic,cell,q_index)
>>>> to this: H_plus(solution_elastic,cell,q_index,*fe_values_damage*)
>>>>
>>>> But I am getting the following error:
>>>> error: use of deleted function ‘dealii::FEValues<2>::FEValues(const
>>>> dealii::FEValues<2>&)’
>>>> 761 | { float H_call =
>>>> H_plus(solution_elastic,cell,q_index,fe_values_damage);
>>>>
>>>>
>>>> On Sunday, January 8, 2023 at 2:32:24 PM UTC+5:30 peterr...@gmail.com
>>>> wrote:
>>>>
>>>>> You are creating a new instance of FEValues at each quadrature point.
>>>>> This is a very expensive operation, since there all the shape functions
>>>>> are
>>>>> evaluated. Try to reuse that by passing it to the function as a parameter.
>>>>>
>>>>> Hope that helps!
>>>>> Peter
>>>>>
>>>>> On Sunday, 8 January 2023 at 09:58:41 UTC+1 ce21...@smail.iitm.ac.in
>>>>> wrote:
>>>>>
>>>>>> Following is my H_plus function:
>>>>>>
>>>>>> float PhaseField::H_plus(Vector<double> solution_elastic
>>>>>> , const auto cell,const unsigned int q_point)
>>>>>> {
>>>>>> QGauss<2> quadrature_formula_damage(fe_damage.degree + 1);
>>>>>>
>>>>>> FEValues<2> fe_values_damage(fe_damage,
>>>>>> quadrature_formula_damage,
>>>>>> update_gradients |
>>>>>> update_quadrature_points );
>>>>>>
>>>>>> fe_values_damage.reinit(cell);
>>>>>>
>>>>>> int node = 0;
>>>>>>
>>>>>> /* Initialising all strains as zero */
>>>>>> float e_xx = 0.000;
>>>>>> float e_yy = 0.000;
>>>>>> float e_xy = 0.000;
>>>>>>
>>>>>> /*calculating strains*/
>>>>>> for (const auto vertex : cell->vertex_indices())
>>>>>> {
>>>>>> int a = (cell->vertex_dof_index(vertex, 0));
>>>>>> e_xx = e_xx +
>>>>>> solution_elastic[a*2]*fe_values_damage.shape_grad(node, q_point)[0];
>>>>>> e_yy = e_yy +
>>>>>> solution_elastic[a*2+1]*fe_values_damage.shape_grad(node, q_point)[1];
>>>>>> e_xy = e_xy
>>>>>> +0.5*(solution_elastic[a*2]*fe_values_damage.shape_grad(node, q_point)[1]
>>>>>>
>>>>>>
>>>>>> +solution_elastic[a*2+1]*fe_values_damage.shape_grad(node, q_point)[0]);
>>>>>> node = node +1;
>>>>>> }
>>>>>>
>>>>>> const auto &x_q = fe_values_damage.quadrature_point(q_point);
>>>>>> float H_plus_val;
>>>>>>
>>>>>> /*This is the actual quantity I want to evaluate for each quadrature
>>>>>> point*/
>>>>>> H_plus_val = 0.5*lambda(x_q)*(pow((e_xx+e_yy),2))
>>>>>> +
>>>>>> mu(x_q)*((pow(e_xx,2))+2*(pow(e_xy,2)) + (pow(e_yy,2)));
>>>>>>
>>>>>> return H_plus_val;
>>>>>> }
>>>>>>
>>>>>> H_plus function is called in assemble_damage for each quadrature point
>>>>>>
>>>>>> I am also attaching some lines of code from assemble_damage where
>>>>>> H_plus is being called
>>>>>>
>>>>>> for (const auto &cell : dof_handler_damage.active_cell_iterators())
>>>>>> {
>>>>>>
>>>>>> fe_values_damage.reinit(cell);
>>>>>> cell_matrix_damage = 0;
>>>>>> cell_rhs_damage = 0;
>>>>>>
>>>>>>
>>>>>> float gc = 0.000027; //energy release rate
>>>>>> float l = 0.015;
>>>>>> float H;
>>>>>>
>>>>>> for (const unsigned int q_index :
>>>>>> fe_values_damage.quadrature_point_indices())
>>>>>> { float H_call = H_plus(solution_elastic,cell,q_index);
>>>>>>
>>>>>>
>>>>>> if (H_call > H_vector[4*cell_number + q_index])
>>>>>> {
>>>>>> H = H_call;
>>>>>> }
>>>>>> else
>>>>>> {
>>>>>> H = H_vector[4*cell_number + q_index];
>>>>>> }
>>>>>> H_vector_new[4*cell_number + q_index] = H;
>>>>>> for (const unsigned int i : fe_values_damage.dof_indices())
>>>>>> {
>>>>>>
>>>>>>
>>>>>> for (const unsigned int j : fe_values_damage.dof_indices())
>>>>>> {
>>>>>>
>>>>>> const auto &x_q =
>>>>>> fe_values_damage.quadrature_point(q_index);
>>>>>>
>>>>>>
>>>>>> cell_matrix_damage(i, j) +=
>>>>>> // contribution to stiffness from -laplace u term
>>>>>>
>>>>>> Conductivity_damage(x_q)*fe_values_damage.shape_grad(i, q_index) * //
>>>>>> kappa*grad phi_i(x_q)
>>>>>> fe_values_damage.shape_grad(j, q_index) * // grad
>>>>>> phi_j(x_q)
>>>>>> fe_values_damage.JxW(q_index) // dx
>>>>>> +
>>>>>> // Contribution to stiffness from u term
>>>>>>
>>>>>>
>>>>>> ((1+(2*l*H)/gc)*(1/pow(l,2))*fe_values_damage.shape_value(i, q_index) *
>>>>>> // phi_i(x_q)
>>>>>> fe_values_damage.shape_value(j, q_index)
>>>>>> * // phi_j(x_q)
>>>>>> fe_values_damage.JxW(q_index));
>>>>>> // dx
>>>>>> }
>>>>>>
>>>>>> cell_rhs_damage(i) += (fe_values_damage.shape_value(i,
>>>>>> q_index) * // phi_i(x_q)
>>>>>> (2/(l*gc))* H*
>>>>>> fe_values_damage.JxW(q_index)); // dx
>>>>>> }
>>>>>> }
>>>>>>
>>>>>> /*Adding the local k and local f to global k and global f*/
>>>>>> cell->get_dof_indices(local_dof_indices);
>>>>>> for (const unsigned int i : fe_values_damage.dof_indices())
>>>>>> {
>>>>>> for (const unsigned int j : fe_values_damage.dof_indices())
>>>>>> system_matrix_damage.add(local_dof_indices[i],
>>>>>> local_dof_indices[j],
>>>>>> cell_matrix_damage(i, j));
>>>>>>
>>>>>> system_rhs_damage(local_dof_indices[i]) += cell_rhs_damage(i);
>>>>>> }
>>>>>> cell_number = cell_number + 1;
>>>>>>
>>>>>> }
>>>>>>
>>>>>> Thanks and regards
>>>>>> Wasim
>>>>>>
>>>>>> On Saturday, January 7, 2023 at 11:59:49 PM UTC+5:30
>>>>>> blais...@gmail.com wrote:
>>>>>>
>>>>>>> There might be many things that can be done to improve the speed of
>>>>>>> this function. You can ask yourselve the following question as guidance:
>>>>>>> - Does the function allocate memory?
>>>>>>> - Could it be inlined?
>>>>>>> - Are you calling the function inside the DOF loops or inside the
>>>>>>> quadrature loop?
>>>>>>>
>>>>>>> Then I would time the function to measure if this is actually the
>>>>>>> real culprit or if it could be something else.
>>>>>>> If you copy/paste the content of your assembly code and the
>>>>>>> function, I would be glad to give it a look (and I am sure others here
>>>>>>> will
>>>>>>> help you too).
>>>>>>>
>>>>>>>
>>>>>>> On Saturday, January 7, 2023 at 12:02:13 a.m. UTC-5
>>>>>>> ce21...@smail.iitm.ac.in wrote:
>>>>>>>
>>>>>>>> Sorry for the confusion. I think I made a mistake while writing the
>>>>>>>> first email.
>>>>>>>>
>>>>>>>> H_plus is being called in Assemble_damage and not assemble_elastic.
>>>>>>>> It uses elastic solution, cell and gauss point to evaluate strain at a
>>>>>>>> gauss point. Then some quantity is evaluated based on the strain.
>>>>>>>>
>>>>>>>> Similarly I have another function damage_gauss which is being
>>>>>>>> called in assemble_elastic that evaluates damage at a gauss point
>>>>>>>> using the
>>>>>>>> damage solution, cell and gauss point.
>>>>>>>>
>>>>>>>> Wasim Niyaz
>>>>>>>> Research scholar
>>>>>>>> CE Dept.
>>>>>>>> IITM
>>>>>>>>
>>>>>>>> On Sat, 7 Jan, 2023, 10:15 am Wasim Niyaz Munshi ce21d400, <
>>>>>>>> ce21...@smail.iitm.ac.in> wrote:
>>>>>>>>
>>>>>>>>> I use it to evaluate strain at Gauss points. Then, i evaluate some
>>>>>>>>> quantity which is a function of this strain.
>>>>>>>>>
>>>>>>>>> Wasim Niyaz
>>>>>>>>> Research scholar
>>>>>>>>> CE Dept.
>>>>>>>>> IITM
>>>>>>>>>
>>>>>>>>> On Sat, 7 Jan, 2023, 3:09 am Wolfgang Bangerth, <
>>>>>>>>> bang...@colostate.edu> wrote:
>>>>>>>>>
>>>>>>>>>> On 1/6/23 13:53, Wasim Niyaz Munshi ce21d400 wrote:
>>>>>>>>>> > I am using 65536 elements. For step-8 the assembly takes very
>>>>>>>>>> less time
>>>>>>>>>> > (around 0.15second) while for my assemble_elastic, it takes
>>>>>>>>>> around 5 seconds.
>>>>>>>>>> > The only difference between my assemble_elastic function and
>>>>>>>>>> the assemble
>>>>>>>>>> > function of step-8 is that for each Gauss point, I
>>>>>>>>>> additionally call a
>>>>>>>>>> > function(H_plus) that takes the laplace solution, the current
>>>>>>>>>> cell and Gauss
>>>>>>>>>> > point as input and evaluates some quantity using this
>>>>>>>>>> information.
>>>>>>>>>> > The H_plus function is called 4*65536 times but the function is
>>>>>>>>>> very simple.
>>>>>>>>>> > My doubt is whether such a huge increase in cost (from 0.15 sec
>>>>>>>>>> to 5 sec) is
>>>>>>>>>> > expected for this problem or is there something that I am doing
>>>>>>>>>> that is
>>>>>>>>>> > increasing the cost so much.
>>>>>>>>>>
>>>>>>>>>> Wasim, the question is what you do with "the laplace solution,
>>>>>>>>>> the current
>>>>>>>>>> cell and Gauss point". If you show us what H_plus does, we may be
>>>>>>>>>> able to advise.
>>>>>>>>>>
>>>>>>>>>> Best
>>>>>>>>>> W.
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>>
>>>>>>>>>> ------------------------------------------------------------------------
>>>>>>>>>> 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 a topic
>>>>>>>>>> in the Google Groups "deal.II User Group" group.
>>>>>>>>>> To unsubscribe from this topic, visit
>>>>>>>>>> https://groups.google.com/d/topic/dealii/-6ndTW_k5fQ/unsubscribe.
>>>>>>>>>> To unsubscribe from this group and all its topics, send an email
>>>>>>>>>> to dealii+un...@googlegroups.com.
>>>>>>>>>> To view this discussion on the web visit
>>>>>>>>>> https://groups.google.com/d/msgid/dealii/895a079c-2b85-b14f-94ee-b4b78336884d%40colostate.edu
>>>>>>>>>> .
>>>>>>>>>>
>>>>>>>>> --
>>> 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 a topic in the
>>> Google Groups "deal.II User Group" group.
>>> To unsubscribe from this topic, visit
>>> https://groups.google.com/d/topic/dealii/-6ndTW_k5fQ/unsubscribe.
>>> To unsubscribe from this group and all its topics, send an email to
>>> dealii+un...@googlegroups.com.
>>>
>> To view this discussion on the web visit
>>> https://groups.google.com/d/msgid/dealii/deb31ad3-7a15-438b-bba4-bdb79d17a18bn%40googlegroups.com
>>> <https://groups.google.com/d/msgid/dealii/deb31ad3-7a15-438b-bba4-bdb79d17a18bn%40googlegroups.com?utm_medium=email&utm_source=footer>
>>> .
>>>
>> --
> 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/7e840958-f3af-49a1-9b53-4bea537a13d8n%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/7e840958-f3af-49a1-9b53-4bea537a13d8n%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
--
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/CAMvdQMCu402WOZMxxN7hdvOpdNLjvxTxEg34F_%3Dm7Us%3Dt%3D3HcQ%40mail.gmail.com.
