libMesh SVD convergence failure

[Illy Perl]

Hello, 

I'm running a flapping wings simulation using IBFE, in IBAMR v0.8.0

After quit a few timesteps, I get the following error:

"...

At beginning of timestep # 115418

Simulation time is 1.15418

The algorithm for computing the SVD failed to converge!

Stack frames: 11

0: libMesh::print_trace(std::ostream&)

1: libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*)

2: libMesh::DenseMatrix<double>::_svd_solve_lapack(libMesh::DenseVector<double> const&, libMesh::DenseVector<double>&, double) const

3: ./main3d() [0x938765]

4: ./main3d() [0x9265ab]

5: ./main3d() [0x660325]

6: ./main3d() [0x4f82c3]

7: ./main3d() [0x87f9c4]

8: ./main3d() [0x429c95]

9: __libc_start_main

10: ./main3d() [0x43e69f]

[0] ../LIBMESH/src/numerics/dense_matrix_blas_lapack.C, line 660, compiled Feb 10 2021 at 07:22:54

Abort(1) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0

finished IBAMR simulation

..."

I've tried changing the solid mesh, and changing the kinematics. I keep getting the same error,  though at different times.  

I couldn't find any info of this kind of libMesh error online. 

I would greatly appreciate any kind of help to understand and/or solve that error. 

Thanks, 

Illy

[Boyce Griffith]

We have switched to using Eigen for this solve in IBAMR’s master branch. Can you see if the error persists if you use IBAMR master?

Thanks!

— Boyce

I would greatly appreciate any kind of help to understand and/or solve that error. 

Thanks, 

Illy

--
You received this message because you are subscribed to the Google Groups "IBAMR Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ibamr-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/71c7ff5f-47dd-4528-95d5-c56ae831b902n%40googlegroups.com.

[Illy Perl]

Thanks Boyce! for the quick response

I'll do that and update. 

Might take a couple of days because the error occurs after several days of computation.

[Illy Perl]

Hi again, 

it took me quite a long time to get the system admin to help me with updating to the master version, but we eventually did it.

Now something very weird happened:

when I run the same code using the master Ibamr version, it doesn't crush, but instead, right when the older version crushes, all the forces kind of 'turn off' and the mesh just deforms with the flow field. 

It starts exactly at the moment where the simulation crashed with the previous version. 

Have you seen this kind of behavior? or perhaps you can direct me what king of things I should be looking for, in the code or elsewhere? 

I attach 2 snapshots of the mesh, before and after it starts deforming.

 

Thanks again, 

Illy 

[Wells, David]

Hi Illy,

This probably indicates that something unexpected happened inside the finite element solver, like an inverted element, which caused those linear solvers to no longer converge. If the solvers don't converge then the forces will remain at default values (zero).

I added a more careful check for this behavior last week in

https://github.com/IBAMR/IBAMR/pull/1276

but, in principle, the SVD should always converge as long as the matrix contains valid values. If it does not contain valid values that means that the norm of either a RHS or solution vector is infinite or NaN - i.e., the SVD can fail, but only if something else failed first.

A few questions:

1) Do you have a log file from this simulation? If so, would you please either send it to me or attach it here?

2) How long does it take the code to reach this point? Do you have restart files nearby? If so, it would help if you added

-ksp_monitor

to the command line arguments so that we can see how well the finite element linear solvers are performing.

3) If it is a finite element problem, I recommend adding some code to check for element inversion - if you are computing a dilational stress (or have det(FF) computed for some other reason) then it might be useful to check if det(FF) gets very small (like less than 0.01). You'll have to write some MPI code to do the min reduction - I can throw something together if you don't know how to do it.

Best,

Davi

---

From: ibamr...@googlegroups.com <ibamr...@googlegroups.com> on behalf of Illy Perl <illy...@mail.huji.ac.il>
Sent: Sunday, July 4, 2021 8:57 AM
To: IBAMR Users <ibamr...@googlegroups.com>
Subject: Re: [ibamr-users] libMesh SVD convergence failure

 

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/fd7fbf14-3906-4360-9d82-f3e269e40da3n%40googlegroups.com.

[Illy Perl]

