Error attempting Bone Remodelling Analysis

[laure...@gmail.com]

I am unsure why I am encountering an error when running the bone remodelling analysis.

I have copied the code from bitbucket and adapted the parameters to match the Kuhl and Steinmann paper.

However the analysis will not run past the first time step.

I have attached an image of my code and the error message for reference.

Thanks

[Lukasz Kaczmraczyk]

Hello Lauren,

Can you attach your mesh file, history file and command line how you run the code.

Regards,

Lukasz

[Karol Lewandowski]

Your mesh size is suspeciously large, I can see from the mass. Check your units

[laurenboag1]

Hello, I have attached the requested below.

Thanks,

Lauren

mpirun -np 4 ./bone_adaptation -my_file bar_1D_on_4parts.h5m -my_order 2 -my_output_prt 1 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package mumps -ksp_atol 1e-12 -ksp_rtol 1e-12 -snes_monitor -snes_type newtonls -snes_linesearch_type basic -snes_max_it 1000 -snes_atol 1e-6 -snes_rtol 1e-8 -ts_type beuler -ts_max_snes_failures 1 -ts_dt 0.1 -ts_final_time 8 -my_load_history 1D_load_history.txt -mass_postproc -ts_monitor -young_modulus 1 -poisson_ratio 0 -c 1 -rho_ref 1 -psi_ref 1 -m 3 -n 2 | tee log

[laurenboag1]

My bar is 1 x 1 x 10, and the mesh size is 8.

Should I change the units of the cubit file?

Thanks,

Lauren

[Karol Lewandowski]

Yes, change the dimensions. And your load_history should look like this:

0.0             0.0

0.1             0.5

2.0             0.5

2.1             1.0

4.0             1.0

4.1             1.5

6.0             1.5

6.1             2.0
8.0                           2.0

[laurenboag1]

I have now reduced my bar to 1 x 1 x 5 and corrected my load history file but I am still having problems.

The analysis does not run past the first time step when the Young's modulus is 1.

I have tried using E = 1000, but the analysis crashes around time = 3s.

[Karol Lewandowski]

It's crashing because the elements are getting squashed too much. The original example is in 1D and the length of the bar is 1 (not 5 or 10). You can also change the direction of the loading. 

[laurenboag1]

Hi Karol,

I now have a new error message relating to file accessibility.

[Karol Lewandowski]

use mofem_part tool instead of cp

[laurenboag1]

I have used the mofem part command and I'm still getting the same error.

[Karol Lewandowski]

Now you used cubit file for some reason. Use the command like this: 

../tools/mofem_part -my_file 1D_bar.cub  -my_nparts 1 && mpirun -np 1 ./bone_adaptation -my_file out.h5m -my_order 2 -my_output_prt 1 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package mumps -ksp_atol 1e-12 -ksp_rtol 1e-12 -snes_monitor -snes_type newtonls -snes_linesearch_type basic -snes_max_it 10 -snes_atol 1e-8 -snes_rtol 1e-8 -ts_type beuler -ts_dt 0.1 -ts_final_time 8 -ts_max_snes_failures -1 -my_load_history load_history.in  -young_modulus 1 -poisson_ratio 0 -c 1 -rho_ref 1 -psi_ref 1 -m 3 -n 2  -mass_postproc -ts_monitor | tee log

[laurenboag1]

I have managed to get the MoFEM analysis operating.

The graphs seem to follow the same trend as the Kuhl and Steinmann paper however the numbers are not a match.

Any suggestions as to why?

Thanks,

Lauren

[Lukasz Kaczmraczyk]

Lauren,

You asking why displacements are so large? Psi and rho do nit seem to be wrong?

You used cube from 1D_bar.cub, if I understand correctly. Can you explain how you applied force and how you constrain displacements? 

Do you have any kinematic/displacements constrains?

Kind regards,

Lukasz

[Lukasz Kaczmraczyk]

Lauren,

I've changed boundary conditions, such that rigid motion is not possible. And changed forces to pressures, and apply pressure to make tension not compression. See attached.

Regards,

Lukasz

[Karol Lewandowski]

Lauren,

here is the mesh that will work. 

Also, remember if you want to compare the results directly with the paper, some of the quantities that are presented there have to be calculated from the results, like this: 

[laurenboag1]

Hi Lukasz,

I have tried to run the analysis using the model you attached, using the code below.

mpirun -np 4 ./bone_adaptation -my_file out.h5m -my_order 2 -my_output_prt 1 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package mumps -ksp_atol 1e-12 -ksp_rtol 1e-12 -snes_monitor -snes_type newtonls -snes_linesearch_type basic -snes_max_it 1000 -snes_atol 1e-6 -snes_rtol 1e-8 -ts_type beuler -ts_max_snes_failures 1 -ts_monitor -ts_dt 0.1 -ts_final_time 8 -my_load_history 1D_load_history.txt -mass_postproc -ts_monitor -young_modulus 1 -poisson_ratio 0 -c 1 -rho_ref 1 -psi_ref 1 -m 3 -n 2 | tee log

The code is crashing at time = 2.9 or 3.0 seconds. 

If I change the code such that  ' ts_max_snes_failures -1', the code runs, however at from time 2.9, the time step begins to get much smaller than the prescribed 0.1s.

[Karol Lewandowski]

Try my mesh. 

[laurenboag1]

Hi Karol,

I have tried your mesh, and the mofem analysis does not run.

I do not receive an error message, but nothing happens and I have to terminate the job (see attached image).

Kind regards,

Lauren

[Karol Lewandowski]

Try with one partition. 

[laurenboag1]

This now runs with one partition. 

What effect does the number of partitions have on the solutions?

[Karol Lewandowski]

I general it doesn't. This is a special case, because the problem has only 12 elements and you tried to partition it to 4 parts. For such a small problem there is no gain from running on multiple cores and as it turns out it can cause some problems.