Subroutine Type

Modeling

Definition

The TUNSUB allows you to modify certain solver parameters during a simulation.

Use

<Param_Simulation

    usrsub_param_string = "USER(9998, 0.1)"

    usrsub_dll_name     = "NULL"

    usrsub_fnc_name     = "TUNSUB"

 />

Calling Syntax

Fortran

SUBROUTINE TUNSUB (TIME, PAR, NPAR, IFLAG, ISTATE, H, IORDER)

C

void TUNSUB (double *time, double *par, int *npar, int *iflag, int *istate, double *h, int *iorder)

Python

def TUNSUB(time, par, npar, iflag, istate, h, iorder)

MATLAB

function TUNSUB(time, par, npar, iflag, istate, h, iorder)

Input Arguments

[double precision] TIME

Current simulation time.

[double precision] PAR

An array that contains the constant arguments from the list provided in the user-defined statement.

[integer] NPAR

Number of entries in the PAR array.

[integer] IFLAG

An integer variable that MotionSolve sets to 1 when the user subroutine is called during MotionSolve initialization.  The IFLAG will be set to zero at all other times.

[integer] ISTATE

An integer variable that MotionSolve sets to 1 when the integrator step is a successful step.  If the integrator step is not successful, ISTATE will be set to zero.

[double precision] H

Current integrator step size (h < h_max, specified in Param_Transient).

[integer] IORDER

Current integrator order (IORDER < max_order, specified in Param_Transient).

Output Values

None

Example

The following TUNSUB example modifies the absolute and relative integrator error tolerance to the value passed in as the second argument in the USER(…) in the model when the angular velocity of the rigid bodies exceed 20 radians/second.

from math import fabs

g_num_rigid_body = 0

g_cg_id_vec = []

g_body_id_vec = []

#define the TUNSUB function

def TUNSUB(time, par, npar, iflag, istate, h, iorder):

    w = 3*[0.0]

    global g_num_rigid_body

    global g_cg_id_vec

    global g_body_id_vec

    #set solver parameters

    py_set_dae_error(0.01)

    py_set_dae_hmax(0.01)

    py_set_dae_maxord(5)

    py_set_dae_maxit(4)

    t0 = py_get_starting_time()

    te = py_get_end_time()

    #define dependencies during iflag = 1

    if iflag==1:        

        tmp_num = 0

        tmp_num = py_getnumid("PART")

        tmp_iarray1 = tmp_num*[0]

        tmp_iarray2 = tmp_num*[0]

        [tmp_iarray1, istat] = py_getidlist("PART", tmp_num)

        g_num_rigid_body = 0 

        for i in xrange(tmp_num):

            [tmp_str, istat] = py_modfnc("Body_Rigid", tmp_iarray1[i], "cg_id")

            cg_id = int(tmp_str)

            [tmp_str, istat] = py_modfnc("Body_Rigid", tmp_iarray1[i], "IsGround")

            if tmp_str=='TRUE':

                continue

            tmp_iarray2[g_num_rigid_body] = tmp_iarray1[i]

            tmp_iarray1[g_num_rigid_body] = cg_id

            g_num_rigid_body +=1 

        g_cg_id_vec = tmp_iarray1

        g_body_id_vec = tmp_iarray2

        for i in xrange(g_num_rigid_body):

            cg_id = g_cg_id_vec[i];

            [w[0], istat] = py_sysfnc("WX", cg_id)

            [w[1], istat] = py_sysfnc("WY", cg_id)

            [w[2], istat] = py_sysfnc("WZ", cg_id)

        return

#check if current step is converged i.e. istate = 1

    if istate==1:

        for i in xrange(g_num_rigid_body):

            cg_id = g_cg_id_vec[i];

            #obtain angular velocities for all parts

            [w[0], istat] = py_sysfnc("WX", cg_id)

            [w[1], istat] = py_sysfnc("WY", cg_id)

            [w[2], istat] = py_sysfnc("WZ", cg_id)

#check if the angular velocities are above some threshold

#value

            if (fabs(w[0])>20) | (fabs(w[1])>20) | (fabs(w[2])>20):

#change tolerance values

                istat = py_set_tol_factor("Body_Rigid", g_body_id_vec[i], "ATOL", par[1])

                istat = py_set_tol_factor("Body_Rigid", g_body_id_vec[i], "RTOL", par[1])

Comments

set_tol_factor(char *entity_type, int entity_id, char *tol_type, double factor, int *errflg)

set_dae_error(double error);

set_dae_hmax(double h_max);

set_dae_maxord(int maxord);

set_dae_maxit(int maxit);

set_dae_jacobeval(int jacob_eval);

set_dae_jacobinit(int jacob_init);