Personal tools
You are here: Home / cmgui / Wiki / Convergence tests
Navigation
Log in


Forgot your password?
 

Convergence tests

PAGE IN DEVLOPMENT

To test the validity of your solution it is often a good idea to perform a convergence test.

cmiss is not always particularly friendly when it comes to performing convergence tests and to date I find that the best option is to write a set of perl scripts and shell scripts to do most of the work for you.

Below I will detail the files that I created and commands that I used to perform a convergence test using coupled electro mechanics. This simulation involves cellml, mono domain solve, finite deformation and time. Although similar techniques could be used for other 'lung' based solutions.

DISCLAIMER This is probably not the best method for performing a convergance test but it works (ask Karl if you need advice form some one smarter).

When solving cmiss problems you will find yourself with multiple versions of .com, .pimate .exnode ..... etc

To avoid a wall of confusing filenames I tend to set up problems with the following directories: geom (for ipnode, ipelem, ipelfb ipfibr) enviro (for ipbase, ippara) mech (relevant ipmate, ipsolve, ipequa) cell (relevant ipcell, .cml, ipequa, ipmate, ipmatc, ...) init (for all initial conditions becomes important for mechanics)

for the problem output I create a directory: output/

containing elem node gauss image strain tension E principal_strains

these allow you to store data as well as generate and store additional fields calculated using cmgui.

Having set up your problem directories you require a '.com' file that you wish to perform you convergence test on.

First we will deal with mesh refinement.

For a mesh convergence problem you may have a number of elements in the 3 RC directions say a,b and c for the height, width and length of a cube or the height, circumference and radius of a cylinder.

We will set these as: $NE_a $NE_b $NE_c

We now wish to define a set of perl scripts that will recalculate your ipnode, ipelfb, ipelem, ipfibr, ipgrid, ipinit files for a given combination of $NE_a $NE_b $NE_c. these will take the form:

sub create_mesh {

$NE_a = shift(@_); $NE_b = shift(@_); $NE_c = shift(@_);

...

}

and we will assume saved in the file mesh_sub.pm. This routine will be called from your '.com' file by adding : use mesh_sub at the top of the '.com' and called in the file by: create_mesh($NE_a,$NE_b,$NE_c);

This exercise can be repeated for grid_sub.pm and init_sub.pm (feel free to choose your own more informative names)

At this stage you will have a '.com' file which will solve for any given combination of $NE_a $NE_b $NE_c.

Now we wish to set up a script that will run our problem a number of times changing one or some of $NE_a $NE_b $NE_c. For the sake of this example if you wished to change $NE_c you would replace the definition of $NE_c in you '.com' file with: $NE_c = a_replacebale_number (do not have any other variable called a_replacebale_number in your com file)

As cmiss does not always like solving a problem with different meshses, cell models etc twice it is often best to restart cmiss for each solve. To do this add a 'q' at the bottom of your '.com' file and create a script.

your script will look something like:

i=1 while $i -lt 16 ; do echo "solving for $i elements"; file_out="solve_ref_$i.com" sed s/a_replacable_number/$i/ <solve_temp.com> $file_out /home/users/physiome/san/mycmiss/cm/bin/x86_64-linux/cm-mt $file_out rm -f $file_out i=$((i*2))

this script runs for i = 1,2,4,8,16 and your '.co,' template is in solve_temp.com. sed swaps a_replacable_number in the '.com' file with the value of $i and you perl scripts should then generate the mesh and grid required.

So now we have a means to run a set of problems with increasing mesh refinement. Now we need to compare the results.

A challenge with this is that you wish to get data (say the strain, tension etc) at points defined by xi coordinates within an element. You could always just use node points but may be you want MORE.

So given that you have exported you deformed and unreformed node and elements points and if you are using a cell model to calc active tension then that you have used this to out put an active tension field.

you are able to calculate the large strain and principal strain field using example a/large_strain from the cmgui example tree.

you now need to define a set of xi points that you wish to obtain data from. Once again I generated a perl script to do this so that I could calculate the equivalent xi points for each refinement.

For a set of xi points space evenly along a rod at xi = 0,0.25,0.5,0.75 and 1.0 in one elem the file to create the xi points looked like:

$Ne = 1; $Nn = 4000000;

$filename = "output/xi_points_$NE_c.xi"; open(OUTFID, ">$filename") || die "Couldn't open $filename"; print OUTFID " Group name: xitestn"; print OUTFID " #Fields=1n"; print OUTFID " 1) element_xi, field, element_xi, #Components=1n"; print OUTFID " 1. Value index=1, #Derivatives=0n"; print OUTFID " Node: $Nnn"; print OUTFID " E $Ne 3 1.00000000000000 1.00000000000000 0.00000000000000n"; $delta_xi = $NE_c/4; $xi = 0; for ($points = 1; $points < 5; $points++) { $Nn = $Nn+1; $xi =$xi+$delta_xi; while($xi>1) { $xi = $xi-1; $Ne = $Ne +4; }; print OUTFID " Node: $Nnn"; print OUTFID " E $Ne 3 1.00000000000000 1.00000000000000 $xin"; };

Having generated a set of xi points at for each elem you can then use cmgui to evaluate the desired fields at each of these points.

It is important to remember that xi points are considered embedded fields in cmgui so you must convert your strain/tension fields into embedded fields before calculating the strain/tension values at the defined xi points

Read in the xi points: gfx read n xi_points_$NE_c.xi;

Define a set of fields gfx define field E_xi finite_element num 9 # create field E_xi with 9 components gfx modify node define E_xi group xitest # add field to xitest group

gfx define field pot_xi finite_element num 1 #pot (potential) holds the tension from cell gfx modify node define pot_xi group xitest

gfx define field p_strain_xi finite_element #(principal strains do not need to define num as defualt is 3) gfx modify node define p_strain_xi group xitest

for the calcualted fields E, principal_strains and potential (tension from cell model) you calculate the embeded fields: element_xi_E, element_xi_principal_strains and element_xi_potential.

gfx define field element_xi_E embedded field E element_xi element_xi; gfx define field element_xi_principal_strains embedded field principal_strains element_xi element_xi; gfx define field element_xi_potential embedded field potential element_xi element_xi;

you then evalute the each of our newly embedded fields at each of the xi points defined in xi test at the three fields (E_xi pot_xi p_strains) created at xitest points.

gfx evaluate ngroup xitest source element_xi_E destination E_xi gfx evaluate ngroup xitest source element_xi_principal_strains destination p_strain_xi gfx evaluate ngroup xitest source element_xi_potential destination pot_xi

We can then export these to a file $NAM = "output/principal_strains/p_strains_$version.$time.x" gfx write nodes fields p_strain_xi group xitest $NAM

$NAM = "output/tension/tension_$version.$time.x" gfx write nodes fields pot_xi group xitest $NAM

NAM = "output/E/E_$version.$time.x" gfx write nodes fields E_xi group xitest $NAM

whwre $version.$time defiens the tiem and version the level of refinement.

A few more perl scripts can strip out the values from the '.x' files and you can plot and compare these using gnuplot.