0001 function [rt_g m]=cielothermal_core(m,params,thermalparams,inputs,analysis)
0002
0003
0004
0005 if isempty(inputs),
0006 clear inputs;
0007 end
0008
0009
0010 if m.y_s.header.M==0,
0011 flag_y_s= 0;
0012 else
0013 flag_y_s= 1;
0014 end
0015
0016 converged=0;
0017
0018 params.kflag=1; params.pflag=1;
0019 params.ktflag=0;
0020
0021 dt=inputs.dt;
0022
0023 disp(sprintf('%s',' starting direct shooting method'));
0024
0025
0026 pressure=inputs.pressure;
0027
0028
0029 inputs.beta=thermalparams.beta;
0030 inputs.meltingpoint=thermalparams.meltingpoint;
0031 inputs.mintemp=thermalparams.mintemp;
0032
0033
0034 lst=m.lst;
0035
0036 count=1;
0037
0038 while(~converged),
0039
0040
0041 [rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,lst,m.ept,m.mpt,m.geom3,params,inputs,analysis);
0042
0043
0044 [rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s );
0045
0046
0047 rK_gg=IMdb('drop rK_gg');
0048 if ~isempty(rdK_gg), rdK_gg=IMdb('drop rdK_gg');end;
0049
0050
0051 [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s );
0052
0053
0054 [rp_g rK_fs]=IMdb('drop rp_g rK_fs');
0055
0056
0057 [rt_f]=Solver(rK_ff,rp_f,params.solverstring);
0058
0059
0060
0061 [rK_ff rp_f]=IMdb('drop rK_ff rp_f');
0062
0063
0064 rt_g= Mergesolvec( rt_f, m.G_mn, m.y_s );
0065
0066
0067 if (count>2),
0068 rt_f=IMdb('drop rt_f ');
0069 end
0070
0071
0072 [new_lst,converged,num_unstable_constraints]=CieloPenaltyConstraints(m.bgpdt,m.bgpdtb,m.est,lst, struct('pressure',pressure, 'dt',dt, 'temperature',rt_g,'beta',thermalparams.beta,'melting_point',thermalparams.meltingpoint),analysis);
0073
0074 if ~converged,
0075 disp(sprintf(' %s %i','#unstable constraints ',num_unstable_constraints));
0076 if num_unstable_constraints<thermalparams.min_thermal_constraints,
0077 converged=1;
0078 end
0079 if count==2,
0080 disp('Carefull! Screwed up solutoin on purpose!');
0081 converged=1;
0082 end
0083 end
0084
0085
0086 lst=Equiv(lst,new_lst);
0087
0088 count=count+1;
0089 end
0090
0091
0092 m.lst=lst;
0093
0094 end