Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 7847) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 7848) @@ -431,6 +431,11 @@ OptionDoubleEnum, OptionLogicalEnum, - OptionStructEnum - /*}}}*/ + OptionStructEnum, + /*}}}*/ + /*Rheology law (move too Material) {{{1*/ + RheologyLawEnum, + PatersonEnum, + ArrheniusEnum + /*}}}*/ }; Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp (revision 7847) +++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp (revision 7848) @@ -381,4 +381,7 @@ case OptionLogicalEnum : return "OptionLogical"; case OptionStructEnum : return "OptionStruct"; + case RheologyLawEnum : return "RheologyLaw"; + case PatersonEnum : return "Paterson"; + case ArrheniusEnum : return "Arrhenius"; default : return "unknown"; Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp (revision 7847) +++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp (revision 7848) @@ -379,4 +379,7 @@ else if (strcmp(name,"OptionLogical")==0) return OptionLogicalEnum; else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum; + else if (strcmp(name,"RheologyLaw")==0) return RheologyLawEnum; + else if (strcmp(name,"Paterson")==0) return PatersonEnum; + else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; else _error_("Enum %s not found",name); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 7847) +++ /issm/trunk/src/c/Makefile.am (revision 7848) @@ -327,4 +327,5 @@ ./shared/Elements/elements.h\ ./shared/Elements/Paterson.cpp\ + ./shared/Elements/Arrhenius.cpp\ ./shared/Elements/GetVerticesCoordinates.cpp\ ./shared/Elements/GetLocalDofList.cpp\ @@ -962,4 +963,5 @@ ./shared/Elements/elements.h\ ./shared/Elements/Paterson.cpp\ + ./shared/Elements/Arrhenius.cpp\ ./shared/Elements/GetVerticesCoordinates.cpp\ ./shared/Elements/GetLocalDofList.cpp\ Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 7847) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 7848) @@ -36,9 +36,9 @@ parameters = new Parameters(ParametersEnum); - parameters->AddObject(new IntParam(SolutionTypeEnum,solution_type)); - parameters->AddObject(new IntParam(AnalysisTypeEnum,analysis_type)); - parameters->AddObject(new IntParam(AnalysisCounterEnum,analysis_counter)); - if (iomodel->dim==2) parameters->AddObject(new IntParam(DimEnum,2)); - else parameters->AddObject(new IntParam(DimEnum,3)); + parameters->AddObject(new IntParam(SolutionTypeEnum,solution_type)); + parameters->AddObject(new IntParam(AnalysisTypeEnum,analysis_type)); + parameters->AddObject(new IntParam(AnalysisCounterEnum,analysis_counter)); + + parameters->AddObject(new IntParam(DimEnum,iomodel->dim)); parameters->AddObject(new BoolParam(IsHutterEnum,iomodel->ishutter)); parameters->AddObject(new BoolParam(IsMacAyealPattynEnum,iomodel->ismacayealpattyn)); @@ -79,4 +79,5 @@ parameters->AddObject(new BoolParam(IoGatherEnum,iomodel->io_gather)); parameters->AddObject(new IntParam(GroundingLineMigrationEnum,iomodel->gl_migration)); + parameters->AddObject(new IntParam(RheologyLawEnum,iomodel->rheology_law)); /*Deal with more complex parameters*/ Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 7847) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 7848) @@ -5095,8 +5095,9 @@ bool converged; - int i; + int i,rheology_law; + double xyz_list[NUMVERTICES][3]; double values[numdof]; double B[numdof]; - double B_average; + double B_average,s_average; int* doflist=NULL; @@ -5114,13 +5115,36 @@ } + /*Get all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); + this->inputs->GetParameterValue(&converged,ConvergedEnum); if(converged){ this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,values)); - /*Update Rheology only if convreged (we must make sure that the temperature is below melting point + /*Update Rheology only if converged (we must make sure that the temperature is below melting point * otherwise the rheology could be negative*/ - B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0); - for(i=0;imatice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); + this->parameters->FindParam(&rheology_law,RheologyLawEnum); + switch(rheology_law){ + case NoneEnum: + /*Do nothing: B is not temperature dependent*/ + break; + case PatersonEnum: + B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0); + for(i=0;imatice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); + break; + case ArrheniusEnum: + surface_input->GetParameterAverage(&s_average); + B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0, + s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0), + matice->GetN()); + for(i=0;imelting_rate_correction_apply,iomodel_handle,"melting_rate_correction_apply"); IoModelFetchData(&this->gl_melting_rate,iomodel_handle,"gl_melting_rate"); + IoModelFetchData(&this->rheology_law,iomodel_handle,"rheology_law"); /*qmu: */ @@ -225,5 +226,4 @@ IoModelFetchData(&this->qmu_save_femmodel,iomodel_handle,"qmu_save_femmodel"); } - /*i/o: */ @@ -321,4 +321,5 @@ this->rheology_n=NULL; this->rheology_B=NULL; + this->rheology_law=0; /*!solution parameters: */ Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 7847) +++ /issm/trunk/src/c/objects/IoModel.h (revision 7848) @@ -113,4 +113,5 @@ double* rheology_B; double* rheology_n; + int rheology_law; /*numerical parameters: */ Index: /issm/trunk/src/c/shared/Elements/Arrhenius.cpp =================================================================== --- /issm/trunk/src/c/shared/Elements/Arrhenius.cpp (revision 7848) +++ /issm/trunk/src/c/shared/Elements/Arrhenius.cpp (revision 7848) @@ -0,0 +1,50 @@ +/* \file Arrhenius.cpp + * \brief figure out B of ice for a certain temperature + */ + +#include "./elements.h" +#include + +double Arrhenius(double temperature,double depth,double n){ + /*Use EISMINT Parameterization for the rheology: Payne2000 + * + * A(T*) = A0 exp(-Q/RT*) + * + * A0 constant of proportionality + * = 3.61 * 10^-13 if T*<263.15K + * = 1.73 * 10^3 if T*>263.15K + * Q Activation energy for creep + * = 6.0 * 10^4 if T*<263.15K + * = 13.9 * 10^4 if T*>263.15K + * R Universal gas constant + * = 8.314 + * T* Absolute temperature corrected for the dependence of Tpmp on P + * = T - beta (s-z) + * + * Convert A to B : B = A^(-1/n) */ + + /*Some physical constants (Payne2000)*/ + double beta=8.66*pow(10.,-4.); + double R=8.314; + + /*Intermediaries*/ + double A,B,Tstar; + + /*convert temperature to absolute temperature*/ + _assert_(depth); + Tstar=temperature-beta*depth; + _assert_(Tstar); + + /*Get A*/ + if(Tstar<263.15){ + A=3.61*pow(10.,-13.) * exp( -6.*pow(10.,4.)/(R*Tstar)); + } + else{ + A=1.73*pow(10., 3.) * exp(-13.9*pow(10.,4.)/(R*Tstar)); + } + + /*Convert to B*/ + B=pow(A,-1./n); + + return B; +} Index: /issm/trunk/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk/src/c/shared/Elements/elements.h (revision 7847) +++ /issm/trunk/src/c/shared/Elements/elements.h (revision 7848) @@ -12,4 +12,5 @@ double Paterson(double temperature); +double Arrhenius(double temperature,double depth,double n); int GetVerticesCoordinates(double* xyz, Node** nodes, int numvertices); int GetNumberOfDofs( Node** nodes,int numnodes,int setenum,int approximation_enum); Index: /issm/trunk/src/m/classes/model.m =================================================================== --- /issm/trunk/src/m/classes/model.m (revision 7847) +++ /issm/trunk/src/m/classes/model.m (revision 7848) @@ -139,4 +139,5 @@ rheology_B=NaN; rheology_n=NaN; + rheology_law=0; %Geometrical parameters @@ -732,4 +733,7 @@ md.hydro_kn=0; + %Rheology law: what is the temperature dependence of B with T + %available: None, Paterson and Arrhenius + md.rheology_law=PatersonEnum; %i/o: Index: /issm/trunk/src/m/enum/ArrheniusEnum.m =================================================================== --- /issm/trunk/src/m/enum/ArrheniusEnum.m (revision 7848) +++ /issm/trunk/src/m/enum/ArrheniusEnum.m (revision 7848) @@ -0,0 +1,11 @@ +function macro=ArrheniusEnum() +%ARRHENIUSENUM - Enum of Arrhenius +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ArrheniusEnum() + +macro=367; Index: /issm/trunk/src/m/enum/EnumToString.m =================================================================== --- /issm/trunk/src/m/enum/EnumToString.m (revision 7847) +++ /issm/trunk/src/m/enum/EnumToString.m (revision 7848) @@ -376,4 +376,7 @@ case OptionLogicalEnum(), string='OptionLogical'; return case OptionStructEnum(), string='OptionStruct'; return +case RheologyLawEnum(), string='RheologyLaw'; return +case PatersonEnum(), string='Paterson'; return +case ArrheniusEnum(), string='Arrhenius'; return otherwise, error(['Enum ' num2str(enum) ' not found']); Index: /issm/trunk/src/m/enum/PatersonEnum.m =================================================================== --- /issm/trunk/src/m/enum/PatersonEnum.m (revision 7848) +++ /issm/trunk/src/m/enum/PatersonEnum.m (revision 7848) @@ -0,0 +1,11 @@ +function macro=PatersonEnum() +%PATERSONENUM - Enum of Paterson +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=PatersonEnum() + +macro=366; Index: /issm/trunk/src/m/enum/RheologyLawEnum.m =================================================================== --- /issm/trunk/src/m/enum/RheologyLawEnum.m (revision 7848) +++ /issm/trunk/src/m/enum/RheologyLawEnum.m (revision 7848) @@ -0,0 +1,11 @@ +function macro=RheologyLawEnum() +%RHEOLOGYLAWENUM - Enum of RheologyLaw +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=RheologyLawEnum() + +macro=365; Index: /issm/trunk/src/m/enum/StringToEnum.m =================================================================== --- /issm/trunk/src/m/enum/StringToEnum.m (revision 7847) +++ /issm/trunk/src/m/enum/StringToEnum.m (revision 7848) @@ -374,4 +374,7 @@ elseif (strcmpi(name,'OptionLogical')), enum=OptionLogicalEnum(); return elseif (strcmpi(name,'OptionStruct')), enum=OptionStructEnum(); return + elseif (strcmpi(name,'RheologyLaw')), enum=RheologyLawEnum(); return + elseif (strcmpi(name,'Paterson')), enum=PatersonEnum(); return + elseif (strcmpi(name,'Arrhenius')), enum=ArrheniusEnum(); return else error(['Enum ' name ' not found']); Index: /issm/trunk/src/m/model/marshall.m =================================================================== --- /issm/trunk/src/m/model/marshall.m (revision 7847) +++ /issm/trunk/src/m/model/marshall.m (revision 7848) @@ -104,4 +104,5 @@ WriteData(fid,md.rheology_B,'Mat','rheology_B'); WriteData(fid,md.rheology_n,'Mat','rheology_n'); +WriteData(fid,md.rheology_law,'Integer','rheology_law'); %Control methods