Hi David,

First, I really appreciate the help! Thank you.

I'll try to answer your questions the best I can:

1) I do have a log file, hope that is the one you are looking for. I'll attach a google drive link to a folder containing everything, is that OK?

2) It takes quite a long time for the simulation to reach the bug. About 135k steps, which took several days to compute. 

I do have the restart files, I'm not sure where and when to use the CLI argument?

3) you guessed correctly. It is FE, and I do compute dilatational stress. 

and I do not know how to write the MPI code. If it's not too much trouble I'd appreciate some help with it. 

If it is a problem of element inversion, do you have a guess why it is occurring? or how can I try to avoid it? Does it have something to do with the mesh I'm starting with? 

Much thanks! 

Illy 

[Wells, David]

Hi Illy,

I'm happy to help!

The log file isn't that useful on its own, but you are already printing out min and max values of Jacobians (which are useful) - they just aren't getting updated, so we just see

Max Jacobian = 1

Min Jacobian = 1

at every time step. All you need to add is, in the dilational stress function (circa line 278),

J_maximum = std::max(J_maximum, FF.det());

J_minimum = std::min(J_minimum, FF.det());

and then the values printed will be the ones we want to see.

You can use restart files by running your program as

./main3d ./input3d ./restart_IB3d 200000

Here ./restart_IB3d is the directory containing all the restart data and the highest number I see is 200000. If you want something lower do the natural thing (i.e., start from checkpoint 100000 by writing that number instead).

In general, you shouldn't use restart files when you change things in the solver, but since we are just changing some logging it should be okay.

Element inversion in IBFE is tricky - Boyce knows more about this than I do. It typically happens if you don't have strong enough penalties enforcing incompressibility. If your mesh is very low quality then it might be possible for it to get tangledon its own, but the pictures you showed look pretty good so I don't think it's a flaw fundamentally attributable to the grids.

To be clear - element inversion is just my best guess. It still could be something else - we won't know unless we see the Jacobians go negative (or get very small).

Best,

David 

---

From: ibamr...@googlegroups.com <ibamr...@googlegroups.com> on behalf of Illy Perl <illy...@mail.huji.ac.il>

Sent: Monday, July 5, 2021 10:24 AM

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/a4062103-7249-48ee-8def-00de12703bc3n%40googlegroups.com.

[Boyce Griffith]

On Jul 5, 2021, at 10:24 AM, Illy Perl <illy...@mail.huji.ac.il> wrote:

Hi David,

First, I really appreciate the help! Thank you.

I'll try to answer your questions the best I can:

1) I do have a log file, hope that is the one you are looking for. I'll attach a google drive link to a folder containing everything, is that OK?

2) It takes quite a long time for the simulation to reach the bug. About 135k steps, which took several days to compute. 

I do have the restart files, I'm not sure where and when to use the CLI argument?

3) you guessed correctly. It is FE, and I do compute dilatational stress. 

and I do not know how to write the MPI code. If it's not too much trouble I'd appreciate some help with it. 

If it is a problem of element inversion, do you have a guess why it is occurring? or how can I try to avoid it? Does it have something to do with the mesh I'm starting with? 

Depending on the compiler and compiler flags, the solvers might not catch “inf” or “nan” values that may be computed if something goes wrong with the structural model. My experience is that if this happens, then the structural forces can wind up being set to 0. And if that happens, the structure “drifts away” with the fluid. This seems to be what you are describing.

Can you give some additional details on the state of the model when this happens? If your elastic energy or stress include terms with J^{-1} or log(J), then bad things will happen if you experience any element inversions.

If you are getting element inversion, then one thing that you can do is to use an elastic energy that includes a volumetric term that penalizes compressive deformations. The structural model will wind up looking like a nearly incompressible solid mechanics formulation. This approach is discussed in this paper:

https://doi.org/10.1016/j.cma.2020.112978

https://arxiv.org/abs/1811.06620

and used in, e.g., this paper (with a different volumetric energy):

https://doi.org/10.1007/s10439-020-02466-4

https://doi.org/10.31224/osf.io/8ys2k

