Produced using MATLAB® software. TR_1D_model1_SS\DAE_SS_solver_1 Page 1 of 18 TR_1D_model1_SS\DAE_SS_solver_1.m % TR_1D_model1_ss\DAE_SS_solver_1.m % % function [x_state,iflag_converge,f,f_init] = ... % DAE_SS_solver_1(x_guess,Solver,... % func_calc_A_int,func_calc_b_int, ... % func_implement_BC,epsilon,Param); % % This procedure is a generic solver for a % DAE system arising from a discretized set % of PDE's of the form : % % epsilon(k) * df_dt(k) = b(k) - % \sum_{j} {A(k,j)*x_state(j)}. % % The procedure is passed an initial guess of % the state vector x_state, a structure of the % system parameters, and the names of the % subroutines that calculate from this information % the A matrix, the b vector, and the Jacobian % of the b vector. % % First, if Solver.max_iter_time is not 0, % the solver starts off with a stage of implicit % Euler time integration that robustly approaches % the vicinity of a stable steady state. If the % norm of the function (b-Ax) vector drops below % a magnitude Solver.atol_time, then the % time integration stage ends. % % If the time integration procedure has reduced % the function vector to an acceptable magnitude, % or if Solver.max_iter_time is 0 signifies % that Newton's method is to be performed directly, % then the Newton's method stage begins. If the % function vector norm drops below a magnitude of % Solver.atol_Newton, then convergence to steady % state is deemed to have occurred, and the routine % exists with iflag_converge = 1. Otherwise, the % Newton's method stage ends after a maximum of % Solver.max_iter_Newton number of iterations % with iflag_converge remaining zero or % negative (if an error occurred). % % If Solver.iflag_Adepend is non-zero, then % the value of the A matrix is state dependent and % must be recalculated at each iteration. If % Solver.iflag_nonneg is non-zero, then the 7/16/2002TR_1D_model1_SS\DAE_SS_solver_1 Page 2 of 18 % state variables are enforced to be non-negative % at every iteration. % % INPUT : % ======= % x_guess REAL(num_DOF) % This is the initial guess of the state % vector. % Solver This data structure contains the parameters % that control the operation of the solver. % The fields of this structure are : % .max_iter_time INT % the maximum number of implicit Euler time steps. % If =0, then no time simulation is performed and the % solver goes immediately to Newton's method % .dt REAL % the time step to be used in the implicit % Euler simulation % .atol_time REAL % the norm of the function (time derivative) vector % at which the time integration procedure is deemed % to have been sufficiently converged % .max_iter_Newton INT % the maximum number of Newton's method iterations % .atol_Newton REAL % the norm of the function (time derivative) vector % at which convergence to the steady state solution is % deemed to have been achieved % .iflag_Adepend INT % if this integer flag is non-zero, then the A matrix % is assumed to be state-dependent and so must be % recalculated at every iteration % .iflag_nonneg INT % if this integer flag is non-zero, then the elements % of the state vector are enforced to be non-negative % at every iteration % .iflag_verbose INT % if this integer flag is non-zero, then the solver % routine is instructed to print to the screen the % progress of the solution process; otherwise, it % runs silent % func_calc_A_int FUNCTION NAME % This is the name of the function that takes the % master state vector and the system parameters % and returns the A matrix that discretizes the % transport terms. % func_calc_b_int FUNCTION NAME % This is the name of the function that takes % the master state vector and the system % parameters and returns the b vector that % typically arises from the source term in a PDE. 7/16/2002TR_1D_model1_SS\DAE_SS_solver_1 Page 3 of 18 % It also returns the Jacobian of the b % vector, whose (m,n) element is the % partial of b(m) with respect to x_state(n). % func_implement_BC FUNCTION NAME % This is the name of the function that returns % the values of A, b, and bJac for the boundary % conditions, which are found where epsilon % has values of zero. % epsilon INT(num_DOF = length(x_guess)) % this is a 1-D array of integers of the same size % as x_state. It contains a 1 at position k for % every equation that is an ordinary differential % equation, and a 0 for every algebraic equation, % which in this problem arise from the boundary % conditions % Param % This data structure contains the system parameters % that are passed to the functions for calculating % the DAE system elements. % % OUTPUT : % ======== % x_state REAL(num_DOF) % This is the output estimate of the steady state % solution of the DAE system. % iflag_converge INT % This integer flag is set equal to 1 if the % solution procedure has converged. A value % of 0 means that the method did not converge. % A negative value indicates an error. % f REAL(num_DOF) % This is the time derivative vector (for boundary % points it is a measure of error in the boundary % condition) whose magnitude tells how far the % output estimate is from the steady state. % f_init REAL(num_DOF) % The time derivative vector for State_init % % Kenneth Beers % Massachusetts Institute of Technology % Department of Chemical Engineering % 7/2/2001 % % Version as of 7/23/2001 function [x_state,iflag_converge,f,f_init] = ... DAE_SS_solver_1(x_guess,Solver,... func_calc_A_int,func_calc_b_int, ... func_implement_BC,epsilon,Param); 7/16/2002TR_1D_model1_SS\DAE_SS_solver_1 Page 4 of 18 iflag_converge = 0; func_name = 'DAE_SS_solver'; % This integer flag controls what action to take % in case of an assertion failure. See the % assertion routines for further details. i_error = 2; % check input data % signify not yet completed checking of % initial data iflag_converge = -1; % check data structure that contains % solver parameters. We do this with a % routine that performs the appropriate % assertions on a structure. SolverType.num_fields = 8; % .max_iter_time ifield = 1; FieldType.name = 'max_iter_time'; FieldType.is_numeric = 1; FieldType.num_rows = 1; FieldType.num_columns = 1; FieldType.check_real = 1; FieldType.check_sign = 2; FieldType.check_int = 1; SolverType.field(ifield) = FieldType; % .dt ifield = 2; FieldType.name = 'dt'; FieldType.is_numeric = 1; FieldType.num_rows = 1; FieldType.num_columns = 1;