Assemble system slowing down the program

[Wasim Niyaz Munshi ce21d400]

Hello everyone.

I am trying to solve elasticity and laplace equations one after the other. The elasticity solution feeds as input to the laplace equation (stiffness matrix of laplace equation depends on elasticity solution) and vice versa.

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.

I am running the code in release mode.

Thanks and regards

Wasim

[Wolfgang Bangerth]

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/

[Wasim Niyaz Munshi ce21d400]

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

--
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.

[Wasim Niyaz Munshi ce21d400]

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

[Wolfgang Bangerth]

On 1/6/23 22:01, Wasim Niyaz Munshi ce21d400 wrote:
>
> 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 -- that is not helpful. You complain that the function takes a lot of
time to run, but your description contains nothing that would anyone help
understanding *why* that might be so. For all we know "then some quantity is
evaluated based on the strain" could be an operation that solves a
differential equation in itself.

Show us the function, or better even: Do some timings of that function
yourself to figure out which line or code block within that function is
expensive. Then we can think of improvements.

[blais...@gmail.com]

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).

[Wasim Niyaz Munshi ce21d400]

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

[Peter Munch]

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

[Wasim Niyaz Munshi ce21d400]

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);

[Peter Munch]

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

[Wasim Niyaz Munshi ce21d400]

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

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/deb31ad3-7a15-438b-bba4-bdb79d17a18bn%40googlegroups.com.

[Wasim Niyaz Munshi ce21d400]

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

[Jean-Paul Pelteret]

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.

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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/7e840958-f3af-49a1-9b53-4bea537a13d8n%40googlegroups.com.

[Wasim Niyaz Munshi ce21d400]

Thank you for pointing that out.

After passing my solution vector as reference the assembly time has reduced from 2 seconds to 0.2 seconds.

Thanks and regards

Wasim Niyaz
Research scholar
CE Dept.
IITM

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CAMvdQMCu402WOZMxxN7hdvOpdNLjvxTxEg34F_%3Dm7Us%3Dt%3D3HcQ%40mail.gmail.com.