One issue to keep in mind is that the continuum equations implemented in IBFEMethod correspond to an exactly incompressible structure. This would suggest that the structural mesh should not experience element inversions. However, the discretization of the continuum equations does not guarantee that the structural mesh only experiences incompressible deformations. In the paper that I link to above, we discuss controlling compressible deformations using what we call a numerical bulk modulus or a numerical Poisson ratio. Note that the physical bulk modulus for this formulation is not a well defined quantity because the material model is exactly incompressible — i.e., the physical Poisson ratio is 0.5.

Long story short — including a volumetric energy is one approach that we have found to avoid inversions — albeit potentially at the cost of smaller time step sizes if you wind up needing to use a very large penalty term (i.e. a large numerical bulk modulus). We are working on some new solvers that will treat at least the volumetric energy implicitly, which would allow for the time step sizes to be independent of the numerical bulk modulus, but these are still under development.

Another approach — not currently implemented in IBAMR — is to use divergence-free interpolation. (Although even with divergence-free interpolation, it is still possible for a discrete structural model to experience volume changes.) One approach is described here:

https://doi.org/10.1016/j.jcp.2017.06.041

https://arxiv.org/abs/1701.07169

— Boyce

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/a4062103-7249-48ee-8def-00de12703bc3n%40googlegroups.com.

[Illy Perl]

Thanks, David and Boyce!

OK, you are right. The Jacobian does get very small and eventually negative. So it is element inversion I guess. 

Boyce, thanks for the references. 

Are there any IBFE examples that show how to include the volumetric term in code? 

Thanks again, 

Illy 

[Boyce Griffith]

On Jul 28, 2021, at 7:28 AM, Illy Perl <illy...@mail.huji.ac.il> wrote:

Thanks, David and Boyce!

OK, you are right. The Jacobian does get very small and eventually negative. So it is element inversion I guess. 

Boyce, thanks for the references. 

Are there any IBFE examples that show how to include the volumetric term in code? 

There is not, but the code to write is pretty simple, and I can send a snippet. What kind of material model are you using currently?

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/668b6fd2-dd14-4b35-ae5e-8b7012d38684n%40googlegroups.com.

[Illy Perl]

Thanks, Boyce

The Material model I currently use is:

void PK1_dev_stress_function(

TensorValue<double>& PP,

const TensorValue<double>& FF,

    const libMesh::Point& /*X*/,

    const libMesh::Point& /*s*/,

    Elem* const /*elem*/,

    const std::vector<const std::vector<double>*>& /*var_data*/,

    const std::vector<const std::vector<VectorValue<double> >*>& /*grad_var_data*/,

    double /*time*/,

void* /*ctx*/)

{

double mu;

mu = mu_stem;

PP = (mu)*FF;

return;

}// PK1_dev_stress_function

void PK1_dil_stress_function(

TensorValue<double>& PP,

const TensorValue<double>& FF,

    const libMesh::Point& /*X*/,

    const libMesh::Point& /*s*/,

    Elem* const /*elem*/,

    const std::vector<const std::vector<double>*>& /*var_data*/,

    const std::vector<const std::vector<VectorValue<double> >*>& /*grad_var_data*/,

    double /*time*/,

void* /*ctx*/)

{

double mu;

double beta;

mu = mu_stem;

beta = beta_stem;

PP = (-1.0*(mu) + beta*log(FF.det()))*tensor_inverse_transpose(FF,NDIM);

J_maximum = std::max(J_maximum, FF.det());

J_minimum = std::min(J_minimum, FF.det());

return;

}// PK1_dil_stress_function

I thought that the PK1_dil_stress_function should enforce volume change, but I run tests with larger beta values and it didn't change the rate of change of the Jacobian. 

[Boyce Griffith]

OK, this is a sort of neo-Hookean model. The “dil” stress you were using is from code that was written before the time when we started reinforcing the incompressibility constraint using a volumetric energy.

Here is some updated code that may be more effective:

static double shear_mod, bulk_mod;

void

