Changeset 8365


Ignore:
Timestamp:
05/19/11 21:33:36 (14 years ago)
Author:
Eric.Larour
Message:

trunk: forcings are now completely generic. can specify as many as you want. transient 2d now spits accumulation and melting rates, two most probable forcings

Location:
issm/trunk/src
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp

    r8363 r8365  
    2121        Elements  *elements   = NULL;
    2222        Materials *materials  = NULL;
     23        Parameters *parameters  = NULL;
    2324                       
    2425        /*Create elements, vertices and materials, independent of analysis_type: */
    … …  
    114115        UpdateElementsAndMaterialsControl(elements,materials,iomodel,iomodel_handle);
    115116
    116         /*Update Elements in case we are running a transient solution: */
    117         if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){
    118                 UpdateElementsTransient(elements,iomodel,iomodel_handle,analysis_counter,analysis_type);
    119         }
    120 
    121 
    122117        /*Generate objects that are not dependent on any analysis_type: */
    123118        CreateParameters(pparameters,iomodel,iomodel_handle,solution_type,analysis_type,analysis_counter);
     119
     120        /*Update Elements in case we are running a transient solution: */
     121        parameters=*pparameters;
     122        if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){
     123                UpdateElementsTransient(elements,parameters,iomodel,iomodel_handle,analysis_counter,analysis_type);
     124        }
    124125
    125126        /*Sort datasets: */

  • issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r8330 r8365  
    9393
    9494/*transient: */
    95 void    UpdateElementsTransient(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type);
     95void    UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type);
    9696
    9797/*partitioning: */

  • issm/trunk/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp

    r8363 r8365  
    1313#include "../ModelProcessorx.h"
    1414
    15 void    UpdateElementsTransient(Elements* elements, IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type){
     15void    UpdateElementsTransient(Elements* elements, Parameters* parameters,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type){
    1616
    1717        /*Intermediary*/
    18         int      i,j;
     18        int      i,j,k;
    1919        int      counter;
    2020        Element *element = NULL;
    2121        char     fetchstring[100];
    2222        double   time;
     23        double   forcingenum;
    2324
     25        /*if no forcings, bail out: */
     26        printf("numforcings: %i\n",iomodel->numforcings);
     27        if(!iomodel->numforcings)return;
     28
     29        /*download whatever data will persist through forcing loop: */
     30        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    2431       
    25         /*download whatever data will persist through forcing processing: */
    26         IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
     32        /*download forcing types: */
     33        IoModelFetchData(&iomodel->forcingtypes,NULL,NULL,iomodel_handle,"forcingtypes");
    2734
    28         /*Ok, let's go through forcings: */
    29         /*first, accumulation: */
    30         if(iomodel->forcing_accumulation_num_time_steps==0)_error_("accumulation forcing not available!");
     35        /*loop over forcings: */
     36        for(i=0;i<iomodel->numforcings;i++){
    3137
    32         iomodel->forcing_numtimesteps=iomodel->forcing_accumulation_num_time_steps;
    33         IoModelFetchData(&iomodel->timesteps,NULL,NULL,iomodel_handle,"forcing_accumulation_time_steps");
     38                forcingenum=iomodel->forcingtypes[i];
     39                printf("forcing enum: %g\n",iomodel->forcingtypes[i]);
     40               
     41                sprintf(&fetchstring[0],"forcing_%s_num_time_steps",EnumToStringx(forcingenum));
     42                IoModelFetchData(&iomodel->forcing_numtimesteps,iomodel_handle,fetchstring);
     43                sprintf(&fetchstring[0],"forcing_%s_time_steps",EnumToStringx(forcingenum));
     44                IoModelFetchData(&iomodel->timesteps,NULL,NULL,iomodel_handle,fetchstring);
    3445
    35         for(i=0;i<iomodel->forcing_accumulation_num_time_steps;i++){
    36                 time=iomodel->timesteps[i];
    37                 sprintf(&fetchstring[0],"forcing_accumulation_%i",i+1);
    38                 IoModelFetchData(&iomodel->forcing,NULL,NULL,iomodel_handle,fetchstring);
     46                for(j=0;j<iomodel->forcing_numtimesteps;j++){
     47               
     48                        time=iomodel->timesteps[j];
     49                        sprintf(&fetchstring[0],"forcing_%s_%i",EnumToStringx(forcingenum),j+1);
     50                        IoModelFetchData(&iomodel->forcing,NULL,NULL,iomodel_handle,fetchstring);
    3951
    40 
    41 
    42                 /*we now have the forcing for the corresponding time step, for all the nodes. Use this
    43                  *to write over the existing accumulation input: */
    44 
    45                 counter=0;
    46                 for (j=0;j<iomodel->numberofelements;j++){
    47                         if(iomodel->my_elements[j]){
    48                                 element=(Element*)elements->GetObjectByOffset(counter);
    49                                 element->UpdateForcing(j,iomodel,i,time,AccumulationRateEnum); //we need j to index into elements.
    50                                 counter++;
    51 
     52                        /*we now have the forcing for the corresponding time step, for all the nodes. Use this
     53                         *to write over the existing accumulation input: */
     54                        counter=0;
     55                        for (k=0;k<iomodel->numberofelements;k++){
     56                                if(iomodel->my_elements[k]){
     57                                        element=(Element*)elements->GetObjectByOffset(counter);
     58                                        element->UpdateForcing(k,iomodel,j,time,forcingenum,parameters);
     59                                        counter++;
     60                                }
    5261                        }
     62                        /*Free ressources:*/
     63                        xfree((void**)&iomodel->forcing);
    5364                }
    5465
    55 
    56 
    5766                /*Free ressources:*/
    58                 xfree((void**)&iomodel->forcing);
     67                xfree((void**)&iomodel->timesteps);
    5968        }
    60                
    6169
    6270        /*Free ressources:*/
    … …  
    6472        xfree((void**)&iomodel->timesteps);
    6573        xfree((void**)&iomodel->elements);
    66 
     74        xfree((void**)&iomodel->forcingtypes);
    6775}

  • issm/trunk/src/c/objects/Elements/Element.h

    r8363 r8365  
    5959                virtual void   DeleteResults(void)=0;
    6060                virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0;
    61                 virtual void   UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum)=0;
     61                virtual void   UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters)=0;
    6262                virtual void   InputToResult(int enum_type,int step,double time)=0;
    6363                virtual void   ControlInputGetGradient(Vec gradient,int enum_type)=0;

  • issm/trunk/src/c/objects/Elements/Penta.h

    r8363 r8365  
    129129                double SurfaceArea(void);
    130130                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    131                 void   UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum){_error_("not supported yet!");}
     131                void   UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters){_error_("not supported yet!");}
    132132                int    UpdateShelfStatus(Vec new_shelf_nodes);
    133133                void   UpdateShelfFlags(double* new_shelf_nodes);

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r8363 r8365  
    58685868/*}}}*/
    58695869/*FUNCTION Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){{{1*/
    5870 void Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){
     5870void Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters){
    58715871               
    58725872        /*Intermediaries*/
    … …  
    58805880        /*recover forcing for this time step: */
    58815881        for(i=0;i<3;i++)nodeinputs[i]=iomodel->forcing[tria_vertex_ids[i]-1];
     5882
     5883        /*process units: */
     5884        UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, FieldEnum, parameters);
    58825885
    58835886        if(step==0){

  • issm/trunk/src/c/objects/Elements/Tria.h

    r8363 r8365  
    134134                double SurfaceArea(void);
    135135                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    136                 void   UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum);
     136                void   UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters);
    137137                int    UpdateShelfStatus(Vec new_shelf_nodes);
    138138                int    UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf);

  • issm/trunk/src/c/objects/IoModel.cpp

    r8363 r8365  
    8989        xfree((void**)&this->melting_rate_correction);
    9090        xfree((void**)&this->accumulation_rate);
     91        xfree((void**)&this->forcingtypes);
    9192        xfree((void**)&this->forcing);
    9293        xfree((void**)&this->timesteps);
    … …  
    221222        IoModelFetchData(&this->gl_melting_rate,iomodel_handle,"gl_melting_rate");
    222223        IoModelFetchData(&this->rheology_law,iomodel_handle,"rheology_law");
    223         IoModelFetchData(&this->forcing_accumulation_num_time_steps,iomodel_handle,"forcing_accumulation_num_time_steps");
     224        IoModelFetchData(&this->numforcings,iomodel_handle,"numforcings");
    224225       
    225226        /*qmu: */
    … …  
    402403
    403404        /*forcings: */
    404         this->forcing_accumulation_num_time_steps=0;
     405        this->numforcings=0;
     406        this->forcingtypes=NULL;
    405407        this->forcing_numtimesteps=0;
    406408        this->forcing=NULL;
    407409        this->timesteps=NULL;
    408        
     410
    409411        /*parameter output: */
    410412        this->numoutput=0;

  • issm/trunk/src/c/objects/IoModel.h

    r8363 r8365  
    198198
    199199                /*forcings: */
    200                 int      forcing_accumulation_num_time_steps;
     200                int      numforcings;
     201                double*  forcingtypes;
    201202                int      forcing_numtimesteps;
    202203                double*  forcing;

  • issm/trunk/src/c/shared/Numerics/UnitConversion.cpp

    r6412 r8365  
    6969                case DhDtEnum:        scale=yts;break; //m/yr
    7070                case MeltingRateEnum: scale=yts;break; //m/yr
     71                case AccumulationRateEnum: scale=yts;break; //m/yr
    7172                case MisfitEnum:      scale=pow(yts,2);break; //(m/yr)^2
    7273                case MassFluxEnum:    scale=pow(10,-12)*yts;break; // (GigaTon/year)

  • issm/trunk/src/c/solutions/transient2d_core.cpp

    r8363 r8365  
    8787                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
    8888                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AccumulationRateEnum,step,time);
     89                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,step,time);
    8990                        if(gl_migration!=NoneEnum)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time);
    9091

  • issm/trunk/src/m/classes/model.m

    r8330 r8365  
    183183
    184184                 %Forcings
    185                  forcing_accumulation=NaN;
     185                 forcings=NaN;
    186186
    187187                 %Statics parameters

  • issm/trunk/src/m/model/ismodelselfconsistent.m

    r8352 r8365  
    289289        end
    290290
    291         %Check that forcing has length numberofnodes+1
    292         fields={'forcing_accumulation'};
    293         checklength(md,fields,md.numberofnodes+1)
     291        %Check that all forcings have length numberofnodes+1
     292        forcingnames=fieldnames(md.forcings);
     293        forcingfields={};
     294        for i=1:length(forcingnames),
     295                forcingfields{end+1}=['forcings.' forcingnames{i}];
     296        end
     297        checklength(md,forcingfields,md.numberofnodes+1)
    294298
    295299        %Check that forcing columns are properly ordered
    296         if md.forcing_accumulation(end,:)~=sort(md.forcing_accumulation(end,:)),
    297                 error(['model not consistent: model ' md.name ' forcing_accumulation field columns should be chronological']);
    298         end
    299 
     300        for i=1:length(forcingnames),
     301                if md.forcings.(forcingnames{i})(end,:)~=sort(md.forcings.(forcingnames{i})(end,:)),
     302                        error(['model not consistent: model ' md.name ' forcings.' forcingnames{i} ' columns should be chronological']);
     303                end
     304        end
    300305
    301306end

  • issm/trunk/src/m/model/marshall.m

    r8363 r8365  
    9494WriteData(fid,md.melting_rate_correction_apply,'Integer','melting_rate_correction_apply');
    9595
    96 if ~isnan(md.forcing_accumulation),
    97         WriteData(fid,size(md.forcing_accumulation,2),'Integer','forcing_accumulation_num_time_steps');
    98         WriteData(fid,md.forcing_accumulation(end,:),'Mat','forcing_accumulation_time_steps');
    99         for i=1:size(md.forcing_accumulation,2),
    100                 WriteData(fid,md.forcing_accumulation(1:end-1,i),'Mat',['forcing_accumulation_' num2str(i)]);
     96%deal with forcings
     97if ~isnans(md.forcings),
     98        forcingnames=fieldnames(md.forcings);
     99        numforcings=length(forcingnames);
     100        forcingtypes=zeros(numforcings,1);
     101        for i=1:numforcings,
     102                forcingtypes(i)=StringToEnum(forcingnames{i});
     103        end
     104        WriteData(fid,numforcings,'Integer','numforcings');
     105        WriteData(fid,forcingtypes,'Mat','forcingtypes');
     106
     107        for i=1:numforcings,
     108                forcing=md.forcings.(forcingnames{i});
     109                forcingname=forcingnames{i};
     110
     111                WriteData(fid,size(forcing,2),'Integer',['forcing_' forcingname '_num_time_steps']);
     112                WriteData(fid,forcing(end,:),'Mat',['forcing_' forcingname '_time_steps']);
     113                for j=1:size(forcing,2),
     114                        WriteData(fid,forcing(1:end-1,j),'Mat',['forcing_' forcingname '_' num2str(j)]);
     115                end
    101116        end
    102117else
    103         WriteData(fid,0,'Integer','forcing_accumulation_num_time_steps');
     118        WriteData(fid,0,'Integer','numforcings');
    104119end
    105120

  • issm/trunk/src/m/utils/BC/SetIceSheetBC.m

    r8320 r8365  
    3131if isnan(md.accumulation_rate),
    3232        md.accumulation_rate=zeros(md.numberofnodes,1);
    33         md.forcing_accumulation=[zeros(md.numberofnodes+1,1)];
     33        md.forcings.AccumulationRate=[zeros(md.numberofnodes+1,1)];
    3434        disp('      no accumulation_rate specified: values set as zero');
    3535end

  • issm/trunk/src/m/utils/BC/SetIceShelfBC.m

    r8320 r8365  
    6262if isnan(md.accumulation_rate),
    6363        md.accumulation_rate=zeros(md.numberofnodes,1);
    64         md.forcing_accumulation=zeros(md.numberofnodes+1,1);
     64        md.forcings.AccumulationRate=zeros(md.numberofnodes+1,1);
    6565        disp('      no accumulation_rate specified: values set as zero');
    6666end

  • issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m

    r8320 r8365  
    7272if isnan(md.accumulation_rate),
    7373        md.accumulation_rate=zeros(md.numberofnodes,1);
    74         md.forcing_accumulation=zeros(md.numberofnodes+1,1);
     74        md.forcings.AccumulationRate=zeros(md.numberofnodes+1,1);
    7575        disp('      no accumulation_rate specified: values set as zero');
    7676end
Note: See TracChangeset for help on using the changeset viewer.