Changeset 8363


Ignore:
Timestamp:
05/19/11 18:15:37 (14 years ago)
Author:
Eric.Larour
Message:

trunk: implemented forcing for accumulation in 2D. Verification and validation needed on Nicole's side.

Location:
issm/trunk/src
Files:
3 added
33 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/Inputs.cpp

    r8224 r8363  
    429429}
    430430/*}}}*/
     431/*FUNCTION Inputs::Configure{{{1*/
     432void Inputs::Configure(Parameters* parameters){
     433
     434        vector<Object*>::iterator object;
     435        Input* input=NULL;
     436
     437        for ( object=objects.begin() ; object < objects.end(); object++ ){
     438
     439                input=(Input*)(*object);
     440                input->Configure(parameters);
     441
     442        }
     443
     444}
     445/*}}}*/

  • issm/trunk/src/c/Container/Inputs.h

    r5749 r8363  
    4848                void GetParameterValue(int* pvalue,int enum_type);
    4949                void GetParameterValue(double* pvalue,int enum_type);
     50
     51                void Configure(Parameters* parameters);
    5052                /*}}}*/
    5153

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r8319 r8363  
    135135        PentaVertexInputEnum,
    136136        TriaVertexInputEnum,
     137        TriaVertexForcingEnum,
    137138        ControlInputEnum,
    138139        /*Params: */

  • issm/trunk/src/c/Makefile.am

    r8331 r8363  
    194194                                        ./objects/Inputs/TriaVertexInput.h\
    195195                                        ./objects/Inputs/TriaVertexInput.cpp\
     196                                        ./objects/Inputs/TriaVertexForcing.h\
     197                                        ./objects/Inputs/TriaVertexForcing.cpp\
    196198                                        ./objects/Inputs/PentaVertexInput.h\
    197199                                        ./objects/Inputs/PentaVertexInput.cpp\
    … …  
    839841                                        ./objects/Inputs/TriaVertexInput.h\
    840842                                        ./objects/Inputs/TriaVertexInput.cpp\
     843                                        ./objects/Inputs/TriaVertexForcing.h\
     844                                        ./objects/Inputs/TriaVertexForcing.cpp\
    841845                                        ./objects/Inputs/PentaVertexInput.h\
    842846                                        ./objects/Inputs/PentaVertexInput.cpp\

  • issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r8319 r8363  
    114114                case PentaVertexInputEnum : return "PentaVertexInput";
    115115                case TriaVertexInputEnum : return "TriaVertexInput";
     116                case TriaVertexForcingEnum : return "TriaVertexForcing";
    116117                case ControlInputEnum : return "ControlInput";
    117118                case ParamEnum : return "Param";

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

    r8330 r8363  
    115115
    116116        /*Update Elements in case we are running a transient solution: */
    117         if(analysis_counter==nummodels && (solution_type==Transient2DSolutionEnum | solution_type==Transient3DSolutionEnum)){
     117        if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){
    118118                UpdateElementsTransient(elements,iomodel,iomodel_handle,analysis_counter,analysis_type);
    119119        }

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

    r8331 r8363  
    1616
    1717        /*Intermediary*/
    18         int      i;
     18        int      i,j;
    1919        int      counter;
    2020        Element *element = NULL;
     21        char     fetchstring[100];
     22        double   time;
     23
     24       
     25        /*download whatever data will persist through forcing processing: */
     26        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    2127
    2228        /*Ok, let's go through forcings: */
     29        /*first, accumulation: */
     30        if(iomodel->forcing_accumulation_num_time_steps==0)_error_("accumulation forcing not available!");
     31
     32        iomodel->forcing_numtimesteps=iomodel->forcing_accumulation_num_time_steps;
     33        IoModelFetchData(&iomodel->timesteps,NULL,NULL,iomodel_handle,"forcing_accumulation_time_steps");
     34
     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);
    2339
    2440
     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                        }
     53                }
     54
     55
     56
     57                /*Free ressources:*/
     58                xfree((void**)&iomodel->forcing);
     59        }
     60               
     61
     62        /*Free ressources:*/
     63        xfree((void**)&iomodel->forcing);
     64        xfree((void**)&iomodel->timesteps);
     65        xfree((void**)&iomodel->elements);
     66
    2567}

  • issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r8319 r8363  
    112112        else if (strcmp(name,"PentaVertexInput")==0) return PentaVertexInputEnum;
    113113        else if (strcmp(name,"TriaVertexInput")==0) return TriaVertexInputEnum;
     114        else if (strcmp(name,"TriaVertexForcing")==0) return TriaVertexForcingEnum;
    114115        else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
    115116        else if (strcmp(name,"Param")==0) return ParamEnum;

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

    r7391 r8363  
    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;
    6162                virtual void   InputToResult(int enum_type,int step,double time)=0;
    6263                virtual void   ControlInputGetGradient(Vec gradient,int enum_type)=0;

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

    r8287 r8363  
    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!");}
    131132                int    UpdateShelfStatus(Vec new_shelf_nodes);
    132133                void   UpdateShelfFlags(double* new_shelf_nodes);

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

    r8287 r8363  
    27252725        /*point parameters to real dataset: */
    27262726        this->parameters=parametersin;
     2727
     2728        /*get inputs configured too: */
     2729        this->inputs->Configure(parameters);
     2730
    27272731
    27282732}
    … …  
    58635867}
    58645868/*}}}*/
     5869/*FUNCTION Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){{{1*/
     5870void Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){
     5871               
     5872        /*Intermediaries*/
     5873        int i,j;
     5874        int tria_vertex_ids[3];
     5875        double  nodeinputs[3];
     5876
     5877        /*Recover vertices ids needed to initialize inputs*/
     5878        for(i=0;i<3;i++) tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab
     5879
     5880        /*recover forcing for this time step: */
     5881        for(i=0;i<3;i++)nodeinputs[i]=iomodel->forcing[tria_vertex_ids[i]-1];
     5882
     5883        if(step==0){
     5884                /*This is the first time we are trying to replace the accumulation tria vertex input by an accumulation triav vertex
     5885                 *forcing, which is essentially a tria vertex input that can vary with time. So firt, if there is already
     5886                 an input with enum FieldEnum, squash it: */
     5887                this->inputs->AddInput(new TriaVertexForcing(FieldEnum,nodeinputs,time,iomodel->forcing_numtimesteps,this->parameters));
     5888        }
     5889        else{
     5890
     5891                /*Ok, we can't write over the existing forcing that was created in previous step 0. Find the input, and add
     5892                 *values for this time step: */
     5893                TriaVertexForcing* forcing=(TriaVertexForcing*)inputs->GetInput(FieldEnum); _assert_(forcing);
     5894
     5895                forcing->AddTimeValues(&nodeinputs[0],step,time);
     5896
     5897        }
     5898
     5899
     5900}
     5901/*}}}*/
    58655902/*FUNCTION Tria::UpdateShelfStatus{{{1*/
    58665903int  Tria::UpdateShelfStatus(Vec new_shelf_nodes){

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

    r8287 r8363  
    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);
    136137                int    UpdateShelfStatus(Vec new_shelf_nodes);
    137138                int    UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf);

  • issm/trunk/src/c/objects/Inputs/BoolInput.cpp

    r8224 r8363  
    237237}
    238238/*}}}*/
     239/*FUNCTION BoolInput::Configure{{{1*/
     240void BoolInput::Configure(Parameters* parameters){
     241        /*do nothing: */
     242}
     243/*}}}*/

  • issm/trunk/src/c/objects/Inputs/BoolInput.h

    r8129 r8363  
    4444                Input* PointwiseMax(Input* inputB){_error_("not implemented yet");};
    4545                ElementResult* SpawnResult(int step, double time);
     46                void Configure(Parameters* parameters);
     47                void AddTimeValues(double* values,int step,double time){_error_("not supported yet");};
    4648                /*}}}*/
    4749                /*numerics: {{{1*/

  • issm/trunk/src/c/objects/Inputs/ControlInput.cpp

    r8224 r8363  
    448448        values->VerticallyIntegrate(thickness_input);
    449449}/*}}}*/
     450/*FUNCTION ControlInput::Configure{{{1*/
     451void ControlInput::Configure(Parameters* parameters){
     452        /*do nothing: */
     453}
     454/*}}}*/

  • issm/trunk/src/c/objects/Inputs/ControlInput.h

    r8129 r8363  
    4848                Input* PointwiseMax(Input* inputB){_error_("not implemented yet");};
    4949                ElementResult* SpawnResult(int step, double time);
     50                void AddTimeValues(double* values,int step,double time){_error_("not supported yet");};
     51                void Configure(Parameters* parameters);
    5052                /*}}}*/
    5153                /*numerics: {{{1*/

  • issm/trunk/src/c/objects/Inputs/DoubleInput.cpp

    r8224 r8363  
    368368}
    369369/*}}}*/
     370/*FUNCTION DoubleInput::Configure{{{1*/
     371void DoubleInput::Configure(Parameters* parameters){
     372        /*do nothing: */
     373}
     374/*}}}*/

  • issm/trunk/src/c/objects/Inputs/DoubleInput.h

    r8129 r8363  
    4343                Input* PointwiseMax(Input* inputB);
    4444                ElementResult* SpawnResult(int step, double time);
     45                void AddTimeValues(double* values,int step,double time){_error_("not supported yet");};
     46                void Configure(Parameters* parameters);
    4547                /*}}}*/
    4648                /*numerics: {{{1*/

  • issm/trunk/src/c/objects/Inputs/Input.h

    r8129 r8363  
    1313class ElementResult;
    1414class GaussTria;
     15class Parameters;
    1516/*}}}*/
    1617
    … …  
    3839                virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0;
    3940                virtual void ChangeEnum(int newenumtype)=0;
     41                virtual void Configure(Parameters* parameters)=0;
    4042
    4143                virtual void   SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0;
    … …  
    6062                virtual Input* PointwiseMin(Input* inputmin)=0;
    6163                virtual ElementResult* SpawnResult(int step, double time)=0;
     64                virtual void AddTimeValues(double* values,int step,double time)=0;
    6265
    6366                /*}}}*/

  • issm/trunk/src/c/objects/Inputs/IntInput.cpp

    r8224 r8363  
    242242}
    243243/*}}}*/
     244/*FUNCTION IntInput::Configure{{{1*/
     245void IntInput::Configure(Parameters* parameters){
     246        /*do nothing: */
     247}
     248/*}}}*/

  • issm/trunk/src/c/objects/Inputs/IntInput.h

    r8129 r8363  
    4444                Input* PointwiseMax(Input* inputB){_error_("not implemented yet");};
    4545                ElementResult* SpawnResult(int step, double time);
     46                void AddTimeValues(double* values,int step,double time){_error_("not supported yet");};
     47                void Configure(Parameters* parameters);
    4648                /*}}}*/
    4749                /*numerics: {{{1*/

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp

    r8303 r8363  
    673673}
    674674/*}}}*/
     675/*FUNCTION PentaVertexInput::Configure{{{1*/
     676void PentaVertexInput::Configure(Parameters* parameters){
     677        /*do nothing: */
     678}
     679/*}}}*/

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.h

    r8129 r8363  
    4444                Input* PointwiseMax(Input* inputB);
    4545                ElementResult* SpawnResult(int step, double time);
     46                void AddTimeValues(double* values,int step,double time){_error_("not supported yet");};
     47                void Configure(Parameters* parameters);
    4648                /*}}}*/
    4749                /*numerics: {{{1*/

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

    r8303 r8363  
    465465}
    466466/*}}}*/
     467/*FUNCTION TriaVertexInput::Configure{{{1*/
     468void TriaVertexInput::Configure(Parameters* parameters){
     469        /*do nothing: */
     470}
     471/*}}}*/

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.h

    r8129 r8363  
    4444                Input* PointwiseMax(Input* inputB);
    4545                ElementResult* SpawnResult(int step, double time);
     46                void AddTimeValues(double* values,int step,double time){_error_("not supported yet");};
     47                void Configure(Parameters* parameters);
    4648                /*}}}*/
    4749                /*numerics: {{{1*/

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

    r8330 r8363  
    8989        xfree((void**)&this->melting_rate_correction);
    9090        xfree((void**)&this->accumulation_rate);
    91         xfree((void**)&this->forcing_accumulation);
     91        xfree((void**)&this->forcing);
     92        xfree((void**)&this->timesteps);
    9293        xfree((void**)&this->dhdt);
    9394        xfree((void**)&this->rheology_B);
    … …  
    220221        IoModelFetchData(&this->gl_melting_rate,iomodel_handle,"gl_melting_rate");
    221222        IoModelFetchData(&this->rheology_law,iomodel_handle,"rheology_law");
     223        IoModelFetchData(&this->forcing_accumulation_num_time_steps,iomodel_handle,"forcing_accumulation_num_time_steps");
    222224       
    223225        /*qmu: */
    … …  
    397399        /*!basal: */
    398400        this->accumulation_rate=NULL;
    399         this->forcing_accumulation=NULL;
    400401        this->dhdt=NULL;
     402
     403        /*forcings: */
     404        this->forcing_accumulation_num_time_steps=0;
     405        this->forcing_numtimesteps=0;
     406        this->forcing=NULL;
     407        this->timesteps=NULL;
    401408       
    402409        /*parameter output: */

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

    r8330 r8363  
    195195                int      melting_rate_correction_apply;
    196196                double*  accumulation_rate;
    197                 double*  forcing_accumulation;
    198197                double*  dhdt;
     198
     199                /*forcings: */
     200                int      forcing_accumulation_num_time_steps;
     201                int      forcing_numtimesteps;
     202                double*  forcing;
     203                double*  timesteps;
    199204
    200205                /*parameter output: */

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

    r8206 r8363  
    8080#include "./Inputs/PentaVertexInput.h"
    8181#include "./Inputs/TriaVertexInput.h"
     82#include "./Inputs/TriaVertexForcing.h"
    8283#include "./Inputs/ControlInput.h"
    8384

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

    r8319 r8363  
    5656                }
    5757                time+=dt;
     58                femmodel->parameters->SetParam(time,TimeEnum);
    5859                step+=1;
    5960
    … …  
    8586                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
    8687                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
     88                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AccumulationRateEnum,step,time);
    8789                        if(gl_migration!=NoneEnum)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time);
    8890

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

    r8319 r8363  
    5353                step+=1;
    5454                time+=dt;
     55                femmodel->parameters->SetParam(time,TimeEnum);
    5556
    5657                _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts,"    iteration number: ",step,"/",floor((finaltime-time)/dt)," dt [yr]: ",dt/yts);

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

    r8319 r8363  
    9393WriteData(fid,md.melting_rate,'Mat','melting_rate');
    9494WriteData(fid,md.melting_rate_correction_apply,'Integer','melting_rate_correction_apply');
    95 WriteData(fid,md.forcing_accumulation,'Mat','forcing_accumulation');
     95
     96if ~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)]);
     101        end
     102else
     103        WriteData(fid,0,'Integer','forcing_accumulation_num_time_steps');
     104end
     105
    96106if md.melting_rate_correction_apply,
    97107        WriteData(fid,md.melting_rate_correction,'Mat','melting_rate_correction');

  • issm/trunk/src/m/solutions/transient2d_core.m

    r8319 r8363  
    3333                step=step+1;
    3434                time=time+dt;
     35                femmodel.parameters.Time=time;
    3536
    3637                issmprintf(VerboseSolution,'\n%s%g%s%i%s%g\n','time [yr]: ',time,' iteration number: ',step,'/',floor(ndt/dt));

  • issm/trunk/src/m/solutions/transient3d_core.m

    r8319 r8363  
    3232                step=step+1;
    3333                time=time+dt;
     34                femmodel.parameters.Time=time;
    3435
    3536                issmprintf(VerboseSolution,'\n%s%g%s%i%s%g\n','time [yr] ',time,' iteration number: ',step,'/',floor(ndt/dt));
Note: See TracChangeset for help on using the changeset viewer.