Time: tradeoff
Jump to navigation
Jump to search
Return to GFtbox hints and tips
Choosing a timestep that is short enough to be accurate and yet not take too long to compute
There is a tradeoff between speed and the accuracy with which the equation solver can solve the equations. This is feature of all numerical modelling. Continuous time is approximated by a series of short intervals or steps (dt). With steps that are too long the mesh will grow more than about 2% - an acceptable limit. Very short steps take a long time to compute and can also suffer rounding problems. This is an example of subdividing a rectangular mesh in the region in which curves will develop.
 Simple mesh of 382 elements equally spaced vertices. Superimposed are the polarity arrows (pointing bottom left) and purple and blue factors that control local growth rates.
|
 After one step in which the central region is subdivided twice (the vertices are too dense to see).
|
 As for the previous example but the size of the triangles forming the mesh is only visible along the edge. Red marks regions of increased growth on the top and bottom surfaces.
|
 After growing for 40 steps the excess growth has bent the canvas into an 'S' shape. Note the smooth bend is associated with about 4 elements.
|
 As for the previous example, but there has been no subdivision and the curves are unacceptably incomplete and jagged. There should always be enough elements to permit curves to develop smoothly.
|
Tutorial
Run the project (example project )for 40 steps to see the effect of:
- subdivision
Then change the modelname to 'NOSUBDIVISION', i.e. make a change to the interaction function
% Set up names for variant models. Useful for running multiple models on a cluster.
m.userdata.ranges.modelname.range = { 'NOSUBDIVISION', 'WITHSUBDIVISION' }; % CLUSTER
m.userdata.ranges.modelname.index = 2; % CLUSTER
by setting index to 1. With too few elements, bends are jagged.
The interaction function is shown below, red highlights the region of interest.
function m = gpt_demosubdivision_20121116( m ) %m = gpt_demosubdivision_20121116( m ) % Morphogen interaction function. % Written at 2012-11-16 12:30:06. % GFtbox revision 4351, . % The user may edit any part of this function between delimiters % of the form "USER CODE..." and "END OF USER CODE...". The % delimiters themselves must not be moved, edited, deleted, or added. if isempty(m), return; end fprintf( 1, '%s found in %s\n', mfilename(), which(mfilename()) ); try m = local_setproperties( m ); catch end realtime = m.globalDynamicProps.currenttime; %%% USER CODE: INITIALISATION % In this section you may modify the mesh in any way whatsoever. if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration % Zero out a lot of stuff to create a blank slate. % If no morphogens are set in the GUI it may be useful to % zero some arrays by uncommenting the following. % m.morphogens(:) = 0; % m.morphogenclamp(:) = 0; % m.mgen_production(:) = 0; % m.mgen_absorption(:) = 0; % m.seams(:) = false; % m.mgen_dilution(:) = false; % Set up names for variant models. Useful for running multiple models on a cluster. m.userdata.ranges.modelname.range = { 'NOSUBDIVISION', 'WITHSUBDIVISION' }; % CLUSTER m.userdata.ranges.modelname.index = 1; % CLUSTER end modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index}; % CLUSTER disp(sprintf('\nRunning %s model %s\n',mfilename, modelname)); % to plot polariser on the A side and resultant areal growth rate on the B side: m = leaf_plotoptions( m, 'morphogenA', 'KAPAR', 'morphogenB', 'KBPAR' ); %%% END OF USER CODE: INITIALISATION %%% SECTION 1: ACCESSING MORPHOGENS AND TIME. %%% AUTOMATICALLY GENERATED CODE: DO NOT EDIT. if isempty(m), return; end setGlobals(); global gNEW_KA_PAR gNEW_KA_PER gNEW_KB_PAR gNEW_KB_PER global gNEW_K_NOR gNEW_POLARISER gNEW_STRAINRET gNEW_ARREST dt = m.globalProps.timestep; polariser_i = gNEW_POLARISER; P = m.morphogens(:,polariser_i); [kapar_i,kapar_p,kapar_a,kapar_l] = getMgenLevels( m, 'KAPAR' ); [kaper_i,kaper_p,kaper_a,kaper_l] = getMgenLevels( m, 'KAPER' ); [kbpar_i,kbpar_p,kbpar_a,kbpar_l] = getMgenLevels( m, 'KBPAR' ); [kbper_i,kbper_p,kbper_a,kbper_l] = getMgenLevels( m, 'KBPER' ); [knor_i,knor_p,knor_a,knor_l] = getMgenLevels( m, 'KNOR' ); [strainret_i,strainret_p,strainret_a,strainret_l] = getMgenLevels( m, 'STRAINRET' ); [arrest_i,arrest_p,arrest_a,arrest_l] = getMgenLevels( m, 'ARREST' ); [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' ); [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' ); [id_subdivide_i,id_subdivide_p,id_subdivide_a,id_subdivide_l] = getMgenLevels( m, 'ID_SUBDIVIDE' ); % Mesh type: rectangle % base: 0 % centre: 0 % randomness: 0.1 % version: 1 % xdivs: 16 % xwidth: 16 % ydivs: 8 % ywidth: 8 % Morphogen Diffusion Decay Dilution Mutant % -------------------------------------------------- % KAPAR ---- ---- ---- ---- % KAPER ---- ---- ---- ---- % KBPAR ---- ---- ---- ---- % KBPER ---- ---- ---- ---- % KNOR ---- ---- ---- ---- % POLARISER ---- ---- ---- ---- % STRAINRET ---- ---- ---- ---- % ARREST ---- ---- ---- ---- % ID_A ---- ---- ---- ---- % ID_B ---- ---- ---- ---- % ID_SUBDIVIDE ---- ---- ---- ---- %%% USER CODE: MORPHOGEN INTERACTIONS % In this section you may modify the mesh in any way that does not % alter the set of nodes. if (Steps(m)==0) && m.globalDynamicProps.doinit % Initialisation code. % Put any code here that should only be performed at the start of % the simulation, for example, to set up initial morphogen values. id_a_p((m.nodes(:,1)>0)&(m.nodes(:,1)<=2))=1; id_b_p((m.nodes(:,1)>-2)&(m.nodes(:,1)<=0))=1; id_subdivide_p((m.nodes(:,1)>=-2)&(m.nodes(:,1)<=2))=1; P=m.nodes(:,1); end if realtime<=1 kapar_p(:) = 0; kaper_p(:) = 0; kbpar_p(:) = 0; kbper_p(:) = 0; knor_p(:) = 0; else kapar_p(:) = 0.05*id_a_p; kaper_p(:) = 0; kbpar_p(:) = 0.05*id_b_p; kbper_p(:) = 0; knor_p(:) = 0; end %%% END OF USER CODE: MORPHOGEN INTERACTIONS %%% SECTION 3: INSTALLING MODIFIED VALUES BACK INTO MESH STRUCTURE %%% AUTOMATICALLY GENERATED CODE: DO NOT EDIT. m.morphogens(:,polariser_i) = P; m.morphogens(:,kapar_i) = kapar_p; m.morphogens(:,kaper_i) = kaper_p; m.morphogens(:,kbpar_i) = kbpar_p; m.morphogens(:,kbper_i) = kbper_p; m.morphogens(:,knor_i) = knor_p; m.morphogens(:,strainret_i) = strainret_p; m.morphogens(:,arrest_i) = arrest_p; m.morphogens(:,id_a_i) = id_a_p; m.morphogens(:,id_b_i) = id_b_p; m.morphogens(:,id_subdivide_i) = id_subdivide_p; %%% USER CODE: FINALISATION % In this section you may modify the mesh in any way whatsoever. switch modelname case 'NOSUBDIVISION' % do nothing case 'WITHSUBDIVISION' % subdivide on the first step if realtime>0 && realtime<=0+dt m = leaf_subdivide( m, 'morphogen','id_subdivide',... 'min',0.5,'max',1,... 'mode','mid','levels','all'); end if realtime>0 && realtime<=0+dt m = leaf_subdivide( m, 'morphogen','id_subdivide',... 'min',0.5,'max',1,... 'mode','mid','levels','all'); end otherwise % If this happens, maybe you forgot a model. end %%% END OF USER CODE: FINALISATION end %%% USER CODE: SUBFUNCTIONS function m = local_setproperties( m ) % This function is called at time zero in the INITIALISATION section of the % interaction function. It provides commands to set each of the properties % that are contained in m.globalProps. Uncomment whichever ones you would % like to set yourself, and put in whatever value you want. % % Some of these properties are for internal use only and should never be % set by the user. At some point these will be moved into a different % component of m, but for the present, just don't change anything unless % you know what it is you're changing. % m = leaf_setproperty( m, 'trinodesvalid', true ); % m = leaf_setproperty( m, 'prismnodesvalid', true ); % m = leaf_setproperty( m, 'thicknessRelative', 1.600000 ); % m = leaf_setproperty( m, 'thicknessArea', 0.000000 ); % m = leaf_setproperty( m, 'thicknessMode', 'physical' ); % m = leaf_setproperty( m, 'activeGrowth', 1.000000 ); % m = leaf_setproperty( m, 'displayedGrowth', 1.000000 ); % m = leaf_setproperty( m, 'displayedMulti', [] ); % m = leaf_setproperty( m, 'allowNegativeGrowth', true ); % m = leaf_setproperty( m, 'usePrevDispAsEstimate', true ); % m = leaf_setproperty( m, 'perturbInitGrowthEstimate', 0.000010 ); % m = leaf_setproperty( m, 'perturbRelGrowthEstimate', 0.010000 ); % m = leaf_setproperty( m, 'perturbDiffusionEstimate', 0.000100 ); % m = leaf_setproperty( m, 'resetRand', false ); % m = leaf_setproperty( m, 'mingradient', 0.000000 ); % m = leaf_setproperty( m, 'relativepolgrad', false ); % m = leaf_setproperty( m, 'usefrozengradient', true ); % m = leaf_setproperty( m, 'userpolarisation', false ); % m = leaf_setproperty( m, 'thresholdsq', 4.000000 ); % m = leaf_setproperty( m, 'splitmargin', 1.400000 ); % m = leaf_setproperty( m, 'splitmorphogen', ); % m = leaf_setproperty( m, 'thresholdmgen', 0.500000 ); % m = leaf_setproperty( m, 'bulkmodulus', 1.000000 ); % m = leaf_setproperty( m, 'unitbulkmodulus', true ); % m = leaf_setproperty( m, 'poissonsRatio', 0.300000 ); % m = leaf_setproperty( m, 'starttime', 0.000000 ); % m = leaf_setproperty( m, 'timestep', 0.010000 ); % m = leaf_setproperty( m, 'timeunitname', ); % m = leaf_setproperty( m, 'distunitname', 'mm' ); % m = leaf_setproperty( m, 'scalebarvalue', 0.000000 ); % m = leaf_setproperty( m, 'validateMesh', true ); % m = leaf_setproperty( m, 'rectifyverticals', false ); % m = leaf_setproperty( m, 'allowSplitLongFEM', true ); % m = leaf_setproperty( m, 'longSplitThresholdPower', 0.000000 ); % m = leaf_setproperty( m, 'allowSplitBentFEM', false ); % m = leaf_setproperty( m, 'allowSplitBio', true ); % m = leaf_setproperty( m, 'allowFlipEdges', false ); % m = leaf_setproperty( m, 'allowElideEdges', true ); % m = leaf_setproperty( m, 'mincellangle', 0.200000 ); % m = leaf_setproperty( m, 'alwaysFlat', 0.000000 ); % m = leaf_setproperty( m, 'flattenforceconvex', true ); % m = leaf_setproperty( m, 'flatten', false ); % m = leaf_setproperty( m, 'flattenratio', 1.000000 ); % m = leaf_setproperty( m, 'useGrowthTensors', false ); % m = leaf_setproperty( m, 'plasticGrowth', false ); % m = leaf_setproperty( m, 'totalinternalrotation', 0.000000 ); % m = leaf_setproperty( m, 'stepinternalrotation', 2.000000 ); % m = leaf_setproperty( m, 'showinternalrotation', false ); % m = leaf_setproperty( m, 'performinternalrotation', false ); % m = leaf_setproperty( m, 'internallyrotated', false ); % m = leaf_setproperty( m, 'maxFEcells', 0 ); % m = leaf_setproperty( m, 'inittotalcells', 0 ); % m = leaf_setproperty( m, 'bioApresplitproc', ); % m = leaf_setproperty( m, 'bioApostsplitproc', ); % m = leaf_setproperty( m, 'maxBioAcells', 0 ); % m = leaf_setproperty( m, 'colors', (6 values) ); % m = leaf_setproperty( m, 'colorvariation', 0.050000 ); % m = leaf_setproperty( m, 'colorparams', (12 values) ); % m = leaf_setproperty( m, 'biocolormode', 'auto' ); % m = leaf_setproperty( m, 'freezing', 0.000000 ); % m = leaf_setproperty( m, 'canceldrift', false ); % m = leaf_setproperty( m, 'mgen_interaction', ); % m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_demosubdivision_20121116' ); % m = leaf_setproperty( m, 'allowInteraction', true ); % m = leaf_setproperty( m, 'interactionValid', true ); % m = leaf_setproperty( m, 'gaussInfo', (unknown type struct) ); % m = leaf_setproperty( m, 'stitchDFs', [] ); % m = leaf_setproperty( m, 'D', (36 values) ); % m = leaf_setproperty( m, 'C', (36 values) ); % m = leaf_setproperty( m, 'G', (6 values) ); % m = leaf_setproperty( m, 'solver', 'cgs' ); % m = leaf_setproperty( m, 'solverprecision', 'double' ); % m = leaf_setproperty( m, 'solvertolerance', 0.001000 ); % m = leaf_setproperty( m, 'solvertolerancemethod', 'max' ); % m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 ); % m = leaf_setproperty( m, 'allowsparse', true ); % m = leaf_setproperty( m, 'maxIters', 0 ); % m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 ); % m = leaf_setproperty( m, 'cgiters', 0 ); % m = leaf_setproperty( m, 'simsteps', 0 ); % m = leaf_setproperty( m, 'stepsperrender', 0 ); % m = leaf_setproperty( m, 'growthEnabled', true ); % m = leaf_setproperty( m, 'diffusionEnabled', true ); % m = leaf_setproperty( m, 'flashmovie', false ); % m = leaf_setproperty( m, 'makemovie', false ); % m = leaf_setproperty( m, 'moviefile', ); % m = leaf_setproperty( m, 'codec', 'None' ); % m = leaf_setproperty( m, 'autonamemovie', true ); % m = leaf_setproperty( m, 'overwritemovie', false ); % m = leaf_setproperty( m, 'framesize', [] ); % m = leaf_setproperty( m, 'mov', [] ); % m = leaf_setproperty( m, 'boingNeeded', false ); % m = leaf_setproperty( m, 'initialArea', 128.000000 ); % m = leaf_setproperty( m, 'bendunitlength', 11.313708 ); % m = leaf_setproperty( m, 'targetRelArea', 1.000000 ); % m = leaf_setproperty( m, 'defaultinterp', 'min' ); % m = leaf_setproperty( m, 'readonly', false ); % m = leaf_setproperty( m, 'projectdir', 'D:\ab\Matlab stuff\Growth models' ); % m = leaf_setproperty( m, 'modelname', 'GPT_DemoSubdivision_20121116' ); % m = leaf_setproperty( m, 'allowsave', true ); % m = leaf_setproperty( m, 'addedToPath', false ); % m = leaf_setproperty( m, 'bendsplit', 0.300000 ); % m = leaf_setproperty( m, 'usepolfreezebc', false ); % m = leaf_setproperty( m, 'dorsaltop', true ); % m = leaf_setproperty( m, 'defaultazimuth', -45.000000 ); % m = leaf_setproperty( m, 'defaultelevation', 33.750000 ); % m = leaf_setproperty( m, 'defaultroll', 0.000000 ); % m = leaf_setproperty( m, 'defaultViewParams', (unknown type struct) ); % m = leaf_setproperty( m, 'comment', ); % m = leaf_setproperty( m, 'legendTemplate', '%T: %q\n%m' ); % m = leaf_setproperty( m, 'bioAsplitcells', true ); % m = leaf_setproperty( m, 'bioApullin', 0.142857 ); % m = leaf_setproperty( m, 'bioAfakepull', 0.202073 ); % m = leaf_setproperty( m, 'interactive', false ); % m = leaf_setproperty( m, 'coderevision', 0 ); % m = leaf_setproperty( m, 'coderevisiondate', ); % m = leaf_setproperty( m, 'modelrevision', 0 ); % m = leaf_setproperty( m, 'modelrevisiondate', ); % m = leaf_setproperty( m, 'savedrunname', ); % m = leaf_setproperty( m, 'savedrundesc', ); % m = leaf_setproperty( m, 'vxgrad', (108 values) ); % m = leaf_setproperty( m, 'lengthscale', 16.000000 ); end % Here you may write any functions of your own, that you want to call from % the interaction function, but never need to call from outside it. % Remember that they do not have access to any variables except those % that you pass as parameters, and cannot change anything except by % returning new values as results. % Whichever section they are called from, they must respect the same % restrictions on what modifications they are allowed to make to the mesh. % For example: % function m = do_something( m ) % % Change m in some way. % end % Call it from the main body of the interaction function like this: % m = do_something( m );