PK1_dev_stress_function(TensorValue<double>& PP,

                        const TensorValue<double>& FF,

                        const libMesh::Point& /*x*/,

                        const libMesh::Point& /*X*/,

                        Elem* const /*elem*/,

                        const vector<const vector<double>*>& /*var_data*/,

                        const vector<const vector<VectorValue<double> >*>& /*grad_var_data*/,

                        double /*time*/,

                        void* /*ctx*/)

{

    const RealTensor CC = FF.transpose() * FF;

    const double I1 = CC.tr();

    const double J = FF.det();

    const RealTensor dI1_bar_dFF = pow(J, -2.0 / 3.0) * (FF - I1 / 3.0 * tensor_inverse_transpose(FF, NDIM));

    PP = shear_mod * dI1_bar_dFF;

}

void

PK1_dil_stress_function(TensorValue<double>& PP,

                        const TensorValue<double>& FF,

                        const libMesh::Point& /*x*/,

                        const libMesh::Point& /*X*/,

                        Elem* const /*elem*/,

                        const vector<const vector<double>*>& /*var_data*/,

                        const vector<const vector<VectorValue<double> >*>& /*grad_var_data*/,

                        double /*time*/,

                        void* /*ctx*/)

{

    const double J = FF.det();

    PP = bulk_mod * J * log(J) * tensor_inverse_transpose(FF, NDIM);

}

Notice that dilational stresses will be generated by deformations that deviate from J = 1, and increasing the value of bulk_mod will increase the energetic penalization of such deviations. On the other hand, making bulk_mod larger and larger will eventually require a reduction in the time step size. In essence, you want to choose a value of bulk_mod that is as small as possible while being large enough to get your model to “work”.

I would recommend using what we call the numerical Poisson ratio to determine the bulk modulus from the shear modulus:

bulk_mod = 2.0 * shear_mod * (1.0 + poisson_ratio) / (3.0 * (1.0 - 2.0* poisson_ratio));

Depending on the loading conditions, you may be able to get away with using a modest value of the numerical Poisson ratio = nu ~ 0.4. In more extreme conditions, you might need to increase it. Note that the limit nu —> 0.5 corresponds to an exactly incompressible material response, yielding bulk_mod —> infinity. Values of nu close to 0.5 can result in an effect known as locking. What kinds of structural elements are you using? Locking is less of an issue if you are using second-order elements. It also can be avoided by using smaller values of the Poisson ratio. (This is another reason to try to use a bulk modulus that is as small as possible.)

By the way, in principle, the continuous equations that are discretized by in IBFEMethod correspond to an exactly incompressible material response. This means that, in the continuous equations, J = det(FF) is identically 1. However, exact incompressibility is lost in discretization. We have found that the approach outlined above can be effective in “reinforcing” the incompressibility in the discretized model. There is a discussion of this issue in this paper:

https://doi.org/10.1016/j.cma.2020.112978

https://arxiv.org/abs/1811.06620

Another thing to keep in mind is that mesh quality can impact this as well. “Sliver” elements (e.g., in three spatial dimensions, tetrahedral elements with vertices that are nearly coplanar) will make it harder to avoid large deviations from J ~ 1.

— Boyce

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/5db6a44a-76a1-4226-9062-56d044b09471n%40googlegroups.com.

[Illy Perl]

Thanks, Boyce!

I'm running the simulation with the new stress functions you provided for a couple of days now. I run 2 simulations with different bulk and sear moduli, corresponding to poison ratios of 0.4 & 0.45. 

I'm printing the Jacobians, calculated as David suggested above (J=FF.det()), and I get exactly the same min and max reading for both Poison ratios. It's also the same printout as for the previous simulation, the one with the old stress functions. 

So I assume  I'm doing something wrong. The Jacobian printouts should be at list somewhat different for all those different simulations, right? 

Perhaps I'm not registering the stress functions correctly, so the simulation doesn't even use them? 

In the main file, I have:

// Configure the IBFE solver.

        

IBFEMethod::LagBodyForceFcnData target1_force_data(target1_force_function);

ib_method_ops->registerLagBodyForceFunction(target1_force_data, 0);

        

IBFEMethod::LagBodyForceFcnData target2_force_data(target2_force_function);

ib_method_ops->registerLagBodyForceFunction(target2_force_data, 1);     

