Index: sm/trunk/src/m/solutions/cielo/Builduset.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Builduset.m (revision 30) +++ (revision ) @@ -1,82 +1,0 @@ -function uset=builduset( bgpdt, bgpdtb,analysis); -% BUILDUSET builduset(bgpdt, bgpdtb) -% -% Builds the uset table from the bgpdt table. -% -% Input: bgpdt, bgpdtb pointer to bgpdt table on data base -% analysis: analysis type, string -% -% Output: uset structure of partitioning vectors and its sizes -% in global Matlab workspace -% -% Calls Upartn (mex) -% where ms= m-script, mf= m-function, mex= executable function (c-code) -% -% This function is called after the uset bits in the bgpdt table have been -% determined, e.g. after Gp4 or after Gpsp. -% -% The uset table is a matlab structure containing -% -% uset.pv_g g-set partitioning vector -% uset.gsize size of the g-set -% uset.pv_n n-set partitioning vector with respect to (wrt) the g-set -% uset.nsize size of the n-set -% uset.pv_m m-set partitioning vector wrt the g-set -% uset.msize size of the m-set -% uset.pv_s s-set partitioning vector wrt the n-set -% uset.ssize size of the s-set -% uset.pv_f f-set partitioning vector wrt the n-set -% etc. -% -% The sets have the following hierarchy -% -% p - e p = e + g p= physical dof, e= extra point dofs -% - g - m g= n + m g= grid dofs, m= multi-point constraint dofs -% - n - s n= f + s n= dofs not constr.by m, s= single-point constr.dofs -% - f - o f= a + o f= free dofs, o= omitted dofs -% - a - q a= t + q a= assembled dofs, q= generalized modal dofs -% - t - r t= l + r t= total dofs, r= rigid body dofs -% - l - c l= c + b + lm l= left-over dofs, c= free dofs, -% - b b= fixed dofs for modal reduction, -% - lm lm= Lagrange multipliers. -% -% any combination of sets can be constructed, e.g. -% d= a + e d= dynamic dofs - -% -% Go down the set hierarchy and fill the structure -% note that the order m - s - o - q - r must be kept, the m- and n-set partitioning -% vector are wrt the g-set, the s- and f-set partitioning vectors are wrt the n-set, -% etc. -% - -% split g into m and n - -[uset.pv_m, uset.pv_n]= Upartn( bgpdt, bgpdtb,analysis, 'g', 'm', 'n'); -uset.msize= length(uset.pv_m); -uset.nsize= length(uset.pv_n); -uset.gsize= uset.msize + uset.nsize; - -% split n into s and f - -[uset.pv_s, uset.pv_f]= Upartn( bgpdt, bgpdtb,analysis, 'n', 's', 'f'); -uset.ssize= length(uset.pv_s); -uset.fsize= length(uset.pv_f); - -% split f into o and a - -[uset.pv_o, uset.pv_a]= Upartn( bgpdt, bgpdtb,analysis, 'f', 'o', 'a'); -uset.osize= length(uset.pv_o); -uset.asize= length(uset.pv_a); - -% split a into q and t - -[uset.pv_q, uset.pv_t]= Upartn( bgpdt, bgpdtb,analysis, 'a', 'q', 't'); -uset.qsize= length(uset.pv_q); -uset.tsize= length(uset.pv_t); - -% split t into r and l - -[uset.pv_r, uset.pv_l]= Upartn( bgpdt, bgpdtb,analysis, 't', 'r', 'l'); -uset.rsize= length(uset.pv_r); -uset.lsize= length(uset.pv_l); Index: sm/trunk/src/m/solutions/cielo/CieloHorizontalVelocityDepthAverage.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/CieloHorizontalVelocityDepthAverage.m (revision 30) +++ (revision ) @@ -1,12 +1,0 @@ -function velocity_average_g=CieloHorizontalVelocityDepthAverage(md,velocity_g); -%compute depth averaged horizontal velocity - -vx_average2d=DepthAverage(md,velocity_g(1:6:uset.gsize)); -vy_average2d=DepthAverage(md,velocity_g(2:6:uset.gsize)); - -%project onto layer 1 of 3d mesh -vx_average=project3d(md,vx_average2d,'node',1); -vy_average=project3d(md,vy_average2d,'node',1); - -%output on g-set -velocity_average_g=zeros(uset.gsize,1);velocity_average_g(1:6:uset.gsize)=vx_average;velocity_average_g(2:6:uset.gsize)=vy_average; Index: sm/trunk/src/m/solutions/cielo/CieloVelocityExtrude.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/CieloVelocityExtrude.m (revision 30) +++ (revision ) @@ -1,39 +1,0 @@ -function u_g=CieloVelocityExtrude(md,u_g); -%Extrude collapsed horizontal velocities present in MacAyeal elements - -if strcmpi(md.type,'2d'), - %do nothing, we don't have a 3d mesh anyway. -else - - %Find list of 2d grids that belong to the collapsed macayeal elements. - %Numbering for these 2d grids corresponds to the initial 2d mesh. - grids2d=md.proj2d_coll_grids(1:md.number_grids2d_coll); - vx2d=u_g(1:6:6*md.number_grids2d_coll); - vy2d=u_g(2:6:6*md.number_grids2d_coll); - - %Plug values for these 2d grids into the 2d mesh: - vx=zeros(md.numberofgrids2d,1); - vx(grids2d)=vx2d; - - vy=zeros(md.numberofgrids2d,1); - vy(grids2d)=vy2d; - - %Extrude across the 3d mesh - vx3d=project3d(md,vx,'node'); - vx3d(1:md.number_grids2d_coll)=0; - vy3d=project3d(md,vy,'node'); - vy3d(1:md.number_grids2d_coll)=0; - - %For border grids between MacAyeal and Pattyn, values will double, fix that. - for i=1:length(md.penalties), - vx3d(md.penalties(i).grids3d_noncollapsed)=0; - vy3d(md.penalties(i).grids3d_noncollapsed)=0; - end - - %Plug back into u_g - u_g(1:6:end)=u_g(1:6:end)+vx3d; - u_g(2:6:end)=u_g(2:6:end)+vy3d; - -end - -end %end function Index: sm/trunk/src/m/solutions/cielo/ConstrainControlParameterB.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/ConstrainControlParameterB.m (revision 30) +++ (revision ) @@ -1,7 +1,0 @@ -function B=ConstrainControlParameterB(B,m,params); - -%Ok, simple one: keep B positive. - -pos=find(B<=0); -B(pos)=10; - Index: sm/trunk/src/m/solutions/cielo/ConstrainControlParameterdrag.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/ConstrainControlParameterdrag.m (revision 30) +++ (revision ) @@ -1,5 +1,0 @@ -function drag=ConstrainControlParameterdrag(drag,m,params); -%Ok, simple one: keep drag positive. - -pos=find(drag<=0); -drag(pos)=5; Index: /issm/trunk/src/m/solutions/cielo/ControlOptions.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/ControlOptions.m (revision 30) +++ /issm/trunk/src/m/solutions/cielo/ControlOptions.m (revision 31) @@ -12,4 +12,2 @@ options=optimset(options,'TolX',parameters.tolx); options=optimset(options,'FunValCheck','on'); - - Index: sm/trunk/src/m/solutions/cielo/Mergesolvec.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Mergesolvec.m (revision 30) +++ (revision ) @@ -1,57 +1,0 @@ -function [u_g]= mergesolvec( u_f, G_mn, y_s ,uset); -% MERGESOLVEC [u_g]= mergesolvec ( u_f, G_mn, y_s ) -% Merges solution vectors (displacement/temperature) from f size to g size, -% u_f to u_g, -% the solution vectors can be total in linear or incremental in nonlinear. -% -% g global dof g = m + n -% m dependent dof -% n independent dof n= s + f -% f free dof -% s single point constraint dof -% -% Input: u_f solution vectors (f size) -% in Matlab workspace -% G_mn Reduction matrix from constraints (m x n) -% in Matlab workspace -% y_s vectors of enforced displacements y_s -% in Matlab workspace -% -% Output: u_g solution vectors (g size) -% -% Input from global workspace: -% -% uset uset.pv_* partitioning vectors and sizes -% uset.*size *=n,m etc. -% -% Called by sol101 (ms), Normalmodes (ms), sol159 (ms) (all solution sequences) -% Calls none -% where ms= m-script, mf= m-function, mex= executable function (c-code) -% -% m_set and n_set index vectors are with respect to (wrt) the g- set -% s_set and f_set index vectors are with respect to (wrt) the n- set - -%Merge f_set and s_set to n_set - -ncol= size( u_f, 2); -u_n= sparse(uset.nsize,ncol); -u_n(uset.pv_f,:)= u_f; - - -% Note: we assume u_f and y_s have the same no.of columns -if ~isempty(y_s) - u_n(uset.pv_s,:)= y_s; -end - -% Merge n_set and m_set to g_set - -u_g = sparse(uset.gsize,ncol); -u_g(uset.pv_n,:)= u_n; - - -if (uset.msize>0) - u_m = G_mn * u_n; - u_g(uset.pv_m,:)= u_m; -end - -end Index: sm/trunk/src/m/solutions/cielo/NewCreateFEMModel.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/NewCreateFEMModel.m (revision 30) +++ (revision ) @@ -1,105 +1,0 @@ -%----------------------------------------------------------------------- -% CreateFEMModel.m: from a matlab ISSM @model class object, create a finite element -% compatible with Cielo. -% -%----------------------------------------------------------------------- -function m=CreateFEMModel(md,analysis_type ) - - global ParameterDefaults - -% read input data: - - disp(sprintf('\n reading data from model %s...',md.name)); - [m.ncase, m.dynamic, m.ept, m.geom1, m.geom1b, m.geom2, m.part, m.tpart, m.geom3, ... - m.geom4, m.mpt, m.pvt, m.sets, m.dit] = newImp( md ,analysis_type); - - error; - -% generate list of states from case control (ncase), if any: - - disp(' generating states...'); - [m.states] = GenerateStates( m.ncase, m.dynamic); - - %ncase not needed anymore - m=rmfield(m,'ncase'); - -% generate grid point definition table -% (internal grid records, contain sil and uset data) - - disp(' generating gpdt...'); - [m.gpdt,m.gpdtb] = Gp1( m.geom1, m.geom1b); - - -% generate coordinate system transformation matrices, transform -% grids to basic coordinate system - - disp(' generating updated gpdt...'); - [m.cstm, m.bgpdt, m.bgpdtb] = Gp1b( m.geom1, m.geom1b, m.gpdt, m.gpdtb); - - m=rmfield(m,'gpdt'); - if ~isempty(m.gpdtb), m=rmfield(md,'gpdtb');end; - - %geom1 and geom1b not needed anymore - m=rmfield(m,'geom1'); - if ~isempty(m.geom1b), m=rmfield(m,'geom1b'); end; - - %disp(' generating element objects...'); - [m.est,m.lst] = Gp2( m.bgpdt, m.bgpdtb, m.ept, m.geom2, m.mpt, m.geom3); - - %geom2 not needed anymore - m=rmfield(m,'geom2'); - -% phase 1: -% - -% set parameters using bulk data, nastran case control above subcase level: - current_state = {}; - [next_state] = NextState( m.states, current_state); - m.current_state = next_state; - clear next_state current_state; - - [m.params] = Pvt( m.current_state, m.pvt); - - -% generate usets and constraint matrices: - - disp(' generating uset data and constraints...'); - - [m.updated_bgpdt, m.updated_bgpdtb, m.y_g, m.R_mg] = Gp4( m.states, m.geom4, m.bgpdt, m.bgpdtb,analysis_type); - - - %geom4 not needed anymore - if ~isempty(m.geom4), m=rmfield(m,'geom4');end; - -% replace old bgpdt with new (bgpdt unchanged if nonexistent updated_bgpdt): - m.bgpdt = m.updated_bgpdt;m=rmfield(m,'updated_bgpdt'); - if ~isempty(m.updated_bgpdtb), - m.bgpdtb = m.updated_bgpdtb; - m=rmfield(m,'updated_bgpdtb'); - end; - -% create uset table - m.uset=Builduset(m.bgpdt,m.bgpdtb,analysis_type); - - -% -% phase 1: -% -% autospc operations: -% autospc is done on the g_set, it should be done on the n_set later - -% Partition R_mg into R_mm and R_mn, and -% calculate G_mn = - R_mm^-1 * R_mn (MCE1 type operations) - - [m.G_mn] = Normalizeconstraints( m.R_mg, m.uset); - -% Reduce enforced displacements / temperatures - [m.y_s] = Reducevectorg( m.y_g ,m.uset); - - - if isempty(m.y_s), - flag_y_s= 0; - else - flag_y_s= 1; - end -end Index: sm/trunk/src/m/solutions/cielo/Orth.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Orth.m (revision 30) +++ (revision ) @@ -1,30 +1,0 @@ -function new_gradients=Orth(gradients,oldgradients); - -global client_server_mode IMDataCounter - -%outpu -new_gradients=gradients; - -numparams=length(gradients); -for i=1:numparams, - - %recover gradient #i - gradient=gradients{i}; - oldgradient=oldgradients{i}; - - if isempty(oldgradient), - oldgradient=zeros(length(gradient),1); - end - - % normalize gradient #i to 1 - gradient=gradient/norm(gradient); - - % rough orthagonalization - gradient=gradient-(gradient'*oldgradient)*oldgradient; - - % renormalize gradients - gradient=gradient/norm(gradient); - - %plug back into gradients - new_gradients{i}=gradient; -end Index: sm/trunk/src/m/solutions/cielo/Reducematrix.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Reducematrix.m (revision 30) +++ (revision ) @@ -1,61 +1,0 @@ -function [K_nn]= reducematrix( K_gg, G_mn, pv_m_set, pv_n_set, msize, flag ); -% REDUCEMATRIX [K_nn]= reducematrix( K_gg, G_mn, pv_m_set, pv_n_set, msize, flag ) -% Reduces K_gg to K_nn -% -% K_nn^~ K_nm -% K_gg = [ ] -% K_mn K_mm -% -% if flag == 1 -% K_nn = K_nn^~ + K_nm * G_mn -% else -% K_nn = K_nn^~ + K_nm * G_mn + G_mn^T * K_mn + G_mn^T * K_mm * G_mn -% -% m dof to be condensed -% n dof to be retained -% g total dof g = m + n -% -% Input: K_gg matrix (g x g) -% G_mn Reduction matrix from constraints (m x n) -% pv_m_set partitioning vector for the m dof wrt g -% pv_n_set partitioning vector for the n dof wrt g -% msize number of m-set dof -% flag flag for type of condensation, see above -% -% Output: K_nn reduced matrix (n x n) -% -% Called by Reducematrixfromgtof (mf), Normalmodes (ms) -% Calls none -% where ms= m-script, mf= m-function, mex= executable function (c-code) -% -% pv_m_set and pv_n_set index vectors are with respect to the g- set. -% -% All input and output variables are assumed to be in the workspace. - -if msize > 0 - -% Partition K_gg - - K_mm = K_gg( pv_m_set, pv_m_set); - K_mn = K_gg( pv_m_set, pv_n_set); - K_nm = K_gg( pv_n_set, pv_m_set); - K_nn = K_gg( pv_n_set, pv_n_set); - -% Reduce K_gg to K_nn - - K_nn = K_nn + K_nm * G_mn; - - if flag ~= 1 - - K_mn = K_mn + K_mm * G_mn; - K_nn = K_nn + G_mn' * K_mn; - - end - -else - -% msize=0 just return input - - K_nn = K_gg; - -end Index: sm/trunk/src/m/solutions/cielo/Reducematrixfromgtof.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Reducematrixfromgtof.m (revision 30) +++ (revision ) @@ -1,51 +1,0 @@ -function [ K_ff, K_fs] = reducematrixfromgtof(K_gg, G_mn, flag_u_s,uset); -% REDUCEMATRIXFROMGTOF [K_ff, K_fs] = Reducematrixfromgtof( K_gg, G_mn, flag_u_s ) -% Reduces a matrix from g-size to f-size -% -% Input: -% K_gg Matrix in g-size -% G_mn Reduction matrix -% flag_u_s flag indicating whether K_fs is calculated or not -% =0 K_fs is not calculated, an empty matrix is returned -% =1 K_fs is calculated and returned -% Output: -% K_ff Reduced system tangent matrix f-size x f-size -% K_fs Reduced sytem tangent matrix f-size x s-size -% -% Input from global workspace: -% -% uset uset.pv_* partitioning vector -% uset.*size size -% *=n,m etc., see Builduset -% -% Called by sol101 (ms), Normalmodes (ms), Tangentmatrixheat (mf) -% Calls Reducematrix (mf) -% where ms= m-script, mf= m-function, mex= executable function (c-code) -% -% The g-set is partitioned as follows -% g = n + m -% n = f + s -% g all system degrees of freedom (dof) -% m multi-point constraint dofs to be condensed out -% n remaining dof after m-set is condensed from the g-set -% s boundary condition dofs (single point constraints) -% f remaining free dofs after s-set is eliminated from the n-set -% - -% Reduce matrix from g-size to n-size - -K_nn= Reducematrix( K_gg, G_mn, uset.pv_m, uset.pv_n, uset.msize, 2 ); - -% Reduce matrix from n-size to f-size - -if uset.ssize>0 - K_ff= K_nn(uset.pv_f, uset.pv_f); - if flag_u_s - K_fs= K_nn(uset.pv_f,uset.pv_s); - else - K_fs=[]; - end -else - K_ff= K_nn; - K_fs=[]; -end Index: sm/trunk/src/m/solutions/cielo/Reducerightside.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Reducerightside.m (revision 30) +++ (revision ) @@ -1,111 +1,0 @@ -function [p_f]= reducerightside( p_g, G_mn, K_fs, y_s, flag_y_s, uset ) -% REDUCERIGHTSIDE [p_f]= Reducerightside( p_g, G_mn, K_fs, y_s, flag_y_s ) -% Reduces right hand side vectors from g-size to f-size (p_g to p_f). -% -% Partition p_g into p_n and p_m, then modify p_n -% -% p_n = p_n + G_mn^T * p_m -% -% where G_mn is the reduction matrix from multi-point constraints -% -% Partition p_n into p_f and p_s, then modify p_f if flag_y_s > 0 -% -% p_f = p_f - K_fs * y_s -% -% where y_s are enforced displacements/temperatures, -% -% The g-set is partitioned as follows -% g = n + m -% n = f + s -% g all system degrees of freedom (dof) -% m multi-point constraint dofs to be condensed out -% n remaining dof after m-set is condensed from the g-set -% s boundary condition dofs (single point constraints) -% f remaining free dofs after s-set is eliminated from the n-set -% -% Input: p_g right hand side vectors (g size) -% G_mn Reduction matrix from constraints (m x n) -% K_fs Partitioned stiffness matrix (f x s) -% y_s enforced displacements / temperatures s size -% flag_y_s flag for enforced displacement / temperature -% > 0 K_fs * y_s is subtracted -% <= 0 no action -% -% Output: p_f reduced right hand side vectors (f-size) -% -% Called by sol101 (ms), Normalmodes (ms), Residualheat (mf) -% Calls none -% where ms= m-script, mf= m-function, mex= executable function (c-code) -% -% All input and output variables are assumed to be in the Matlab workspace. -% The load vectors p_g and y_s must have the same number of columns. -% -% Input from global workspace: -% -% uset uset.pv_* partitioning vector -% uset.*size size -% *=n,m etc., see Builduset -% -% uset.pv_m and uset.pv_n index vectors are wrt the g- set -% uset.pv_s and uset.pv_f index vectors are wrt the n- set -% etc. - -% Reduce p_g to p_f if not empty - -if ~isempty(p_g) - - % Reduce p_g to p_n - - if ( uset.msize > 0 ) - p_n = p_g( uset.pv_n, :); - p_m = p_g( uset.pv_m, :); - p_n = p_n + G_mn' * p_m; - else - p_n = p_g; - end - - % Reduce p_n to p_f - - if ( uset.ssize > 0 ) - p_f = p_n( uset.pv_f, :); - else - p_f = p_n; - end - -% Create a p_f if p_g is empty - -else - - if isempty(y_s) - disp('No right hand side found...'); - p_f=[]; - return; - else - ncols= size( y_s, 2); - p_f= sparse([],[],[],uset.fsize, ncols); - end - -end - -% for nonzero boundary conditions, subtract coupling forces, -% this operation is only executed if flag_y_s > 0, -% the flag is set in the calling script for linear analysis -% and the first iteration step in nonlinear analysis - -if ( flag_y_s > 0 ) & ( ~isempty(y_s) ) - % check if p_f and y_s have the same number of columns - - ncolp= size(p_f,2); - ncoly= size(y_s,2); - - if ncolp ~= ncoly - disp('Right hand side can not be calculated...'); - return; - end - - p_f = p_f - K_fs * y_s; - -end - - -end Index: sm/trunk/src/m/solutions/cielo/Reducevectorg.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Reducevectorg.m (revision 30) +++ (revision ) @@ -1,19 +1,0 @@ -function [y_s]= Reducevectorg( y_g,uset ) -% REDUCEVECTORG [y_s ]= reducevectorg( py_g, uset ) -% -% -% Input in global workspace: -% -% uset uset.pv_* partitioning vectors and sizes -% uset.*size *=n,m etc. - -if ( uset.ssize>0 ) - if ~isempty(y_g), - y_n= y_g(uset.pv_n,:); - y_s= y_n(uset.pv_s,:); - else - y_s=[]; - end -else - y_s=[]; -end Index: sm/trunk/src/m/solutions/cielo/Residualice.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/Residualice.m (revision 30) +++ (revision ) @@ -1,3 +1,0 @@ -function [res_g] = Residualice( KH, pnl, u); - res_g = pnl - KH * u; -end Index: sm/trunk/src/m/solutions/cielo/SetParameterDefaults.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/SetParameterDefaults.m (revision 30) +++ (revision ) @@ -1,46 +1,0 @@ -%SETPARAMETERDEFAULTS -% SetParameterDefaults sets parameter defaults for solution sequences. -% In particular, Emg module defaults, as well a Client/Server parameters. -% -% See also Pvt - - -global ParameterDefaults; - -ParameterDefaults.autospc ='yes'; -ParameterDefaults.coupmass = -1; % translated to lumped/consistent option -ParameterDefaults.epzero = 1.e-8; -ParameterDefaults.sigma = 0.; -ParameterDefaults.tabs = 0.; - -% for testing, may be made permanent later: - -ParameterDefaults.linesrch = 1; % 0=no, 1=yes - -% Client/Server mode settings - -ParameterDefaults.tunnel='no'; %enable ssh tunnelling for encrypted communication to Server. -ParameterDefaults.batsiz=1; %size of batch jobs sent to nodes -ParameterDefaults.fdbckout=0; %feedback output level from server (0: no output, 1: basic output, 3: large output) -ParameterDefaults.batch='no'; %run in batch mode. -ParameterDefaults.errchk='yes'; %run with error chechking. - -% Integration rule for sol159 is set to 3 (mid point rule) - - ParameterDefaults.int_rule = 3; - - -% John's Emg.mex, and Plotresults.m are picking up the following parameters - - ParameterDefaults.analysis='structural'; - ParameterDefaults.zgrad ='no'; - ParameterDefaults.imoselem='no'; - ParameterDefaults.wtmass = 1.; - ParameterDefaults.lumpm = 1; - ParameterDefaults.mred = 0; - ParameterDefaults.niref = 1; - ParameterDefaults.njref = 1; - ParameterDefaults.nkref = 1; - ParameterDefaults.cido =-1; - ParameterDefaults.solverstring = ' -mat_type aij -ksp_type cgs -pc_type asm -sub_pc_type lu '; - Index: sm/trunk/src/m/solutions/cielo/SwitchPartitioning.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/SwitchPartitioning.m (revision 30) +++ (revision ) @@ -1,30 +1,0 @@ -function new_u_g=SwitchPartitioning(u_g,destination,part,tpart,doflist); -%INPUT function new_u_g=SwitchPartitioning(u_g,destination,part,tpart,doflist); - -gsize=length(u_g); -if strcmpi(destination,'workspace'), - %g-size vector u_g is ordered according to a parallel partition. REorder it - %so that u_g correspondgs to the workspace partition (ie: the initial grid - %partition) - new_u_g=zeros(gsize,1); - - %go through list of degrees of freedom that need repartitioning - for i=1:length(doflist), - dof=doflist(i); - ind=dof:6:gsize; ind=ind(tpart); - new_u_g(dof:6:gsize)=u_g(ind); - end -elseif strcmpi(destination,'cluster'), - %g size vector u_g is ordered according to the matlab workspace. We need to - %reorder it so that u_g is now ordered in the parallel partition. - - new_u_g=zeros(gsize,1); - - for i=1:length(doflist), - dof=doflist(i); - ind=dof:6:gsize; vector=u_g(ind); - new_u_g(dof:6:gsize)=vector(part); - end -else - error('SwitchPartitioning error message: only types of destination allowed are ''workspace'' or ''cluster''!'); -end Index: /issm/trunk/src/m/solutions/cielo/cielo/CieloHorizontalVelocityDepthAverage.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/CieloHorizontalVelocityDepthAverage.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/CieloHorizontalVelocityDepthAverage.m (revision 31) @@ -0,0 +1,12 @@ +function velocity_average_g=CieloHorizontalVelocityDepthAverage(md,velocity_g); +%compute depth averaged horizontal velocity + +vx_average2d=DepthAverage(md,velocity_g(1:6:uset.gsize)); +vy_average2d=DepthAverage(md,velocity_g(2:6:uset.gsize)); + +%project onto layer 1 of 3d mesh +vx_average=project3d(md,vx_average2d,'node',1); +vy_average=project3d(md,vy_average2d,'node',1); + +%output on g-set +velocity_average_g=zeros(uset.gsize,1);velocity_average_g(1:6:uset.gsize)=vx_average;velocity_average_g(2:6:uset.gsize)=vy_average; Index: /issm/trunk/src/m/solutions/cielo/cielo/CieloVelocityExtrude.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/CieloVelocityExtrude.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/CieloVelocityExtrude.m (revision 31) @@ -0,0 +1,39 @@ +function u_g=CieloVelocityExtrude(md,u_g); +%Extrude collapsed horizontal velocities present in MacAyeal elements + +if strcmpi(md.type,'2d'), + %do nothing, we don't have a 3d mesh anyway. +else + + %Find list of 2d grids that belong to the collapsed macayeal elements. + %Numbering for these 2d grids corresponds to the initial 2d mesh. + grids2d=md.proj2d_coll_grids(1:md.number_grids2d_coll); + vx2d=u_g(1:6:6*md.number_grids2d_coll); + vy2d=u_g(2:6:6*md.number_grids2d_coll); + + %Plug values for these 2d grids into the 2d mesh: + vx=zeros(md.numberofgrids2d,1); + vx(grids2d)=vx2d; + + vy=zeros(md.numberofgrids2d,1); + vy(grids2d)=vy2d; + + %Extrude across the 3d mesh + vx3d=project3d(md,vx,'node'); + vx3d(1:md.number_grids2d_coll)=0; + vy3d=project3d(md,vy,'node'); + vy3d(1:md.number_grids2d_coll)=0; + + %For border grids between MacAyeal and Pattyn, values will double, fix that. + for i=1:length(md.penalties), + vx3d(md.penalties(i).grids3d_noncollapsed)=0; + vy3d(md.penalties(i).grids3d_noncollapsed)=0; + end + + %Plug back into u_g + u_g(1:6:end)=u_g(1:6:end)+vx3d; + u_g(2:6:end)=u_g(2:6:end)+vy3d; + +end + +end %end function Index: /issm/trunk/src/m/solutions/cielo/cielo/ConstrainControlParameterB.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/ConstrainControlParameterB.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/ConstrainControlParameterB.m (revision 31) @@ -0,0 +1,7 @@ +function B=ConstrainControlParameterB(B,m,params); + +%Ok, simple one: keep B positive. + +pos=find(B<=0); +B(pos)=10; + Index: /issm/trunk/src/m/solutions/cielo/cielo/ConstrainControlParameterdrag.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/ConstrainControlParameterdrag.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/ConstrainControlParameterdrag.m (revision 31) @@ -0,0 +1,5 @@ +function drag=ConstrainControlParameterdrag(drag,m,params); +%Ok, simple one: keep drag positive. + +pos=find(drag<=0); +drag(pos)=5; Index: /issm/trunk/src/m/solutions/cielo/cielo/Mergesolvec.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/Mergesolvec.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/Mergesolvec.m (revision 31) @@ -0,0 +1,57 @@ +function [u_g]= mergesolvec( u_f, G_mn, y_s ,uset); +% MERGESOLVEC [u_g]= mergesolvec ( u_f, G_mn, y_s ) +% Merges solution vectors (displacement/temperature) from f size to g size, +% u_f to u_g, +% the solution vectors can be total in linear or incremental in nonlinear. +% +% g global dof g = m + n +% m dependent dof +% n independent dof n= s + f +% f free dof +% s single point constraint dof +% +% Input: u_f solution vectors (f size) +% in Matlab workspace +% G_mn Reduction matrix from constraints (m x n) +% in Matlab workspace +% y_s vectors of enforced displacements y_s +% in Matlab workspace +% +% Output: u_g solution vectors (g size) +% +% Input from global workspace: +% +% uset uset.pv_* partitioning vectors and sizes +% uset.*size *=n,m etc. +% +% Called by sol101 (ms), Normalmodes (ms), sol159 (ms) (all solution sequences) +% Calls none +% where ms= m-script, mf= m-function, mex= executable function (c-code) +% +% m_set and n_set index vectors are with respect to (wrt) the g- set +% s_set and f_set index vectors are with respect to (wrt) the n- set + +%Merge f_set and s_set to n_set + +ncol= size( u_f, 2); +u_n= sparse(uset.nsize,ncol); +u_n(uset.pv_f,:)= u_f; + + +% Note: we assume u_f and y_s have the same no.of columns +if ~isempty(y_s) + u_n(uset.pv_s,:)= y_s; +end + +% Merge n_set and m_set to g_set + +u_g = sparse(uset.gsize,ncol); +u_g(uset.pv_n,:)= u_n; + + +if (uset.msize>0) + u_m = G_mn * u_n; + u_g(uset.pv_m,:)= u_m; +end + +end Index: /issm/trunk/src/m/solutions/cielo/cielo/NewCreateFEMModel.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/NewCreateFEMModel.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/NewCreateFEMModel.m (revision 31) @@ -0,0 +1,105 @@ +%----------------------------------------------------------------------- +% CreateFEMModel.m: from a matlab ISSM @model class object, create a finite element +% compatible with Cielo. +% +%----------------------------------------------------------------------- +function m=CreateFEMModel(md,analysis_type ) + + global ParameterDefaults + +% read input data: + + disp(sprintf('\n reading data from model %s...',md.name)); + [m.ncase, m.dynamic, m.ept, m.geom1, m.geom1b, m.geom2, m.part, m.tpart, m.geom3, ... + m.geom4, m.mpt, m.pvt, m.sets, m.dit] = newImp( md ,analysis_type); + + error; + +% generate list of states from case control (ncase), if any: + + disp(' generating states...'); + [m.states] = GenerateStates( m.ncase, m.dynamic); + + %ncase not needed anymore + m=rmfield(m,'ncase'); + +% generate grid point definition table +% (internal grid records, contain sil and uset data) + + disp(' generating gpdt...'); + [m.gpdt,m.gpdtb] = Gp1( m.geom1, m.geom1b); + + +% generate coordinate system transformation matrices, transform +% grids to basic coordinate system + + disp(' generating updated gpdt...'); + [m.cstm, m.bgpdt, m.bgpdtb] = Gp1b( m.geom1, m.geom1b, m.gpdt, m.gpdtb); + + m=rmfield(m,'gpdt'); + if ~isempty(m.gpdtb), m=rmfield(md,'gpdtb');end; + + %geom1 and geom1b not needed anymore + m=rmfield(m,'geom1'); + if ~isempty(m.geom1b), m=rmfield(m,'geom1b'); end; + + %disp(' generating element objects...'); + [m.est,m.lst] = Gp2( m.bgpdt, m.bgpdtb, m.ept, m.geom2, m.mpt, m.geom3); + + %geom2 not needed anymore + m=rmfield(m,'geom2'); + +% phase 1: +% + +% set parameters using bulk data, nastran case control above subcase level: + current_state = {}; + [next_state] = NextState( m.states, current_state); + m.current_state = next_state; + clear next_state current_state; + + [m.params] = Pvt( m.current_state, m.pvt); + + +% generate usets and constraint matrices: + + disp(' generating uset data and constraints...'); + + [m.updated_bgpdt, m.updated_bgpdtb, m.y_g, m.R_mg] = Gp4( m.states, m.geom4, m.bgpdt, m.bgpdtb,analysis_type); + + + %geom4 not needed anymore + if ~isempty(m.geom4), m=rmfield(m,'geom4');end; + +% replace old bgpdt with new (bgpdt unchanged if nonexistent updated_bgpdt): + m.bgpdt = m.updated_bgpdt;m=rmfield(m,'updated_bgpdt'); + if ~isempty(m.updated_bgpdtb), + m.bgpdtb = m.updated_bgpdtb; + m=rmfield(m,'updated_bgpdtb'); + end; + +% create uset table + m.uset=Builduset(m.bgpdt,m.bgpdtb,analysis_type); + + +% +% phase 1: +% +% autospc operations: +% autospc is done on the g_set, it should be done on the n_set later + +% Partition R_mg into R_mm and R_mn, and +% calculate G_mn = - R_mm^-1 * R_mn (MCE1 type operations) + + [m.G_mn] = Normalizeconstraints( m.R_mg, m.uset); + +% Reduce enforced displacements / temperatures + [m.y_s] = Reducevectorg( m.y_g ,m.uset); + + + if isempty(m.y_s), + flag_y_s= 0; + else + flag_y_s= 1; + end +end Index: /issm/trunk/src/m/solutions/cielo/cielo/Orth.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/Orth.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/Orth.m (revision 31) @@ -0,0 +1,30 @@ +function new_gradients=Orth(gradients,oldgradients); + +global client_server_mode IMDataCounter + +%outpu +new_gradients=gradients; + +numparams=length(gradients); +for i=1:numparams, + + %recover gradient #i + gradient=gradients{i}; + oldgradient=oldgradients{i}; + + if isempty(oldgradient), + oldgradient=zeros(length(gradient),1); + end + + % normalize gradient #i to 1 + gradient=gradient/norm(gradient); + + % rough orthagonalization + gradient=gradient-(gradient'*oldgradient)*oldgradient; + + % renormalize gradients + gradient=gradient/norm(gradient); + + %plug back into gradients + new_gradients{i}=gradient; +end Index: /issm/trunk/src/m/solutions/cielo/cielo/Reducematrix.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/Reducematrix.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/Reducematrix.m (revision 31) @@ -0,0 +1,61 @@ +function [K_nn]= reducematrix( K_gg, G_mn, pv_m_set, pv_n_set, msize, flag ); +% REDUCEMATRIX [K_nn]= reducematrix( K_gg, G_mn, pv_m_set, pv_n_set, msize, flag ) +% Reduces K_gg to K_nn +% +% K_nn^~ K_nm +% K_gg = [ ] +% K_mn K_mm +% +% if flag == 1 +% K_nn = K_nn^~ + K_nm * G_mn +% else +% K_nn = K_nn^~ + K_nm * G_mn + G_mn^T * K_mn + G_mn^T * K_mm * G_mn +% +% m dof to be condensed +% n dof to be retained +% g total dof g = m + n +% +% Input: K_gg matrix (g x g) +% G_mn Reduction matrix from constraints (m x n) +% pv_m_set partitioning vector for the m dof wrt g +% pv_n_set partitioning vector for the n dof wrt g +% msize number of m-set dof +% flag flag for type of condensation, see above +% +% Output: K_nn reduced matrix (n x n) +% +% Called by Reducematrixfromgtof (mf), Normalmodes (ms) +% Calls none +% where ms= m-script, mf= m-function, mex= executable function (c-code) +% +% pv_m_set and pv_n_set index vectors are with respect to the g- set. +% +% All input and output variables are assumed to be in the workspace. + +if msize > 0 + +% Partition K_gg + + K_mm = K_gg( pv_m_set, pv_m_set); + K_mn = K_gg( pv_m_set, pv_n_set); + K_nm = K_gg( pv_n_set, pv_m_set); + K_nn = K_gg( pv_n_set, pv_n_set); + +% Reduce K_gg to K_nn + + K_nn = K_nn + K_nm * G_mn; + + if flag ~= 1 + + K_mn = K_mn + K_mm * G_mn; + K_nn = K_nn + G_mn' * K_mn; + + end + +else + +% msize=0 just return input + + K_nn = K_gg; + +end Index: /issm/trunk/src/m/solutions/cielo/cielo/Reducematrixfromgtof.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/Reducematrixfromgtof.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/Reducematrixfromgtof.m (revision 31) @@ -0,0 +1,51 @@ +function [ K_ff, K_fs] = reducematrixfromgtof(K_gg, G_mn, flag_u_s,uset); +% REDUCEMATRIXFROMGTOF [K_ff, K_fs] = Reducematrixfromgtof( K_gg, G_mn, flag_u_s ) +% Reduces a matrix from g-size to f-size +% +% Input: +% K_gg Matrix in g-size +% G_mn Reduction matrix +% flag_u_s flag indicating whether K_fs is calculated or not +% =0 K_fs is not calculated, an empty matrix is returned +% =1 K_fs is calculated and returned +% Output: +% K_ff Reduced system tangent matrix f-size x f-size +% K_fs Reduced sytem tangent matrix f-size x s-size +% +% Input from global workspace: +% +% uset uset.pv_* partitioning vector +% uset.*size size +% *=n,m etc., see Builduset +% +% Called by sol101 (ms), Normalmodes (ms), Tangentmatrixheat (mf) +% Calls Reducematrix (mf) +% where ms= m-script, mf= m-function, mex= executable function (c-code) +% +% The g-set is partitioned as follows +% g = n + m +% n = f + s +% g all system degrees of freedom (dof) +% m multi-point constraint dofs to be condensed out +% n remaining dof after m-set is condensed from the g-set +% s boundary condition dofs (single point constraints) +% f remaining free dofs after s-set is eliminated from the n-set +% + +% Reduce matrix from g-size to n-size + +K_nn= Reducematrix( K_gg, G_mn, uset.pv_m, uset.pv_n, uset.msize, 2 ); + +% Reduce matrix from n-size to f-size + +if uset.ssize>0 + K_ff= K_nn(uset.pv_f, uset.pv_f); + if flag_u_s + K_fs= K_nn(uset.pv_f,uset.pv_s); + else + K_fs=[]; + end +else + K_ff= K_nn; + K_fs=[]; +end Index: /issm/trunk/src/m/solutions/cielo/cielo/Reducerightside.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/Reducerightside.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/Reducerightside.m (revision 31) @@ -0,0 +1,111 @@ +function [p_f]= reducerightside( p_g, G_mn, K_fs, y_s, flag_y_s, uset ) +% REDUCERIGHTSIDE [p_f]= Reducerightside( p_g, G_mn, K_fs, y_s, flag_y_s ) +% Reduces right hand side vectors from g-size to f-size (p_g to p_f). +% +% Partition p_g into p_n and p_m, then modify p_n +% +% p_n = p_n + G_mn^T * p_m +% +% where G_mn is the reduction matrix from multi-point constraints +% +% Partition p_n into p_f and p_s, then modify p_f if flag_y_s > 0 +% +% p_f = p_f - K_fs * y_s +% +% where y_s are enforced displacements/temperatures, +% +% The g-set is partitioned as follows +% g = n + m +% n = f + s +% g all system degrees of freedom (dof) +% m multi-point constraint dofs to be condensed out +% n remaining dof after m-set is condensed from the g-set +% s boundary condition dofs (single point constraints) +% f remaining free dofs after s-set is eliminated from the n-set +% +% Input: p_g right hand side vectors (g size) +% G_mn Reduction matrix from constraints (m x n) +% K_fs Partitioned stiffness matrix (f x s) +% y_s enforced displacements / temperatures s size +% flag_y_s flag for enforced displacement / temperature +% > 0 K_fs * y_s is subtracted +% <= 0 no action +% +% Output: p_f reduced right hand side vectors (f-size) +% +% Called by sol101 (ms), Normalmodes (ms), Residualheat (mf) +% Calls none +% where ms= m-script, mf= m-function, mex= executable function (c-code) +% +% All input and output variables are assumed to be in the Matlab workspace. +% The load vectors p_g and y_s must have the same number of columns. +% +% Input from global workspace: +% +% uset uset.pv_* partitioning vector +% uset.*size size +% *=n,m etc., see Builduset +% +% uset.pv_m and uset.pv_n index vectors are wrt the g- set +% uset.pv_s and uset.pv_f index vectors are wrt the n- set +% etc. + +% Reduce p_g to p_f if not empty + +if ~isempty(p_g) + + % Reduce p_g to p_n + + if ( uset.msize > 0 ) + p_n = p_g( uset.pv_n, :); + p_m = p_g( uset.pv_m, :); + p_n = p_n + G_mn' * p_m; + else + p_n = p_g; + end + + % Reduce p_n to p_f + + if ( uset.ssize > 0 ) + p_f = p_n( uset.pv_f, :); + else + p_f = p_n; + end + +% Create a p_f if p_g is empty + +else + + if isempty(y_s) + disp('No right hand side found...'); + p_f=[]; + return; + else + ncols= size( y_s, 2); + p_f= sparse([],[],[],uset.fsize, ncols); + end + +end + +% for nonzero boundary conditions, subtract coupling forces, +% this operation is only executed if flag_y_s > 0, +% the flag is set in the calling script for linear analysis +% and the first iteration step in nonlinear analysis + +if ( flag_y_s > 0 ) & ( ~isempty(y_s) ) + % check if p_f and y_s have the same number of columns + + ncolp= size(p_f,2); + ncoly= size(y_s,2); + + if ncolp ~= ncoly + disp('Right hand side can not be calculated...'); + return; + end + + p_f = p_f - K_fs * y_s; + +end + + +end Index: /issm/trunk/src/m/solutions/cielo/cielo/Reducevectorg.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/Reducevectorg.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/Reducevectorg.m (revision 31) @@ -0,0 +1,19 @@ +function [y_s]= Reducevectorg( y_g,uset ) +% REDUCEVECTORG [y_s ]= reducevectorg( py_g, uset ) +% +% +% Input in global workspace: +% +% uset uset.pv_* partitioning vectors and sizes +% uset.*size *=n,m etc. + +if ( uset.ssize>0 ) + if ~isempty(y_g), + y_n= y_g(uset.pv_n,:); + y_s= y_n(uset.pv_s,:); + else + y_s=[]; + end +else + y_s=[]; +end Index: /issm/trunk/src/m/solutions/cielo/cielo/Residualice.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/Residualice.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/Residualice.m (revision 31) @@ -0,0 +1,3 @@ +function [res_g] = Residualice( KH, pnl, u); + res_g = pnl - KH * u; +end Index: /issm/trunk/src/m/solutions/cielo/cielo/SetParameterDefaults.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/SetParameterDefaults.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/SetParameterDefaults.m (revision 31) @@ -0,0 +1,46 @@ +%SETPARAMETERDEFAULTS +% SetParameterDefaults sets parameter defaults for solution sequences. +% In particular, Emg module defaults, as well a Client/Server parameters. +% +% See also Pvt + + +global ParameterDefaults; + +ParameterDefaults.autospc ='yes'; +ParameterDefaults.coupmass = -1; % translated to lumped/consistent option +ParameterDefaults.epzero = 1.e-8; +ParameterDefaults.sigma = 0.; +ParameterDefaults.tabs = 0.; + +% for testing, may be made permanent later: + +ParameterDefaults.linesrch = 1; % 0=no, 1=yes + +% Client/Server mode settings + +ParameterDefaults.tunnel='no'; %enable ssh tunnelling for encrypted communication to Server. +ParameterDefaults.batsiz=1; %size of batch jobs sent to nodes +ParameterDefaults.fdbckout=0; %feedback output level from server (0: no output, 1: basic output, 3: large output) +ParameterDefaults.batch='no'; %run in batch mode. +ParameterDefaults.errchk='yes'; %run with error chechking. + +% Integration rule for sol159 is set to 3 (mid point rule) + + ParameterDefaults.int_rule = 3; + + +% John's Emg.mex, and Plotresults.m are picking up the following parameters + + ParameterDefaults.analysis='structural'; + ParameterDefaults.zgrad ='no'; + ParameterDefaults.imoselem='no'; + ParameterDefaults.wtmass = 1.; + ParameterDefaults.lumpm = 1; + ParameterDefaults.mred = 0; + ParameterDefaults.niref = 1; + ParameterDefaults.njref = 1; + ParameterDefaults.nkref = 1; + ParameterDefaults.cido =-1; + ParameterDefaults.solverstring = ' -mat_type aij -ksp_type cgs -pc_type asm -sub_pc_type lu '; + Index: /issm/trunk/src/m/solutions/cielo/cielo/SwitchPartitioning.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/SwitchPartitioning.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/SwitchPartitioning.m (revision 31) @@ -0,0 +1,30 @@ +function new_u_g=SwitchPartitioning(u_g,destination,part,tpart,doflist); +%INPUT function new_u_g=SwitchPartitioning(u_g,destination,part,tpart,doflist); + +gsize=length(u_g); +if strcmpi(destination,'workspace'), + %g-size vector u_g is ordered according to a parallel partition. REorder it + %so that u_g correspondgs to the workspace partition (ie: the initial grid + %partition) + new_u_g=zeros(gsize,1); + + %go through list of degrees of freedom that need repartitioning + for i=1:length(doflist), + dof=doflist(i); + ind=dof:6:gsize; ind=ind(tpart); + new_u_g(dof:6:gsize)=u_g(ind); + end +elseif strcmpi(destination,'cluster'), + %g size vector u_g is ordered according to the matlab workspace. We need to + %reorder it so that u_g is now ordered in the parallel partition. + + new_u_g=zeros(gsize,1); + + for i=1:length(doflist), + dof=doflist(i); + ind=dof:6:gsize; vector=u_g(ind); + new_u_g(dof:6:gsize)=vector(part); + end +else + error('SwitchPartitioning error message: only types of destination allowed are ''workspace'' or ''cluster''!'); +end Index: /issm/trunk/src/m/solutions/cielo/cielo/cielodiagnostic_core_linear.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cielodiagnostic_core_linear.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cielodiagnostic_core_linear.m (revision 31) @@ -0,0 +1,48 @@ +function [ru_g m]=cielodiagnostic_core_linear(m,params,inputs,analysis) +%INPUT ru_g =cielodiagnostic_core_linear(m,params,inputs,analysis) + + %check inputs + if isempty(inputs), + clear inputs; + end + + %recover flag_y_s from m.y_s + if isempty(m.y_s), + flag_y_s= 0; + else + flag_y_s= 1; + end + + %stiffness and load generation only: + params.kflag=1; params.pflag=1; + params.ktflag=0; + + %system matrices + [rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.lst,m.ept,m.mpt,m.geom3,params,inputs,analysis); + + %Reduce tangent matrix from g size to f size + [rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); + + %Cleanup memory + rK_gg=IMdb('drop rK_gg'); + if ~isempty(rdK_gg), rdK_gg=IMdb('drop rdK_gg');end; + + %Reduce load from g size to f size + [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); + + %Cleanup memory + [rp_g rK_fs]=IMdb('drop rp_g rK_fs'); + + %Solve + [ru_f]=Solver(rK_ff,rp_f,{},params.solverstring); + + %Cleanup memory + [rK_ff rp_f]=IMdb('drop rK_ff rp_f'); + + %Merge back to g set + ru_g= Mergesolvec( ru_f, m.G_mn, m.y_s ); + + %Cleanup memory + ru_f=IMdb('drop ru_f '); + +end %end function Index: /issm/trunk/src/m/solutions/cielo/cielo/cielomelting_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cielomelting_core.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cielomelting_core.m (revision 31) @@ -0,0 +1,50 @@ +function rmelting_g =cielomelting_core_linear(m,params,thermalparams,inputs,analysis) +%INPUT function rmelting_g =cielomelting_core_linear(m,params,thermalparams,inputs,analysis) + + %recover flag_y_s from m.y_s + if m.y_s.header.M==0, + flag_y_s= 0; + else + flag_y_s= 1; + end + + %stiffness and load generation only: + params.kflag=1; params.pflag=1; + params.ktflag=0; + + %optional inputs + inputs.beta=thermalparams.beta; + inputs.meltingpoint=thermalparams.meltingpoint; + inputs.latentheat=thermalparams.latentheat; + inputs.heatcapacity=thermalparams.heatcapacity; + inputs.penalty=thermalparams.penalty; + + %system matrices + [rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.lst,m.ept,m.mpt,m.geom3,params,inputs,analysis); + + %Reduce tangent matrix from g size to f size + [rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); + + if params.lowmem, + [rK_gg,rdK_gg]=IMdb('drop rK_gg rdK_gg'); + end + + %Reduce load from g size to f size + [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); + if params.lowmem, + [rp_g rK_fs]=IMdb('drop rp_g rK_fs'); + end + + %Solve + [rt_f]=Solver(rK_ff,rp_f,{},params.solverstring); + + if params.lowmem, + [rK_ff rp_f]=IMdb('drop rK_ff rp_f'); + end + + %Merge back to g set + rmelting_g= Mergesolvec( rt_f, m.G_mn, m.y_s ); + + if (params.lowmem && count>2), + rt_f=IMdb('drop rt_f '); + end Index: /issm/trunk/src/m/solutions/cielo/cielo/cielomesh.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cielomesh.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cielomesh.m (revision 31) @@ -0,0 +1,20 @@ +function md=cielomesh(md) +%INPUT function md=cielomesh(md) +%mesh solution + + %timing + t1=clock; + + analysis='mesh'; + + %configure, and initialise data model and parameter defaults: + SetParameterDefaults; + + %first build an initial triangulation, using the domain outline and qhull. + [md.elements,md.x,md.y,md.z]=MeshTriangulate(md); + md.numberofgrids=length(md.x); + md.numberofelements=length(md.elements); + + %timing + t2=clock; + disp(sprintf('\n%s\n',[' solution lasted ' num2str(etime(t2,t1)) ' seconds'])); Index: /issm/trunk/src/m/solutions/cielo/cielo/cielomesh2grid.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cielomesh2grid.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cielomesh2grid.m (revision 31) @@ -0,0 +1,37 @@ +function md=cielomesh2grid(md) +%INPUT function md=cielomesh2grid(md) +%mesh2grid solution + + %timing + t1=clock; + + analysis='mesh2grid'; + + %configure, and initialise data model and parameter defaults: + SetParameterDefaults; + + %Initialize results + md.mesh2grid_results=cell(length(md.mesh2grid_parameters),1); + + %Recover domain outline as a contour structure. + writefile('contours.exp',md.domainoutline); + contours=expread('contours.exp',0); + contours=contours(1); %just keep the outer domain outline. + delete('contours.exp'); + + %loop on the number of mesh2grid parameters: + for i=1:length(md.mesh2grid_parameters), + parameter=eval(['md.' md.mesh2grid_parameters{i}]); + + %call GriddataMeshToGrid module: + md.mesh2grid_results{i}=... + GriddataMeshToGrid( md.elements,md.x,md.y,contours,parameter,md.mesh2grid_cornereast,md.mesh2grid_cornernorth,... + md.mesh2grid_xposting,md.mesh2grid_yposting,md.mesh2grid_nlines,md.mesh2grid_ncolumns,md.mesh2grid_interpolation{i},... + md.mesh2grid_filter{i},md.mesh2grid_windowsize); + + end + + + %timing + t2=clock; + disp(sprintf('\n%s\n',[' solution lasted ' num2str(etime(t2,t1)) ' seconds'])); Index: /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic.m (revision 31) @@ -0,0 +1,51 @@ +function md=cieloprognostic(md) +%diagnostic solution for horizontal velocity + + %timing + t1=clock; + + analysis_p='prognostic'; + + %configure, and initialise data model and parameter defaults: + SetParameterDefaults; + + %Build model requested for prognostic simulation + m=CreateFEMModel(md,analysis_p); + + % figure out number of dof: just for information purposes. + md.dof=m.uset.fsize; + + %Build partitioning vectors. We only need one partitioning for all our models, we made sure of it in Imp! + indx=1:6:m.uset.gsize; indx=indx(tpart); + + %build some parameters needed in core solution. + m.rifts=md.rifts; + + %compute depth averaged horizontal velocity in Matlab workspace partitioning + u_g=zeros(m.uset.gsize,1); u_g(1:6:end)=md.vx/md.yts; u_g(2:6:end)=md.vy/md.yts; + velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g); + + %move surface and bed velocities to first lower layer, in workspace partitioning + ws=ShiftLayers(md,md.vz/md.yts,md.numlayers,1); + wb=ShiftLayers(md,md.vz/md.yts,1,1); + melting=ShiftLayers(md,md.melting,1,1); + accumulation=ShiftLayers(md,md.accumulation,md.numlayers,1); + + %plug into g set, with cluster partitioning + ws_g=zeros(m.uset.gsize,1);ws_g(1:6:m.uset.gsize)=ws(part); + wb_g=zeros(m.uset.gsize,1);wb_g(1:6:m.uset.gsize)=wb(part); + melting_g=zeros(m.uset.gsize,1);melting_g(1:6:m.uset.gsize)=melting(part); + accumulation_g=zeros(m.uset.gsize,1);accumulation_g(1:6:m.uset.gsize)=accumulation(part); + + %initialize thickness + thickness=zeros(m.uset.gsize,1);thickness(1:6:end)=md.thickness(part); + + %prepare inputs strcuture + inputs=struct('thickness',thickness,'melting',melting_g,'accumulation',accumulation_g,... + 'surface_vertical_velocity',ws_g,'basal_vertical_velocity',wb_g,'dt',md.dt,'velocity_average',velocity_average); + + %Run core solution + h_g=cieloprognostic_core(m,m.params,inputs,analysis_p); + + %recover thickness and project onto every grid of the mesh + md.new_thickness=project3d(md,project2d(md,h_g(indx),1),'node'); Index: /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic_core.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic_core.m (revision 31) @@ -0,0 +1,35 @@ +function h_g =cieloprognostic_core_linear(m,params,inputs,analysis) +%INPUT function rh_g =cieloprognostic_core_linear(m,params,inputs,analysis) + + %recover flag_y_s from m.y_s + if size(m.y_s,1)==0, + flag_y_s= 0; + else + flag_y_s= 1; + end + + %stiffness and load generation only: + params.kflag=1; params.pflag=1; + params.ktflag=0; + + %system matrices + [K_gg , p_g , dK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.lst,m.ept,m.mpt,m.geom3,params,inputs,analysis); + + %Reduce tangent matrix from g size to f size + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.G_mn, flag_y_s ); + clear K_gg dK_gg + + %Reduce load from g size to f size + [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); + [clear p_g K_fs + + %Solve + [h_f]=Solver(K_ff,p_f,params.solverstring); + clear K_ff p_f + + %Merge back to g set + h_g= Mergesolvec( h_f, m.G_mn, m.y_s, m.uset); + + if count>2, + clear h_f + end Index: /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic_func.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic_func.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cieloprognostic_func.m (revision 31) @@ -0,0 +1,62 @@ +function [rh_g]=cieloprognostic_func(m,params,inputs) +%INPUT function [rh_g]=cieloprognostic_func(m,params,inputs) + + analysis='prognostic'; + + %Check arguments + if nargin==0, + cieloprognostic_func_usage; + error(' '); + end + + options_string=''; + + %recover flag_y_s from m.y_s + if m.y_s.header.M==0, + flag_y_s= 0; + else + flag_y_s= 1; + end + + % stiffness and load generation only: + kflag=1; pflag=1; ktflag=0; + + % disp('stiffness, damping, and mass matrix generation...'); + [rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.ept,m.mpt,m.geom3,kflag,params.sparsity,params.penalty,pflag,ktflag, {},{},inputs,analysis); + + % Reduce tangent matrix from g size to f size + + [rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); + + [rK_gg,rdK_gg]=IMdb('drop rK_gg rdK_gg'); + + % Reduce load from g size to f size + + [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); + + [rp_g rK_fs]=IMdb('drop rp_g rK_fs'); + + + % Solve + if(uset.fsize>0) + if strcmp(client_server_mode,'yes'), + %options_string='-ksp_type gmres or cg -pc_type bjacobi or jacobi -ksp_view -pc_factor_zeropivot 1e-30 -ksp_singmonitor '; + %options_string='-ksp_type cg -ksp_singmonitor -pc_type bjacobi -ksp_rtol 1e-11 -ksp_singmonitor '; + options_string='-mat_type aijmumps -mat_mumps_sym 1 -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 40'; + end + rh_f=Solver(rK_ff,rp_f,{},options_string); + else + disp('All dof are constrained, f_set is empty...'); + end + [rK_ff rp_f]=IMdb('drop rK_ff rp_f'); + + + % Merge back to g set + rh_g= Mergesolvec( rh_f, m.G_mn, m.y_s ); + [rh_f]=IMdb('drop rh_f '); + +function cieloprognostic_func_usage(), + +disp('function [rh_g]=cieloprognostic_func(m,params,inputs)'); +disp('%INPUT function [rh_g]=cieloprognostic_func(m,params,inputs)'); + Index: /issm/trunk/src/m/solutions/cielo/cielo/cielothermal.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cielothermal.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cielothermal.m (revision 31) @@ -0,0 +1,98 @@ +function md=cielothermal(md,solutiontype) +%thermal solution: steady state and transient + + analysis_t=solutiontype; + analysis_m='melting'; + + %configure, and initialise data model and parameter defaults: + SetParameterDefaults; + + %Build thermal and melting models + m_t=CreateFEMModel(md,analysis_t); + m_m=CreateFEMModel(md,analysis_m); + + + %Build partitioning vectors. + indx=1:6:mt.uset.gsize; indx=indx(tpart); + + % figure out number of dof: just for information purposes. + md.dof=m_t.uset.fsize; %biggest dof number + + %initialize velocity,pressure and material properties + velocity=zeros(m_t.uset.gsize,1); + velocity(1:6:m_t.uset.gsize)=md.vx(part)/md.yts; + velocity(2:6:m_t.uset.gsize)=md.vy(part)/md.yts; + velocity(3:6:m_t.uset.gsize)=md.vz(part)/md.yts; + + pressure=zeros(m_t.uset.gsize,1); + pressure(1:6:m_t.uset.gsize)=md.pressure(part); + + thermalparams.beta=md.beta; + thermalparams.meltingpoint=md.meltingpoint; + thermalparams.latentheat=md.latentheat; + thermalparams.heatcapacity=md.heatcapacity; + thermalparams.penalty=md.penalty_thermal; + thermalparams.min_thermal_constraints=md.min_thermal_constraints; + thermalparams.mintemp=min(md.dirichletvalues_thermal(find(md.gridondirichlet_thermal))); + + if strcmpi(analysis_t,'thermalsteady'), + + inputs=struct('pressure',pressure,'velocity',velocity,'dt',md.dt); + + %Call core thermal computation + disp(sprintf('\n%s\n','computing temperature...')); + [t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t); + + %Call core melting computation: + %disp(sprintf('\n%s\n','computing melting...')); + %SetUset(m_m); + %inputs=struct('pressure',pressure,'temperature',rt_g,'dt',md.dt); + %rmelting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m); + + %load results onto model + md.temperature=t_g(indx); + %melting_g=IMdb('select matrix from rmelting_g'); + %md.melting=melting_g(indx)*md.yts; %from m/s to m/a + else + + nsteps=md.ndt/md.dt; + + %initialize temperature and melting + temperature=zeros(m_t.uset.gsize,1); + temperature(1:6:m_t.uset.gsize)=md.temperature(part); + soln.t_g=temperature; + %melting=zeros(uset.gsize,1); + %melting(1:6:uset.gsize)=md.melting(part)/md.yts; %in m/s + %soln.melting_g=melting; + + for n=1:nsteps, + + disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt)); + %Call core thermal computation + disp(' computing temperature...'); + inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'velocity',velocity,'dt',md.dt); + [soln(n+1).t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t); + + %Call core melting computation + %disp(' computing melting...'); + %SetUset(m_m); + %inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'dt',md.dt); + %rmelting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m); + %soln(n+1).melting_g=IMdb('select matrix from rmelting_g'); + + end + + %Wrap up + solution_temperature=struct('temperature',{}); + for n=1:nsteps+1, + solution_temperature(n).temperature=soln(n).t_g(indx); + end + md.temperature=solution_temperature; + + %solution_melting=struct('melting',{}); + %for n=1:nsteps+1, + % solution_melting(n).melting=soln(n).melting_g(indx)*md.yts; %in m/a + %end + %md.melting=solution_melting; + end + Index: /issm/trunk/src/m/solutions/cielo/cielo/cielothermal_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielo/cielothermal_core.m (revision 31) +++ /issm/trunk/src/m/solutions/cielo/cielo/cielothermal_core.m (revision 31) @@ -0,0 +1,94 @@ +function [rt_g m]=cielothermal_core(m,params,thermalparams,inputs,analysis) +%INPUT function rt_g =cielothermal_core(m,params,thermalparams,inputs,analysis) + + %check inputs + if isempty(inputs), + clear inputs; + end + + %recover flag_y_s from m.y_s + if m.y_s.header.M==0, + flag_y_s= 0; + else + flag_y_s= 1; + end + + converged=0; + %stiffness and load generation only: + params.kflag=1; params.pflag=1; + params.ktflag=0; + + dt=inputs.dt; + + disp(sprintf('%s',' starting direct shooting method')); + + %Recover pressure from inputs + pressure=inputs.pressure; + + %Optional inputs + inputs.beta=thermalparams.beta; + inputs.meltingpoint=thermalparams.meltingpoint; + inputs.mintemp=thermalparams.mintemp; + + %intialize lst + lst=m.lst; + + count=1; + + while(~converged), + + %system matrices + [rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,lst,m.ept,m.mpt,m.geom3,params,inputs,analysis); + + %Reduce tangent matrix from g size to f size + [rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); + + %Cleanup memory + rK_gg=IMdb('drop rK_gg'); + if ~isempty(rdK_gg), rdK_gg=IMdb('drop rdK_gg');end; + + %Reduce load from g size to f size + [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); + + %Cleanup memory + [rp_g rK_fs]=IMdb('drop rp_g rK_fs'); + + %Solve + [rt_f]=Solver(rK_ff,rp_f,params.solverstring); + + + %Cleanup memory + [rK_ff rp_f]=IMdb('drop rK_ff rp_f'); + + %Merge back to g set + rt_g= Mergesolvec( rt_f, m.G_mn, m.y_s ); + + %Cleanup memory + if (count>2), + rt_f=IMdb('drop rt_f '); + end + + %Deal with penalty loads + [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); + + if ~converged, + disp(sprintf(' %s %i','#unstable constraints ',num_unstable_constraints)); + if num_unstable_constraints2), - rt_f=IMdb('drop rt_f '); - end Index: sm/trunk/src/m/solutions/cielo/cielomesh.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielomesh.m (revision 30) +++ (revision ) @@ -1,20 +1,0 @@ -function md=cielomesh(md) -%INPUT function md=cielomesh(md) -%mesh solution - - %timing - t1=clock; - - analysis='mesh'; - - %configure, and initialise data model and parameter defaults: - SetParameterDefaults; - - %first build an initial triangulation, using the domain outline and qhull. - [md.elements,md.x,md.y,md.z]=MeshTriangulate(md); - md.numberofgrids=length(md.x); - md.numberofelements=length(md.elements); - - %timing - t2=clock; - disp(sprintf('\n%s\n',[' solution lasted ' num2str(etime(t2,t1)) ' seconds'])); Index: sm/trunk/src/m/solutions/cielo/cielomesh2grid.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielomesh2grid.m (revision 30) +++ (revision ) @@ -1,37 +1,0 @@ -function md=cielomesh2grid(md) -%INPUT function md=cielomesh2grid(md) -%mesh2grid solution - - %timing - t1=clock; - - analysis='mesh2grid'; - - %configure, and initialise data model and parameter defaults: - SetParameterDefaults; - - %Initialize results - md.mesh2grid_results=cell(length(md.mesh2grid_parameters),1); - - %Recover domain outline as a contour structure. - writefile('contours.exp',md.domainoutline); - contours=expread('contours.exp',0); - contours=contours(1); %just keep the outer domain outline. - delete('contours.exp'); - - %loop on the number of mesh2grid parameters: - for i=1:length(md.mesh2grid_parameters), - parameter=eval(['md.' md.mesh2grid_parameters{i}]); - - %call GriddataMeshToGrid module: - md.mesh2grid_results{i}=... - GriddataMeshToGrid( md.elements,md.x,md.y,contours,parameter,md.mesh2grid_cornereast,md.mesh2grid_cornernorth,... - md.mesh2grid_xposting,md.mesh2grid_yposting,md.mesh2grid_nlines,md.mesh2grid_ncolumns,md.mesh2grid_interpolation{i},... - md.mesh2grid_filter{i},md.mesh2grid_windowsize); - - end - - - %timing - t2=clock; - disp(sprintf('\n%s\n',[' solution lasted ' num2str(etime(t2,t1)) ' seconds'])); Index: sm/trunk/src/m/solutions/cielo/cieloprognostic.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cieloprognostic.m (revision 30) +++ (revision ) @@ -1,51 +1,0 @@ -function md=cieloprognostic(md) -%diagnostic solution for horizontal velocity - - %timing - t1=clock; - - analysis_p='prognostic'; - - %configure, and initialise data model and parameter defaults: - SetParameterDefaults; - - %Build model requested for prognostic simulation - m=CreateFEMModel(md,analysis_p); - - % figure out number of dof: just for information purposes. - md.dof=m.uset.fsize; - - %Build partitioning vectors. We only need one partitioning for all our models, we made sure of it in Imp! - indx=1:6:m.uset.gsize; indx=indx(tpart); - - %build some parameters needed in core solution. - m.rifts=md.rifts; - - %compute depth averaged horizontal velocity in Matlab workspace partitioning - u_g=zeros(m.uset.gsize,1); u_g(1:6:end)=md.vx/md.yts; u_g(2:6:end)=md.vy/md.yts; - velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g); - - %move surface and bed velocities to first lower layer, in workspace partitioning - ws=ShiftLayers(md,md.vz/md.yts,md.numlayers,1); - wb=ShiftLayers(md,md.vz/md.yts,1,1); - melting=ShiftLayers(md,md.melting,1,1); - accumulation=ShiftLayers(md,md.accumulation,md.numlayers,1); - - %plug into g set, with cluster partitioning - ws_g=zeros(m.uset.gsize,1);ws_g(1:6:m.uset.gsize)=ws(part); - wb_g=zeros(m.uset.gsize,1);wb_g(1:6:m.uset.gsize)=wb(part); - melting_g=zeros(m.uset.gsize,1);melting_g(1:6:m.uset.gsize)=melting(part); - accumulation_g=zeros(m.uset.gsize,1);accumulation_g(1:6:m.uset.gsize)=accumulation(part); - - %initialize thickness - thickness=zeros(m.uset.gsize,1);thickness(1:6:end)=md.thickness(part); - - %prepare inputs strcuture - inputs=struct('thickness',thickness,'melting',melting_g,'accumulation',accumulation_g,... - 'surface_vertical_velocity',ws_g,'basal_vertical_velocity',wb_g,'dt',md.dt,'velocity_average',velocity_average); - - %Run core solution - h_g=cieloprognostic_core(m,m.params,inputs,analysis_p); - - %recover thickness and project onto every grid of the mesh - md.new_thickness=project3d(md,project2d(md,h_g(indx),1),'node'); Index: sm/trunk/src/m/solutions/cielo/cieloprognostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cieloprognostic_core.m (revision 30) +++ (revision ) @@ -1,35 +1,0 @@ -function h_g =cieloprognostic_core_linear(m,params,inputs,analysis) -%INPUT function rh_g =cieloprognostic_core_linear(m,params,inputs,analysis) - - %recover flag_y_s from m.y_s - if size(m.y_s,1)==0, - flag_y_s= 0; - else - flag_y_s= 1; - end - - %stiffness and load generation only: - params.kflag=1; params.pflag=1; - params.ktflag=0; - - %system matrices - [K_gg , p_g , dK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.lst,m.ept,m.mpt,m.geom3,params,inputs,analysis); - - %Reduce tangent matrix from g size to f size - [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.G_mn, flag_y_s ); - clear K_gg dK_gg - - %Reduce load from g size to f size - [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); - [clear p_g K_fs - - %Solve - [h_f]=Solver(K_ff,p_f,params.solverstring); - clear K_ff p_f - - %Merge back to g set - h_g= Mergesolvec( h_f, m.G_mn, m.y_s, m.uset); - - if count>2, - clear h_f - end Index: sm/trunk/src/m/solutions/cielo/cieloprognostic_func.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cieloprognostic_func.m (revision 30) +++ (revision ) @@ -1,62 +1,0 @@ -function [rh_g]=cieloprognostic_func(m,params,inputs) -%INPUT function [rh_g]=cieloprognostic_func(m,params,inputs) - - analysis='prognostic'; - - %Check arguments - if nargin==0, - cieloprognostic_func_usage; - error(' '); - end - - options_string=''; - - %recover flag_y_s from m.y_s - if m.y_s.header.M==0, - flag_y_s= 0; - else - flag_y_s= 1; - end - - % stiffness and load generation only: - kflag=1; pflag=1; ktflag=0; - - % disp('stiffness, damping, and mass matrix generation...'); - [rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.ept,m.mpt,m.geom3,kflag,params.sparsity,params.penalty,pflag,ktflag, {},{},inputs,analysis); - - % Reduce tangent matrix from g size to f size - - [rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); - - [rK_gg,rdK_gg]=IMdb('drop rK_gg rdK_gg'); - - % Reduce load from g size to f size - - [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); - - [rp_g rK_fs]=IMdb('drop rp_g rK_fs'); - - - % Solve - if(uset.fsize>0) - if strcmp(client_server_mode,'yes'), - %options_string='-ksp_type gmres or cg -pc_type bjacobi or jacobi -ksp_view -pc_factor_zeropivot 1e-30 -ksp_singmonitor '; - %options_string='-ksp_type cg -ksp_singmonitor -pc_type bjacobi -ksp_rtol 1e-11 -ksp_singmonitor '; - options_string='-mat_type aijmumps -mat_mumps_sym 1 -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 40'; - end - rh_f=Solver(rK_ff,rp_f,{},options_string); - else - disp('All dof are constrained, f_set is empty...'); - end - [rK_ff rp_f]=IMdb('drop rK_ff rp_f'); - - - % Merge back to g set - rh_g= Mergesolvec( rh_f, m.G_mn, m.y_s ); - [rh_f]=IMdb('drop rh_f '); - -function cieloprognostic_func_usage(), - -disp('function [rh_g]=cieloprognostic_func(m,params,inputs)'); -disp('%INPUT function [rh_g]=cieloprognostic_func(m,params,inputs)'); - Index: sm/trunk/src/m/solutions/cielo/cielothermal.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielothermal.m (revision 30) +++ (revision ) @@ -1,98 +1,0 @@ -function md=cielothermal(md,solutiontype) -%thermal solution: steady state and transient - - analysis_t=solutiontype; - analysis_m='melting'; - - %configure, and initialise data model and parameter defaults: - SetParameterDefaults; - - %Build thermal and melting models - m_t=CreateFEMModel(md,analysis_t); - m_m=CreateFEMModel(md,analysis_m); - - - %Build partitioning vectors. - indx=1:6:mt.uset.gsize; indx=indx(tpart); - - % figure out number of dof: just for information purposes. - md.dof=m_t.uset.fsize; %biggest dof number - - %initialize velocity,pressure and material properties - velocity=zeros(m_t.uset.gsize,1); - velocity(1:6:m_t.uset.gsize)=md.vx(part)/md.yts; - velocity(2:6:m_t.uset.gsize)=md.vy(part)/md.yts; - velocity(3:6:m_t.uset.gsize)=md.vz(part)/md.yts; - - pressure=zeros(m_t.uset.gsize,1); - pressure(1:6:m_t.uset.gsize)=md.pressure(part); - - thermalparams.beta=md.beta; - thermalparams.meltingpoint=md.meltingpoint; - thermalparams.latentheat=md.latentheat; - thermalparams.heatcapacity=md.heatcapacity; - thermalparams.penalty=md.penalty_thermal; - thermalparams.min_thermal_constraints=md.min_thermal_constraints; - thermalparams.mintemp=min(md.dirichletvalues_thermal(find(md.gridondirichlet_thermal))); - - if strcmpi(analysis_t,'thermalsteady'), - - inputs=struct('pressure',pressure,'velocity',velocity,'dt',md.dt); - - %Call core thermal computation - disp(sprintf('\n%s\n','computing temperature...')); - [t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t); - - %Call core melting computation: - %disp(sprintf('\n%s\n','computing melting...')); - %SetUset(m_m); - %inputs=struct('pressure',pressure,'temperature',rt_g,'dt',md.dt); - %rmelting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m); - - %load results onto model - md.temperature=t_g(indx); - %melting_g=IMdb('select matrix from rmelting_g'); - %md.melting=melting_g(indx)*md.yts; %from m/s to m/a - else - - nsteps=md.ndt/md.dt; - - %initialize temperature and melting - temperature=zeros(m_t.uset.gsize,1); - temperature(1:6:m_t.uset.gsize)=md.temperature(part); - soln.t_g=temperature; - %melting=zeros(uset.gsize,1); - %melting(1:6:uset.gsize)=md.melting(part)/md.yts; %in m/s - %soln.melting_g=melting; - - for n=1:nsteps, - - disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt)); - %Call core thermal computation - disp(' computing temperature...'); - inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'velocity',velocity,'dt',md.dt); - [soln(n+1).t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t); - - %Call core melting computation - %disp(' computing melting...'); - %SetUset(m_m); - %inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'dt',md.dt); - %rmelting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m); - %soln(n+1).melting_g=IMdb('select matrix from rmelting_g'); - - end - - %Wrap up - solution_temperature=struct('temperature',{}); - for n=1:nsteps+1, - solution_temperature(n).temperature=soln(n).t_g(indx); - end - md.temperature=solution_temperature; - - %solution_melting=struct('melting',{}); - %for n=1:nsteps+1, - % solution_melting(n).melting=soln(n).melting_g(indx)*md.yts; %in m/a - %end - %md.melting=solution_melting; - end - Index: sm/trunk/src/m/solutions/cielo/cielothermal_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/cielothermal_core.m (revision 30) +++ (revision ) @@ -1,94 +1,0 @@ -function [rt_g m]=cielothermal_core(m,params,thermalparams,inputs,analysis) -%INPUT function rt_g =cielothermal_core(m,params,thermalparams,inputs,analysis) - - %check inputs - if isempty(inputs), - clear inputs; - end - - %recover flag_y_s from m.y_s - if m.y_s.header.M==0, - flag_y_s= 0; - else - flag_y_s= 1; - end - - converged=0; - %stiffness and load generation only: - params.kflag=1; params.pflag=1; - params.ktflag=0; - - dt=inputs.dt; - - disp(sprintf('%s',' starting direct shooting method')); - - %Recover pressure from inputs - pressure=inputs.pressure; - - %Optional inputs - inputs.beta=thermalparams.beta; - inputs.meltingpoint=thermalparams.meltingpoint; - inputs.mintemp=thermalparams.mintemp; - - %intialize lst - lst=m.lst; - - count=1; - - while(~converged), - - %system matrices - [rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,lst,m.ept,m.mpt,m.geom3,params,inputs,analysis); - - %Reduce tangent matrix from g size to f size - [rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); - - %Cleanup memory - rK_gg=IMdb('drop rK_gg'); - if ~isempty(rdK_gg), rdK_gg=IMdb('drop rdK_gg');end; - - %Reduce load from g size to f size - [rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s ); - - %Cleanup memory - [rp_g rK_fs]=IMdb('drop rp_g rK_fs'); - - %Solve - [rt_f]=Solver(rK_ff,rp_f,params.solverstring); - - - %Cleanup memory - [rK_ff rp_f]=IMdb('drop rK_ff rp_f'); - - %Merge back to g set - rt_g= Mergesolvec( rt_f, m.G_mn, m.y_s ); - - %Cleanup memory - if (count>2), - rt_f=IMdb('drop rt_f '); - end - - %Deal with penalty loads - [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); - - if ~converged, - disp(sprintf(' %s %i','#unstable constraints ',num_unstable_constraints)); - if num_unstable_constraints