Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 14733) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 14734) @@ -509,5 +509,6 @@ ./modules/ModelProcessorx/Gia/CreateLoadsGia.cpp\ ./modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp\ - ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h + ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h\ + ./classes/GiaDeflectionCoreArgs.h #}}} Index: /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h =================================================================== --- /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h (revision 14734) +++ /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h (revision 14734) @@ -0,0 +1,31 @@ +/*!\file: GiaDeflectionCoreArgs.h + * \brief place holder for arguments to the GiaDeflectionCoreArgs routine + */ + +#ifndef _GIADEFLECTIONCOREARGS_H_ +#define _GIADEFLECTIONCOREARGS_H_ + +struct GiaDeflectionCoreArgs{ + + /*inputs: */ + IssmDouble ri; //radial distance from center of disk to vertex i + IssmDouble re; //radius of disk + IssmDouble* hes; //loading history (in ice thickness) + IssmDouble* times; //loading history times + int numtimes; //loading history length + IssmDouble currenttime; + + /*gia material parameters: */ + IssmDouble lithosphere_shear_modulus; + IssmDouble lithosphere_density; + IssmDouble mantle_shear_modulus; + IssmDouble mantle_viscosity; + IssmDouble mantle_density; + IssmDouble lithosphere_thickness; + + /*ice properties: */ + IssmDouble rho_ice; + +}; + +#endif Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 14733) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 14734) @@ -28,4 +28,5 @@ #include "./Update.h" #include "./FemModel.h" +#include "./GiaDeflectionCoreArgs.h" #include "./OptArgs.h" #include "./OptPars.h" Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14733) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14734) @@ -3057,19 +3057,48 @@ IssmDouble currenttime; int numtimes; + Input* thickness_input=NULL; + + /*gia material parameters: */ + IssmDouble lithosphere_shear_modulus; + IssmDouble lithosphere_density; + IssmDouble mantle_shear_modulus; + IssmDouble mantle_viscosity; + IssmDouble mantle_density; + Input* lithosphere_thickness_input=NULL; + IssmDouble lithosphere_thickness; + + /*ice properties: */ + IssmDouble rho_ice; /*output: */ IssmDouble wi; + /*arguments to GiaDeflectionCorex: */ + GiaDeflectionCoreArgs arguments; + /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - + /*what time is it? :*/ this->parameters->FindParam(¤ttime,TimeEnum); + /*recover material parameters: */ + lithosphere_shear_modulus=matpar->GetLithosphereShearModulus(); + lithosphere_density=matpar->GetLithosphereDensity(); + mantle_shear_modulus=matpar->GetMantleShearModulus(); + mantle_viscosity=matpar->GetMantleViscosity(); + mantle_density=matpar->GetMantleDensity(); + rho_ice=matpar->GetRhoIce(); + /*pull thickness averages: */ - Input* thickness_input=inputs->GetInput(ThicknessEnum); + thickness_input=inputs->GetInput(ThicknessEnum); if (!thickness_input)_error_("thickness input needed to compute gia deflection!"); thickness_input->GetInputAllTimeAverages(&hes,×,&numtimes); - + + /*recover lithosphere thickness: */ + lithosphere_thickness_input=inputs->GetInput(GiaLithosphereThicknessEnum); + if (!lithosphere_thickness_input)_error_("lithosphere thickness input needed to compute gia deflection!"); + lithosphere_thickness_input->GetInputAverage(&lithosphere_thickness); + /*pull area of this Tria: */ area=this->GetArea(); @@ -3083,14 +3112,32 @@ y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + /*start loading GiaDeflectionCore arguments: */ + arguments.re=re; + arguments.hes=hes; + arguments.times=times; + arguments.numtimes=numtimes; + arguments.currenttime=currenttime; + arguments.lithosphere_shear_modulus=lithosphere_shear_modulus; + arguments.lithosphere_density=lithosphere_density; + arguments.mantle_shear_modulus=mantle_shear_modulus; + arguments.mantle_viscosity=mantle_viscosity; + arguments.mantle_density=mantle_density; + arguments.lithosphere_thickness=lithosphere_thickness; + arguments.rho_ice=rho_ice; + for(i=0;iSetValue(i,wi,ADD_VAL); + } Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp (revision 14733) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp (revision 14734) @@ -60,5 +60,14 @@ _error_("Hydrology model "<Constant(&this->lithosphere_shear_modulus,MaterialsLithosphereShearModulusEnum); + iomodel->Constant(&this->lithosphere_density,MaterialsLithosphereDensityEnum); + iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum); + iomodel->Constant(&this->mantle_viscosity,MaterialsMantleViscosityEnum); + iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum); + + /*Unit conversion: */ + this->UnitConversion(); this->inputs=NULL; /*not used here*/ @@ -383,2 +392,38 @@ } /*}}}*/ + +/*FUNCTION Matpar::GetLithosphereShearModulus {{{*/ +IssmDouble Matpar::GetLithosphereShearModulus(){ + return lithosphere_shear_modulus; +} +/*}}}*/ +/*FUNCTION Matpar::GetLithosphereDensity {{{*/ +IssmDouble Matpar::GetLithosphereDensity(){ + return lithosphere_density; +} +/*}}}*/ +/*FUNCTION Matpar::GetMantleDensity {{{*/ +IssmDouble Matpar::GetMantleDensity(){ + return mantle_density; +} +/*}}}*/ +/*FUNCTION Matpar::GetMantleShearModulus {{{*/ +IssmDouble Matpar::GetMantleShearModulus(){ + return mantle_shear_modulus; +} +/*}}}*/ +/*FUNCTION Matpar::GetMantleViscosity {{{*/ +IssmDouble Matpar::GetMantleViscosity(){ + return mantle_viscosity; +} +/*}}}*/ +/*FUNCTION Matpar::UnitConversion {{{*/ +void Matpar::UnitConversion(void){ + + /*convert units of fields that were allocated using the Ext unit system, into the Iu unit system: */ + ::UnitConversion(&this->lithosphere_density,1,ExtToIuEnum,MaterialsLithosphereDensityEnum); + ::UnitConversion(&this->mantle_density,1,ExtToIuEnum,MaterialsMantleDensityEnum); + +} +/*}}}*/ + Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h (revision 14733) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h (revision 14734) @@ -43,4 +43,13 @@ IssmDouble sediment_transmitivity; IssmDouble water_compressibility; + + /*gia: */ + IssmDouble lithosphere_shear_modulus; + IssmDouble lithosphere_density; + IssmDouble mantle_shear_modulus; + IssmDouble mantle_viscosity; + IssmDouble mantle_density; + + public: Matpar(); @@ -111,6 +120,12 @@ IssmDouble PureIceEnthalpy(IssmDouble pressure); IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); + IssmDouble GetLithosphereShearModulus(); + IssmDouble GetLithosphereDensity(); + IssmDouble GetMantleShearModulus(); + IssmDouble GetMantleViscosity(); + IssmDouble GetMantleDensity(); void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure); + void UnitConversion(void); /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 14733) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 14734) @@ -7,4 +7,5 @@ #include "./GiaDeflectionCorex.h" +#include "../../classes/classes.h" #include "../../shared/shared.h" #include "../../include/include.h" @@ -13,9 +14,46 @@ #include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h" -void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble ri, IssmDouble re, IssmDouble* hes, IssmDouble* times, IssmDouble currentime,int numtimes){ +void GiaDeflectionCorex( IssmDouble* pwi, GiaDeflectionCoreArgs* arguments){ /*output: */ IssmDouble wi=1; + /*inputs: {{{*/ + IssmDouble ri; //radial distance from center of disk to vertex i + IssmDouble re; //radius of disk + IssmDouble* hes; //loading history (in ice thickness) + IssmDouble* times; //loading history times + int numtimes; //loading history length + IssmDouble currenttime; + + /*gia material parameters: */ + IssmDouble lithosphere_shear_modulus; + IssmDouble lithosphere_density; + IssmDouble mantle_shear_modulus; + IssmDouble mantle_viscosity; + IssmDouble mantle_density; + IssmDouble lithosphere_thickness; + + /*ice properties: */ + IssmDouble rho_ice; + /*}}}*/ + /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/ + ri=arguments->ri; + re=arguments->re; + hes=arguments->hes; + times=arguments->times; + numtimes=arguments->numtimes; + currenttime=arguments->currenttime; + lithosphere_shear_modulus=arguments->lithosphere_shear_modulus; + lithosphere_density=arguments->lithosphere_density; + mantle_shear_modulus=arguments->mantle_shear_modulus; + mantle_viscosity=arguments->mantle_viscosity; + mantle_density=arguments->mantle_density; + lithosphere_thickness=arguments->lithosphere_thickness; + rho_ice=arguments->rho_ice; + /*}}}*/ + + + /*allocate output pointer: */ *pwi=wi; Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (revision 14733) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (revision 14734) @@ -6,9 +6,9 @@ #define _GIADEFLECTIONCOREX_H -#include "../../classes/objects/objects.h" +#include "../../classes/classes.h" #include "../../Container/Container.h" /* local prototypes: */ -void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble ri, IssmDouble re, IssmDouble* hes, IssmDouble* times, IssmDouble currentime,int numtimes); +void GiaDeflectionCorex( IssmDouble* pwi, GiaDeflectionCoreArgs* arguments); #endif /* _GIADEFLECTIONCOREX_H */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Gia/UpdateElementsGia.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Gia/UpdateElementsGia.cpp (revision 14733) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Gia/UpdateElementsGia.cpp (revision 14734) @@ -22,5 +22,5 @@ iomodel->Constant(&dim,MeshDimensionEnum); iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); - iomodel->FetchData(1,MeshElementsEnum); + iomodel->FetchData(2,MeshElementsEnum,GiaLithosphereThicknessEnum); /*Update elements: */ @@ -35,6 +35,7 @@ iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,GiaLithosphereThicknessEnum); /*Free data: */ - iomodel->DeleteData(1,MeshElementsEnum); + iomodel->DeleteData(2,MeshElementsEnum,GiaLithosphereThicknessEnum); } Index: /issm/trunk-jpl/src/m/classes/matdamageice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 14733) +++ /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 14734) @@ -21,4 +21,13 @@ rheology_Z = NaN; rheology_law = ''; + + %gia: + lithosphere_shear_modulus = 0.; + lithosphere_density = 0.; + mantle_shear_modulus = 0.; + mantle_viscosity = 0.; + mantle_density = 0.; + + end methods @@ -79,4 +88,12 @@ %available: none, paterson and arrhenius obj.rheology_law='Paterson'; + + %GIA: + obj.lithosphere_shear_modulus = 6.7*10^10; %(Pa) + obj.lithosphere_density = 3.32; %(g/cm^-3) + obj.mantle_shear_modulus = 1.45*10^11; %(Pa) + obj.mantle_viscosity = 10^21; %(Pa.s) + obj.mantle_density = 3.34; %(g/cm^-3) + end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ @@ -89,4 +106,10 @@ md = checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'}); + md = checkfield(md,'materials.lithosphere_shear_modulus','>',0); + md = checkfield(md,'materials.lithosphere_density','>',0); + md = checkfield(md,'materials.mantle_shear_modulus','>',0); + md = checkfield(md,'materials.mantle_viscosity','>',0); + md = checkfield(md,'materials.mantle_density','>',0); + end % }}} function disp(obj) % {{{ @@ -108,4 +131,10 @@ fielddisplay(obj,'rheology_Z','rheology multiplier'); fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius'''); + fielddisplay(obj,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'); + fielddisplay(obj,'lithosphere_density','Lithosphere density [g/cm^-3]'); + fielddisplay(obj,'mantle_shear_modulus','Mantle shear modulus [Pa]'); + fielddisplay(obj,'mantle_viscosity','Mantle viscosity [Pa.s]'); + fielddisplay(obj,'mantle_density','Mantle density [g/cm^-3]'); + end % }}} function marshall(obj,fid) % {{{ @@ -126,4 +155,11 @@ WriteData(fid,'object',obj,'class','materials','fieldname','rheology_Z','format','DoubleMat','mattype',1); WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer'); + + WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_shear_modulus','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_density','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mantle_shear_modulus','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mantle_viscosity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mantle_density','format','Double'); + end % }}} end Index: /issm/trunk-jpl/test/NightlyRun/test330.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test330.m (revision 14733) +++ /issm/trunk-jpl/test/NightlyRun/test330.m (revision 14734) @@ -4,6 +4,6 @@ md=parameterize(md,'../Par/SquareSheetConstrained.par'); md=setflowequation(md,'macayeal','all'); -md.cluster=generic('name',oshostname(),'np',3); -%md.verbose=verbose('convergence',true,'solution',true,'module',true,'solver',true,'mprocessor',true); +md.cluster=generic('name',oshostname(),'np',8); +md.verbose=verbose('1111111'); md=solve(md,GiaSolutionEnum());