IBFEMethod::PK1StressFcnData PK1_dev_stress_data(PK1_dev_stress_function);

IBFEMethod::PK1StressFcnData PK1_dil_stress_data(PK1_dil_stress_function);

ib_method_ops->registerPK1StressFunction(PK1_dev_stress_data, 0);

ib_method_ops->registerPK1StressFunction(PK1_dev_stress_data, 1);

ib_method_ops->registerPK1StressFunction(PK1_dil_stress_data, 0);

ib_method_ops->registerPK1StressFunction(PK1_dil_stress_data, 1);       

        

ib_method_ops->initializeFEEquationSystems();

EquationSystems* mesh1_equation_systems = ib_method_ops->getFEDataManager(0)->getEquationSystems();

EquationSystems* mesh2_equation_systems = ib_method_ops->getFEDataManager(1)->getEquationSystems();

        

PK1_dev_stress_data.quad_order = Utility::string_to_enum<libMeshEnums::Order>(input_db->getStringWithDefault("PK1_DEV_QUAD_ORDER","SEVENTH"));

PK1_dil_stress_data.quad_order = Utility::string_to_enum<libMeshEnums::Order>(input_db->getStringWithDefault("PK1_DIL_QUAD_ORDER","SEVENTH"));

Or perhaps I have something laking in the input3d file? 

Maybe you have an idea what it might be? 

thanks, 

illy  

[Yang Zhang]

Dear Prof. Boyce,

Is there a model for calculating pure tensile stress in IBFE?

Thanks,

Yang

[Boyce Griffith]

On May 31, 2022, at 10:07 PM, Yang Zhang <yang.zha...@gmail.com> wrote:

Dear Prof. Boyce,

Is there a model for calculating pure tensile stress in IBFE?

What exactly do you have in mind?

If you mean something like elastic stresses that are active only in tension but not in compression, then if you can formulate the material response using anisotropic (pseudo)invariants, it is pretty easy to add a “tension/compression switch”. E.g. write the energy in term of

I4_star = min(1, I4)

with I4 = |lambda|^2 = square of the stretch in the preferred material direction.

Please let me know if this isn’t the kind of thing you are thinking of.

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/d9834cf0-46d8-4e89-87b4-bf50434b2b73n%40googlegroups.com.

[Yang Zhang]

Thanks, Professor. That's what I am thinking of.

Best,

Yang 

[Boyce Griffith]

On Jun 6, 2022, at 9:45 PM, Yang Zhang <yang.zha...@gmail.com> wrote:

Thanks, Professor. That's what I am thinking of.

Note that you also could use a superposition of multiple energies/stresses associated with different material directions.

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/16596822-9809-4aa2-83d4-6b02bd528c45n%40googlegroups.com.

[Matea Santiago]

Hello,

I have a quick question about this response. In the PK1_dil_stress_function you have  PP = bulk_mod * J * log(J) * tensor_inverse_transpose(FF, NDIM); In the paper you give, I only see PP=bulk_mod*log(J)FF^{-T} (no additional multiplication of J). I didn't see the model bulk_mod*J*log(J)FF^{-T} anywhere else in the paper. Is it a different model? If so, is it better/ more desirable for some reason?

Thank you,

Matea

[Boyce Griffith]

On Mar 27, 2023, at 8:27 PM, Matea Santiago <mateas...@math.arizona.edu> wrote:

Hello,

I have a quick question about this response. In the PK1_dil_stress_function you have  PP = bulk_mod * J * log(J) * tensor_inverse_transpose(FF, NDIM); In the paper you give, I only see PP=bulk_mod*log(J)FF^{-T} (no additional multiplication of J). I didn't see the model bulk_mod*J*log(J)FF^{-T} anywhere else in the paper. Is it a different model? If so, is it better/ more desirable for some reason?

There are different expressions that one could use for U(J). Sec. 4 of this paper discusses some choices: https://www.sciencedirect.com/science/article/abs/pii/S0045782516307563. We have not systematically compared different options, but have been mostly using the form listed in 3.4.6 in large models.

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/3ea18dc0-58d2-404f-a046-678acf62c446n%40googlegroups.com.

[Matea Santiago]

Thank you very much!