Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8364) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8365) @@ -21,4 +21,5 @@ Elements *elements = NULL; Materials *materials = NULL; + Parameters *parameters = NULL; /*Create elements, vertices and materials, independent of analysis_type: */ @@ -114,12 +115,12 @@ UpdateElementsAndMaterialsControl(elements,materials,iomodel,iomodel_handle); - /*Update Elements in case we are running a transient solution: */ - if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){ - UpdateElementsTransient(elements,iomodel,iomodel_handle,analysis_counter,analysis_type); - } - - /*Generate objects that are not dependent on any analysis_type: */ CreateParameters(pparameters,iomodel,iomodel_handle,solution_type,analysis_type,analysis_counter); + + /*Update Elements in case we are running a transient solution: */ + parameters=*pparameters; + if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){ + UpdateElementsTransient(elements,parameters,iomodel,iomodel_handle,analysis_counter,analysis_type); + } /*Sort datasets: */ Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 8364) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 8365) @@ -93,5 +93,5 @@ /*transient: */ -void UpdateElementsTransient(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type); +void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type); /*partitioning: */ Index: /issm/trunk/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp (revision 8364) +++ /issm/trunk/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp (revision 8365) @@ -13,50 +13,58 @@ #include "../ModelProcessorx.h" -void UpdateElementsTransient(Elements* elements, IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type){ +void UpdateElementsTransient(Elements* elements, Parameters* parameters,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type){ /*Intermediary*/ - int i,j; + int i,j,k; int counter; Element *element = NULL; char fetchstring[100]; double time; + double forcingenum; + /*if no forcings, bail out: */ + printf("numforcings: %i\n",iomodel->numforcings); + if(!iomodel->numforcings)return; + + /*download whatever data will persist through forcing loop: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - /*download whatever data will persist through forcing processing: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + /*download forcing types: */ + IoModelFetchData(&iomodel->forcingtypes,NULL,NULL,iomodel_handle,"forcingtypes"); - /*Ok, let's go through forcings: */ - /*first, accumulation: */ - if(iomodel->forcing_accumulation_num_time_steps==0)_error_("accumulation forcing not available!"); + /*loop over forcings: */ + for(i=0;inumforcings;i++){ - iomodel->forcing_numtimesteps=iomodel->forcing_accumulation_num_time_steps; - IoModelFetchData(&iomodel->timesteps,NULL,NULL,iomodel_handle,"forcing_accumulation_time_steps"); + forcingenum=iomodel->forcingtypes[i]; + printf("forcing enum: %g\n",iomodel->forcingtypes[i]); + + sprintf(&fetchstring[0],"forcing_%s_num_time_steps",EnumToStringx(forcingenum)); + IoModelFetchData(&iomodel->forcing_numtimesteps,iomodel_handle,fetchstring); + sprintf(&fetchstring[0],"forcing_%s_time_steps",EnumToStringx(forcingenum)); + IoModelFetchData(&iomodel->timesteps,NULL,NULL,iomodel_handle,fetchstring); - for(i=0;iforcing_accumulation_num_time_steps;i++){ - time=iomodel->timesteps[i]; - sprintf(&fetchstring[0],"forcing_accumulation_%i",i+1); - IoModelFetchData(&iomodel->forcing,NULL,NULL,iomodel_handle,fetchstring); + for(j=0;jforcing_numtimesteps;j++){ + + time=iomodel->timesteps[j]; + sprintf(&fetchstring[0],"forcing_%s_%i",EnumToStringx(forcingenum),j+1); + IoModelFetchData(&iomodel->forcing,NULL,NULL,iomodel_handle,fetchstring); - - - /*we now have the forcing for the corresponding time step, for all the nodes. Use this - *to write over the existing accumulation input: */ - - counter=0; - for (j=0;jnumberofelements;j++){ - if(iomodel->my_elements[j]){ - element=(Element*)elements->GetObjectByOffset(counter); - element->UpdateForcing(j,iomodel,i,time,AccumulationRateEnum); //we need j to index into elements. - counter++; - + /*we now have the forcing for the corresponding time step, for all the nodes. Use this + *to write over the existing accumulation input: */ + counter=0; + for (k=0;knumberofelements;k++){ + if(iomodel->my_elements[k]){ + element=(Element*)elements->GetObjectByOffset(counter); + element->UpdateForcing(k,iomodel,j,time,forcingenum,parameters); + counter++; + } } + /*Free ressources:*/ + xfree((void**)&iomodel->forcing); } - - /*Free ressources:*/ - xfree((void**)&iomodel->forcing); + xfree((void**)&iomodel->timesteps); } - /*Free ressources:*/ @@ -64,4 +72,4 @@ xfree((void**)&iomodel->timesteps); xfree((void**)&iomodel->elements); - + xfree((void**)&iomodel->forcingtypes); } Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 8364) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 8365) @@ -59,5 +59,5 @@ virtual void DeleteResults(void)=0; virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; - virtual void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum)=0; + virtual void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters)=0; virtual void InputToResult(int enum_type,int step,double time)=0; virtual void ControlInputGetGradient(Vec gradient,int enum_type)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8364) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8365) @@ -129,5 +129,5 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum){_error_("not supported yet!");} + void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters){_error_("not supported yet!");} int UpdateShelfStatus(Vec new_shelf_nodes); void UpdateShelfFlags(double* new_shelf_nodes); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8364) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8365) @@ -5868,5 +5868,5 @@ /*}}}*/ /*FUNCTION Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){{{1*/ -void Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){ +void Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters){ /*Intermediaries*/ @@ -5880,4 +5880,7 @@ /*recover forcing for this time step: */ for(i=0;i<3;i++)nodeinputs[i]=iomodel->forcing[tria_vertex_ids[i]-1]; + + /*process units: */ + UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, FieldEnum, parameters); if(step==0){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 8364) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 8365) @@ -134,5 +134,5 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum); + void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters); int UpdateShelfStatus(Vec new_shelf_nodes); int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf); Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 8364) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 8365) @@ -89,4 +89,5 @@ xfree((void**)&this->melting_rate_correction); xfree((void**)&this->accumulation_rate); + xfree((void**)&this->forcingtypes); xfree((void**)&this->forcing); xfree((void**)&this->timesteps); @@ -221,5 +222,5 @@ IoModelFetchData(&this->gl_melting_rate,iomodel_handle,"gl_melting_rate"); IoModelFetchData(&this->rheology_law,iomodel_handle,"rheology_law"); - IoModelFetchData(&this->forcing_accumulation_num_time_steps,iomodel_handle,"forcing_accumulation_num_time_steps"); + IoModelFetchData(&this->numforcings,iomodel_handle,"numforcings"); /*qmu: */ @@ -402,9 +403,10 @@ /*forcings: */ - this->forcing_accumulation_num_time_steps=0; + this->numforcings=0; + this->forcingtypes=NULL; this->forcing_numtimesteps=0; this->forcing=NULL; this->timesteps=NULL; - + /*parameter output: */ this->numoutput=0; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 8364) +++ /issm/trunk/src/c/objects/IoModel.h (revision 8365) @@ -198,5 +198,6 @@ /*forcings: */ - int forcing_accumulation_num_time_steps; + int numforcings; + double* forcingtypes; int forcing_numtimesteps; double* forcing; Index: /issm/trunk/src/c/shared/Numerics/UnitConversion.cpp =================================================================== --- /issm/trunk/src/c/shared/Numerics/UnitConversion.cpp (revision 8364) +++ /issm/trunk/src/c/shared/Numerics/UnitConversion.cpp (revision 8365) @@ -69,4 +69,5 @@ case DhDtEnum: scale=yts;break; //m/yr case MeltingRateEnum: scale=yts;break; //m/yr + case AccumulationRateEnum: scale=yts;break; //m/yr case MisfitEnum: scale=pow(yts,2);break; //(m/yr)^2 case MassFluxEnum: scale=pow(10,-12)*yts;break; // (GigaTon/year) Index: /issm/trunk/src/c/solutions/transient2d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 8364) +++ /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 8365) @@ -87,4 +87,5 @@ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AccumulationRateEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,step,time); if(gl_migration!=NoneEnum)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time); Index: /issm/trunk/src/m/classes/model.m =================================================================== --- /issm/trunk/src/m/classes/model.m (revision 8364) +++ /issm/trunk/src/m/classes/model.m (revision 8365) @@ -183,5 +183,5 @@ %Forcings - forcing_accumulation=NaN; + forcings=NaN; %Statics parameters Index: /issm/trunk/src/m/model/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 8364) +++ /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 8365) @@ -289,13 +289,18 @@ end - %Check that forcing has length numberofnodes+1 - fields={'forcing_accumulation'}; - checklength(md,fields,md.numberofnodes+1) + %Check that all forcings have length numberofnodes+1 + forcingnames=fieldnames(md.forcings); + forcingfields={}; + for i=1:length(forcingnames), + forcingfields{end+1}=['forcings.' forcingnames{i}]; + end + checklength(md,forcingfields,md.numberofnodes+1) %Check that forcing columns are properly ordered - if md.forcing_accumulation(end,:)~=sort(md.forcing_accumulation(end,:)), - error(['model not consistent: model ' md.name ' forcing_accumulation field columns should be chronological']); - end - + for i=1:length(forcingnames), + if md.forcings.(forcingnames{i})(end,:)~=sort(md.forcings.(forcingnames{i})(end,:)), + error(['model not consistent: model ' md.name ' forcings.' forcingnames{i} ' columns should be chronological']); + end + end end Index: /issm/trunk/src/m/model/marshall.m =================================================================== --- /issm/trunk/src/m/model/marshall.m (revision 8364) +++ /issm/trunk/src/m/model/marshall.m (revision 8365) @@ -94,12 +94,27 @@ WriteData(fid,md.melting_rate_correction_apply,'Integer','melting_rate_correction_apply'); -if ~isnan(md.forcing_accumulation), - WriteData(fid,size(md.forcing_accumulation,2),'Integer','forcing_accumulation_num_time_steps'); - WriteData(fid,md.forcing_accumulation(end,:),'Mat','forcing_accumulation_time_steps'); - for i=1:size(md.forcing_accumulation,2), - WriteData(fid,md.forcing_accumulation(1:end-1,i),'Mat',['forcing_accumulation_' num2str(i)]); +%deal with forcings +if ~isnans(md.forcings), + forcingnames=fieldnames(md.forcings); + numforcings=length(forcingnames); + forcingtypes=zeros(numforcings,1); + for i=1:numforcings, + forcingtypes(i)=StringToEnum(forcingnames{i}); + end + WriteData(fid,numforcings,'Integer','numforcings'); + WriteData(fid,forcingtypes,'Mat','forcingtypes'); + + for i=1:numforcings, + forcing=md.forcings.(forcingnames{i}); + forcingname=forcingnames{i}; + + WriteData(fid,size(forcing,2),'Integer',['forcing_' forcingname '_num_time_steps']); + WriteData(fid,forcing(end,:),'Mat',['forcing_' forcingname '_time_steps']); + for j=1:size(forcing,2), + WriteData(fid,forcing(1:end-1,j),'Mat',['forcing_' forcingname '_' num2str(j)]); + end end else - WriteData(fid,0,'Integer','forcing_accumulation_num_time_steps'); + WriteData(fid,0,'Integer','numforcings'); end Index: /issm/trunk/src/m/utils/BC/SetIceSheetBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetIceSheetBC.m (revision 8364) +++ /issm/trunk/src/m/utils/BC/SetIceSheetBC.m (revision 8365) @@ -31,5 +31,5 @@ if isnan(md.accumulation_rate), md.accumulation_rate=zeros(md.numberofnodes,1); - md.forcing_accumulation=[zeros(md.numberofnodes+1,1)]; + md.forcings.AccumulationRate=[zeros(md.numberofnodes+1,1)]; disp(' no accumulation_rate specified: values set as zero'); end Index: /issm/trunk/src/m/utils/BC/SetIceShelfBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetIceShelfBC.m (revision 8364) +++ /issm/trunk/src/m/utils/BC/SetIceShelfBC.m (revision 8365) @@ -62,5 +62,5 @@ if isnan(md.accumulation_rate), md.accumulation_rate=zeros(md.numberofnodes,1); - md.forcing_accumulation=zeros(md.numberofnodes+1,1); + md.forcings.AccumulationRate=zeros(md.numberofnodes+1,1); disp(' no accumulation_rate specified: values set as zero'); end Index: /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 8364) +++ /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 8365) @@ -72,5 +72,5 @@ if isnan(md.accumulation_rate), md.accumulation_rate=zeros(md.numberofnodes,1); - md.forcing_accumulation=zeros(md.numberofnodes+1,1); + md.forcings.AccumulationRate=zeros(md.numberofnodes+1,1); disp(' no accumulation_rate specified: values set as zero'); end