Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 22003) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 22004) @@ -99,4 +99,5 @@ ./classes/Materials/Materials.cpp\ ./classes/Materials/Matice.cpp\ + ./classes/Materials/Matlitho.cpp\ ./classes/Materials/Matestar.cpp\ ./classes/Materials/Matpar.cpp\ @@ -471,4 +472,16 @@ endif #}}} +#Love sources {{{ +if LOVE +issm_sources += ./cores/love_core.cpp\ + ./analyses/LoveAnalysis.cpp\ + ./modules/FourierLoveCorex/FourierLoveCorex.cpp\ + ./modules/FourierLoveCorex/lnb_param.f90\ + ./modules/FourierLoveCorex/model.f90\ + ./modules/FourierLoveCorex/util.f90\ + ./modules/FourierLoveCorex/lovenb_sub.f90\ + ./modules/FourierLoveCorex/love_numbers.f90 +endif +#}}} #Esa sources {{{ if ESA Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 22004) @@ -1525,5 +1525,5 @@ /*Update Rheology only if converged (we must make sure that the temperature is below melting point * otherwise the rheology could be negative*/ - element->FindParam(&rheology_law,MaterialsRheologyLawEnum); + rheology_law=element->material->GetRheologyLaw(); element->GetInputListOnNodes(&surface[0],SurfaceEnum); switch(rheology_law){ Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 22004) @@ -111,5 +111,5 @@ #endif #ifdef _HAVE_LOVE_ - //case LoveAnalysisEnum : return new LoveAnalysis(); + case LoveAnalysisEnum : return new LoveAnalysis(); #endif #ifdef _HAVE_LEVELSET_ Index: /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp (revision 22004) +++ /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp (revision 22004) @@ -0,0 +1,57 @@ +#include "./LoveAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +void LoveAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void LoveAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void LoveAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +int LoveAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + _error_("not needed!"); +}/*}}}*/ +void LoveAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + +}/*}}}*/ +void LoveAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void LoveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ + +/*Finite Element Analysis*/ +void LoveAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_(" not needed!"); +}/*}}}*/ +ElementVector* LoveAnalysis::CreateDVector(Element* element){/*{{{*/ + _error_(" not needed!"); +}/*}}}*/ +ElementMatrix* LoveAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ + _error_(" not needed!"); +}/*}}}*/ +ElementMatrix* LoveAnalysis::CreateKMatrix(Element* element){/*{{{*/ + _error_(" not needed!"); +}/*}}}*/ +ElementVector* LoveAnalysis::CreatePVector(Element* element){/*{{{*/ + _error_("not supported"); +}/*}}}*/ +void LoveAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + _error_("not supported"); +}/*}}}*/ +void LoveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + _error_("not supported"); +}/*}}}*/ +void LoveAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void LoveAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void LoveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + _error_("not supported"); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/LoveAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/LoveAnalysis.h (revision 22004) +++ /issm/trunk-jpl/src/c/analyses/LoveAnalysis.h (revision 22004) @@ -0,0 +1,35 @@ +/*! \file LoveAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _LoveAnalysis_ +#define _LoveAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class LoveAnalysis: public Analysis{ + + public: + /*Model processing*/ + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); +}; +#endif Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22004) @@ -809,5 +809,5 @@ /*Update Rheology only if converged (we must make sure that the temperature is below melting point * otherwise the rheology could be negative*/ - element->FindParam(&rheology_law,MaterialsRheologyLawEnum); + rheology_law=element->material->GetRheologyLaw(); element->GetInputListOnNodes(&surface[0],SurfaceEnum); switch(rheology_law){ Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 22003) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 22004) @@ -24,4 +24,5 @@ #include "./FreeSurfaceTopAnalysis.h" #include "./GiaIvinsAnalysis.h" +#include "./LoveAnalysis.h" #include "./EsaAnalysis.h" #include "./HydrologyDCEfficientAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 22004) @@ -46,5 +46,5 @@ /*retrieve material_id: */ - matpar_id = iomodel->numberofelements+1; + matpar_id = iomodel->matparcounter; /*retrieve material_id*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22004) @@ -678,5 +678,9 @@ analyses_temp[numanalyses++]=GiaIvinsAnalysisEnum; break; - + + case LoveSolutionEnum: + analyses_temp[numanalyses++]=LoveAnalysisEnum; + break; + case EsaSolutionEnum: analyses_temp[numanalyses++]=EsaAnalysisEnum; Index: /issm/trunk-jpl/src/c/classes/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.h (revision 22003) +++ /issm/trunk-jpl/src/c/classes/IoModel.h (revision 22004) @@ -87,4 +87,5 @@ int loadcounter; //keep track of how many loads are being created in each analysis int nodecounter; //keep track of how many nodes are being created in each analysis + int matparcounter; //keep track of where we end up putting this object accessed by every element. /*Methods*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 22004) @@ -42,4 +42,5 @@ virtual void GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; + virtual int GetRheologyLaw()=0; virtual bool IsDamage()=0; virtual bool IsEnhanced()=0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 22004) @@ -41,4 +41,7 @@ this->mid=matestar_mid; + /*rheology law:*/ + iomodel->FetchData(&this->rheology_law,"md.materials.rheology_law"); + /*Hooks: */ matestar_eid=index+1; @@ -367,4 +370,8 @@ void Matestar::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/ _error_("not implemented yet"); +} +/*}}}*/ +int Matestar::GetRheologyLaw(){/*{{{*/ + return this->rheology_law; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 22004) @@ -28,4 +28,5 @@ Hook *helement; Element *element; + int rheology_law; public: @@ -63,4 +64,5 @@ void GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff); void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + int GetRheologyLaw(); IssmDouble GetA(); IssmDouble GetAbar(); Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 22004) @@ -40,4 +40,5 @@ int materialtype; iomodel->FindConstant(&materialtype,"md.materials.type"); + iomodel->FindConstant(&this->rheology_law,"md.materials.rheology_law"); this->Init(matice_mid,index,materialtype); @@ -639,4 +640,8 @@ } /*}}}*/ +int Matice::GetRheologyLaw(){/*{{{*/ + return this->rheology_law; +} +/*}}}*/ void Matice::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/ /*Nothing updated yet*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 22004) @@ -28,4 +28,5 @@ bool isdamaged; bool isenhanced; + int rheology_law; Hook *helement; Element *element; @@ -67,4 +68,5 @@ void GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff); void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + int GetRheologyLaw(); IssmDouble GetA(); IssmDouble GetAbar(); Index: /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 22004) +++ /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 22004) @@ -0,0 +1,197 @@ +/*!\file Matlitho.c + * \brief: implementation of the Matlitho object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../classes.h" +#include "../../shared/shared.h" + +/*Matlitho constructors and destructor*/ +Matlitho::Matlitho(){/*{{{*/ + this->numlayers=0; + this->radius=NULL; + this->viscosity=NULL; + this->lame_lambda=NULL; + this->lame_mu=NULL; + this->burgers_viscosity=NULL; + this->burgers_mu=NULL; + this->density=NULL; + this->isburgers=NULL; + this->issolid=NULL; + return; +} +/*}}}*/ +Matlitho::Matlitho(int matlitho_mid, IoModel* iomodel){/*{{{*/ + + IssmDouble* isburgersd=NULL; + IssmDouble* issolidd=NULL; + + this->mid=matlitho_mid; + iomodel->FindConstant(&this->numlayers,"md.materials.numlayers"); + + this->radius=xNew(this->numlayers+1); + xMemCpy(this->radius, iomodel->Data("md.materials.radius"),this->numlayers+1); + + this->viscosity=xNew(this->numlayers); + xMemCpy(this->viscosity, iomodel->Data("md.materials.viscosity"),this->numlayers); + + this->lame_lambda=xNew(this->numlayers); + xMemCpy(this->lame_lambda, iomodel->Data("md.materials.lame_lambda"),this->numlayers); + + this->lame_mu=xNew(this->numlayers); + xMemCpy(this->lame_mu, iomodel->Data("md.materials.lame_mu"),this->numlayers); + + this->burgers_viscosity=xNew(this->numlayers); + xMemCpy(this->burgers_viscosity, iomodel->Data("md.materials.burgers_viscosity"),this->numlayers); + + this->burgers_mu=xNew(this->numlayers); + xMemCpy(this->burgers_mu, iomodel->Data("md.materials.burgers_mu"),this->numlayers); + + this->density=xNew(this->numlayers); + xMemCpy(this->density, iomodel->Data("md.materials.density"),this->numlayers); + + this->isburgers=xNew(this->numlayers); + xMemCpy(this->isburgers, iomodel->Data("md.materials.isburgers"),this->numlayers); + + this->issolid=xNew(this->numlayers); + xMemCpy(this->issolid, iomodel->Data("md.materials.issolid"),this->numlayers); + + /*isburgersd= xNew(this->numlayers); + this->isburgers=xNew(this->numlayers); + xMemCpy(isburgersd, iomodel->Data("md.materials.isburgers"),this->numlayers); + for (int i=0;inumlayers;i++)this->isburgers[i]=reCast(isburgersd[i]); + + issolidd= xNew(this->numlayers); + this->issolid=xNew(this->numlayers); + xMemCpy(issolidd, iomodel->Data("md.materials.issolid"),this->numlayers); + for (int i=0;inumlayers;i++)this->issolid[i]=reCast(issolidd[i]);*/ + + /*free ressources: */ + xDelete(isburgersd); + xDelete(issolidd); +} +/*}}}*/ +Matlitho::~Matlitho(){/*{{{*/ + + xDelete(radius); + xDelete(viscosity); + xDelete(lame_lambda); + xDelete(lame_mu); + xDelete(burgers_viscosity); + xDelete(burgers_mu); + xDelete(density); + xDelete(isburgers); + xDelete(issolid); + + return; +} +/*}}}*/ +void Matlitho::SetMid(int matlitho_mid){/*{{{*/ + this->mid=matlitho_mid; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +Object* Matlitho::copy() {/*{{{*/ + + /*Output*/ + Matlitho* matlitho; + + /*Initialize output*/ + matlitho=new Matlitho(*this); + + /*copy fields: */ + matlitho->mid=this->mid; + matlitho->numlayers=this->numlayers; + if(matlitho->numlayers){ + matlitho->radius=xNew(this->numlayers+1); xMemCpy(matlitho->radius, this->radius,this->numlayers+1); + matlitho->viscosity=xNew(this->numlayers); xMemCpy(matlitho->viscosity, this->viscosity,this->numlayers); + matlitho->lame_lambda=xNew(this->numlayers); xMemCpy(matlitho->lame_lambda, this->lame_lambda,this->numlayers); + matlitho->lame_mu=xNew(this->numlayers); xMemCpy(matlitho->lame_mu, this->lame_mu,this->numlayers); + matlitho->burgers_viscosity=xNew(this->numlayers); xMemCpy(matlitho->burgers_viscosity, this->burgers_viscosity,this->numlayers); + matlitho->burgers_mu=xNew(this->numlayers); xMemCpy(matlitho->burgers_mu, this->burgers_mu,this->numlayers); + matlitho->density=xNew(this->numlayers); xMemCpy(matlitho->density, this->density,this->numlayers); + matlitho->isburgers=xNew(this->numlayers); xMemCpy(matlitho->isburgers, this->isburgers,this->numlayers); + matlitho->issolid=xNew(this->numlayers); xMemCpy(matlitho->issolid, this->issolid,this->numlayers); + /*matlitho->isburgers=xNew(this->numlayers); for(int i=0;inumlayers;i++)matlitho->isburgers[i]=this->isburgers[i]; + matlitho->issolid=xNew(this->numlayers); for(int i=0;inumlayers;i++)matlitho->issolid[i]=this->issolid[i];*/ + } + + return matlitho; +} +/*}}}*/ +void Matlitho::DeepEcho(void){/*{{{*/ + + this->Echo(); +} +/*}}}*/ +void Matlitho::Echo(void){/*{{{*/ + + _printf_("Matlitho:\n"); + _printf_(" mid: " << mid << "\n"); + _printf_(" numlayers: " << numlayers << "\n"); + _printf_("layer radius viscosity lame_lambda lame_mu burgers_viscosity burgers_mu density isburgers issolid\n"); + for (int i=0;iFindConstant(&materials_type,"md.materials.type"); - this->mid = matpar_mid; + this->mid = iomodel->matparcounter; switch(materials_type){ @@ -169,4 +173,114 @@ /*slr:*/ iomodel->FindConstant(&this->earth_density,"md.materials.earth_density"); + + break; + case MaterialsEnum: + //we have several types of materials. Retrieve this info first: + iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature"); + + //go through list of materials, and create constant parameters accordingly: + for(int i=0;iFindConstant(&this->rho_ice,"md.materials.rho_ice"); + iomodel->FindConstant(&this->rho_water,"md.materials.rho_water"); + iomodel->FindConstant(&this->rho_freshwater,"md.materials.rho_freshwater"); + iomodel->FindConstant(&this->mu_water,"md.materials.mu_water"); + iomodel->FindConstant(&this->heatcapacity,"md.materials.heatcapacity"); + iomodel->FindConstant(&this->thermalconductivity,"md.materials.thermalconductivity"); + iomodel->FindConstant(&this->temperateiceconductivity,"md.materials.temperateiceconductivity"); + iomodel->FindConstant(&this->latentheat,"md.materials.latentheat"); + iomodel->FindConstant(&this->beta,"md.materials.beta"); + iomodel->FindConstant(&this->meltingpoint,"md.materials.meltingpoint"); + iomodel->FindConstant(&this->referencetemperature,"md.constants.referencetemperature"); + iomodel->FindConstant(&this->mixed_layer_capacity,"md.materials.mixed_layer_capacity"); + iomodel->FindConstant(&this->thermal_exchange_velocity,"md.materials.thermal_exchange_velocity"); + iomodel->FindConstant(&this->g,"md.constants.g"); + + switch(smb_model){ //{{{ + case SMBforcingEnum: + /*Nothing to add*/ + break; + case SMBgembEnum: + iomodel->FindConstant(&this->albedo_ice,"md.smb.aIce"); + iomodel->FindConstant(&this->albedo_snow,"md.smb.aSnow"); + break; + case SMBpddEnum: + iomodel->FindConstant(&this->desfac,"md.smb.desfac"); + iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); + iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm"); + break; + case SMBd18opddEnum: + iomodel->FindConstant(&this->desfac,"md.smb.desfac"); + iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); + iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm"); + iomodel->FindConstant(&this->dpermil,"md.smb.dpermil"); + case SMBgradientsEnum: + /*Nothing to add*/ + break; + case SMBgradientselaEnum: + /*Nothing to add*/ + break; + case SMBhenningEnum: + /*Nothing to add*/ + break; + case SMBcomponentsEnum: + /*Nothing to add*/ + break; + case SMBmeltcomponentsEnum: + /*Nothing to add*/ + break; + default: + _error_("Surface mass balance model "<FindConstant(&this->sediment_compressibility,"md.hydrology.sediment_compressibility"); + iomodel->FindConstant(&this->sediment_porosity,"md.hydrology.sediment_porosity"); + iomodel->FindConstant(&this->sediment_thickness,"md.hydrology.sediment_thickness"); + iomodel->FindConstant(&this->water_compressibility,"md.hydrology.water_compressibility"); + iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer"); + + if(isefficientlayer){ + iomodel->FindConstant(&this->epl_compressibility,"md.hydrology.epl_compressibility"); + iomodel->FindConstant(&this->epl_porosity,"md.hydrology.epl_porosity"); + iomodel->FindConstant(&this->epl_init_thickness,"md.hydrology.epl_initial_thickness"); + iomodel->FindConstant(&this->epl_colapse_thickness,"md.hydrology.epl_colapse_thickness"); + iomodel->FindConstant(&this->epl_max_thickness,"md.hydrology.epl_max_thickness"); + iomodel->FindConstant(&this->epl_conductivity,"md.hydrology.epl_conductivity"); + } + } + else if(hydrology_model==HydrologyshreveEnum){ + /*Nothing to add*/ + } + else if(hydrology_model==HydrologysommersEnum){ + /*Nothing to add*/ + } + else{ + _error_("Hydrology model "<FindConstant(&this->lithosphere_shear_modulus,"md.materials.lithosphere_shear_modulus"); + iomodel->FindConstant(&this->lithosphere_density,"md.materials.lithosphere_density"); + iomodel->FindConstant(&this->mantle_shear_modulus,"md.materials.mantle_shear_modulus"); + iomodel->FindConstant(&this->mantle_density,"md.materials.mantle_density"); + + /*slr:*/ + iomodel->FindConstant(&this->earth_density,"md.materials.earth_density"); + //}}} + break; + default: + _error_("Materials "<(nature); + break; break; Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 22004) @@ -75,5 +75,5 @@ public: Matpar(); - Matpar(int matpar_id, IoModel* iomodel); + Matpar(IoModel* iomodel); ~Matpar(); void SetMid(int matpar_mid); @@ -107,4 +107,5 @@ void GetViscosity_D(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; + int GetRheologyLaw(){_error_("not supported");}; IssmDouble GetA(){_error_("not supported");}; IssmDouble GetAbar(){_error_("not supported");}; Index: /issm/trunk-jpl/src/c/classes/Vertices.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Vertices.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/classes/Vertices.cpp (revision 22004) @@ -144,4 +144,6 @@ int vertex_max_sid; + if (this->Size()==0)return 0; + for(i=0;iSize();i++){ Vertex* vertex=xDynamicCast(this->GetObjectByOffset(i)); Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 22003) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 22004) @@ -84,4 +84,5 @@ #include "./Materials/Material.h" #include "./Materials/Matice.h" +#include "./Materials/Matlitho.h" #include "./Materials/Matestar.h" #include "./Materials/Matpar.h" Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 22004) @@ -75,4 +75,8 @@ solutioncore=&damage_core; break; + case LoveSolutionEnum: + solutioncore=&love_core; + break; + default: _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!"); Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 22003) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 22004) @@ -47,4 +47,5 @@ void dummy_core(FemModel* femmodel); void gia_core(FemModel* femmodel); +void love_core(FemModel* femmodel); void esa_core(FemModel* femmodel); void smb_core(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/cores/love_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 22004) +++ /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 22004) @@ -0,0 +1,91 @@ +/*!\file: love_core.cpp + * \brief: core of the LOVE numbers solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" +void love_core(FemModel* femmodel){ + + Vector *wg = NULL; + Vector *dwdtg = NULL; + IssmDouble *x = NULL; + IssmDouble *y = NULL; + + /*love parameters: */ + IssmDouble* frequencies=NULL; + int nfreq,dummy; + int sh_nmin,sh_nmax; + IssmDouble g0,r0,mu0; + bool allow_layer_deletion; + int forcing_type; + bool verbosemod =(int) VerboseModule(); + + /*parameters: */ + bool save_results; + + if(VerboseSolution()) _printf0_(" computing LOVE numbers\n"); + + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + /*recover love number parameters: */ + femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum); + femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy); + femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum); + femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum); + femmodel->parameters->FindParam(&g0,LoveG0Enum); + femmodel->parameters->FindParam(&r0,LoveR0Enum); + femmodel->parameters->FindParam(&mu0,LoveMu0Enum); + femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum); + femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum); + + /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/ + Matlitho* matlitho=NULL; + for (int i=femmodel->materials->Size()-1;i>=0;i--){ + Material* material=xDynamicCast(femmodel->materials->GetObjectByOffset(i)); + if (material->ObjectEnum()==MatlithoEnum)matlitho=xDynamicCast(material); + } + _assert_(matlitho); + + /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */ + IssmDouble* LoveKr = xNewZeroInit(nfreq*(sh_nmax+1)); + IssmDouble* LoveHr = xNewZeroInit(nfreq*(sh_nmax+1)); + IssmDouble* LoveLr = xNewZeroInit(nfreq*(sh_nmax+1)); + IssmDouble* LoveKi = xNewZeroInit(nfreq*(sh_nmax+1)); + IssmDouble* LoveHi = xNewZeroInit(nfreq*(sh_nmax+1)); + IssmDouble* LoveLi = xNewZeroInit(nfreq*(sh_nmax+1)); + + + + /*call the main module: */ + FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi, //output + nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs + matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu, + matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->isburgers, matlitho->issolid //matlitho inputs + ); + + /*Add love matrices to results:*/ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0)); + + /*Free ressources:*/ + xDelete(frequencies); + xDelete(LoveKr); + xDelete(LoveHr); + xDelete(LoveLr); + xDelete(LoveKi); + xDelete(LoveHi); + xDelete(LoveLi); + + + +} Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.cpp (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.cpp (revision 22004) @@ -0,0 +1,38 @@ +/*!\file FourierLoveCorex + * \brief: core solution module for love analysis. Based on Caron's fortran core. + */ + +#include "./FourierLoveCorex.h" + +#include "../../classes/classes.h" +#include "../../shared/shared.h" +#include "../../toolkits/toolkits.h" +#include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h" + +extern "C" { + int lnb_setup_( IssmDouble* LoveKr, IssmDouble* LoveKi, IssmDouble* LoveHr, IssmDouble* LoveHi, IssmDouble* LoveLr, + IssmDouble* LoveLi, int* pnfreq, IssmDouble* frequencies, int* psh_nmax, int* psh_nmin, + IssmDouble* pg0, IssmDouble* pr0, IssmDouble* pmu0,bool* pallow_layer_deletion, int* pforcing_type, bool* pverbosemod, + int* pnumlayers, IssmDouble* radius, IssmDouble* viscosity, IssmDouble* lame_lambda, IssmDouble* lame_mu, + IssmDouble* burgers_viscosity, IssmDouble* burgers_mu, IssmDouble* density, IssmDouble* isburgers, IssmDouble* issolid); + +} + +/*}}}*/ + +void FourierLoveCorex( IssmDouble* LoveKr, IssmDouble* LoveKi, IssmDouble* LoveHr, IssmDouble* LoveHi, IssmDouble* LoveLr, + IssmDouble* LoveLi, int nfreq, IssmDouble* frequencies, int sh_nmax, int sh_nmin, + IssmDouble g0, IssmDouble r0, IssmDouble mu0,bool allow_layer_deletion, int forcing_type, bool verbosemod, + int numlayers, IssmDouble* radius, IssmDouble* viscosity, IssmDouble* lame_lambda, IssmDouble* lame_mu, + IssmDouble* burgers_viscosity, IssmDouble* burgers_mu, IssmDouble* density, IssmDouble* isburgers, IssmDouble* issolid){ + + /*Call fortran driver: */ + lnb_setup_( LoveKr, LoveKi, LoveHr, LoveHi, LoveLr, + LoveLi, &nfreq, frequencies, &sh_nmax, &sh_nmin, + &g0, &r0, &mu0, &allow_layer_deletion, &forcing_type, &verbosemod, + &numlayers, radius, viscosity, lame_lambda, lame_mu, + burgers_viscosity, burgers_mu, density, isburgers, issolid); + + + +} Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.h (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.h (revision 22004) @@ -0,0 +1,17 @@ +/*!\file: FourierLoveCorex.h + * \brief header file for ... + */ + +#ifndef _FOURIERLOVECOREX_H +#define _FOURIERLOVECOREX_H + +#include "../../classes/classes.h" + +/* local prototypes: */ +void FourierLoveCorex( IssmDouble* LoveKr, IssmDouble* LoveKi, IssmDouble* LoveHr, IssmDouble* LoveHi, IssmDouble* LoveLr, + IssmDouble* LoveLi,int nfreq, IssmDouble* frequencies, int sh_nmax, int sh_nmin, IssmDouble g0, IssmDouble r0, + IssmDouble mu0,bool allow_layer_deletion, int forcing_type, bool verbosemod, int numlayers, IssmDouble* radius, + IssmDouble* viscosity,IssmDouble* lame_lambda, IssmDouble* lame_mu, IssmDouble* burgers_viscosity, IssmDouble* burgers_mu, + IssmDouble* density,IssmDouble* isburgers, IssmDouble* issolid); + +#endif /* _FOURIERLOVECOREX_H */ Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lnb_param.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lnb_param.f90 (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lnb_param.f90 (revision 22004) @@ -0,0 +1,35 @@ +module lnb_param + + implicit none + + public + + !Physical Constants + double precision, parameter :: GG = 6.67259e-11! gravitational constant + double precision :: ra,rb,rc ! radii of the surface, CMB and ICB + double precision :: pi,ro_mean,go_surf,go_icb,go_cmb !some constants + ! GG --> gravitational constant + ! ra, rb, rc --> Surface, CMB and ICB radius + ! ro_mean --> mean density of the planet + ! go_surf --> mean surface gravity + ! go_icb/cmb --> gravity at the Inner Core/Core-Mantle Boundary + double precision :: mu0,go0,r0 !3 normalization constants + +!Numerical Constants + double precision :: epsdb, layerrap !used in detecting underflows in the lovenb integration + double complex :: valini,freq ! integration constant and angular frequency being processed + logical :: allow_layer_del,display, benchmark_spada !do we allow deletion of layers in case of underflow, enable some prints in terminal + +!GIA problem setup + integer :: deg,ifmin,ifmax ! deg=SH degree, ifmin/ifmax=which kind of forcing is applied (usually we want loading) + logical :: incompressible !choose which mode of calculation is activated + +!Earth Model + integer :: nfreq, degmax, degmin + integer :: nlayer! number of layers + double precision, dimension(:), pointer :: roc,lac,muc,vic,vic2,muc2 ! Earth layers parameters + double precision, dimension(:), pointer :: radius ! Earth layers parameters + double precision :: vicratio,mucratio !2 Burgers long-term to short-term ratios + logical, dimension(:), pointer :: soliddim,radbc,burgers ! is this layer solid (if not it's liquid), does it have specific boundary conditions and does it have a Burgers rheology (if not, it uses Maxwell) + +end module lnb_param Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 (revision 22004) @@ -0,0 +1,362 @@ + subroutine lnb_setup( ksr, ksi, hsr, hsi, lsr, lsi, dnfreq, frequencies, ddegmax, ddegmin, & + dgo0, dr0, dmu0, dallow_layer_del, source_type, ddisplay, dnlayer, dradius, & + dvic, dlac, dmuc, dvic2, dmuc2, droc, dburgers, dsoliddim) + use lnb_param + use util + implicit none + + integer :: source_type, fr, i, IERR + double precision :: T,cst, fluid_tscale,g + double precision, dimension(dnfreq,ddegmax+1) :: hsr,hsi,lsr,lsi,ksr,ksi + double complex, dimension(dnfreq,ddegmax+1) :: hs,ls,ks + double precision, dimension(dnfreq) :: frequencies + + double precision :: dmu0,dgo0,dr0 + logical :: dallow_layer_del, ddisplay + integer :: dnfreq, ddegmax, ddegmin + integer :: dnlayer + double precision, dimension(dnlayer) :: droc,dlac,dmuc,dvic,dvic2,dmuc2 + double precision, dimension(dnlayer+1) :: dradius + double precision, dimension(dnlayer) :: dsoliddim,dburgers + + !assign local variables to public variables + mu0=dmu0; go0=dgo0; r0=dr0; allow_layer_del=dallow_layer_del; + nfreq=dnfreq;degmax=ddegmax;degmin=ddegmin;nlayer=dnlayer; + allocate(roc(nlayer),lac(nlayer),muc(nlayer),vic(nlayer),vic2(nlayer),muc2(nlayer)) + roc=droc; lac=dlac; muc=dmuc; vic=dvic; vic2=dvic2; muc2=dmuc2; + allocate(radius(nlayer+1),soliddim(nlayer),burgers(nlayer),radbc(nlayer+1)) + radius=dradius; + + display=ddisplay + + do i=1,nlayer + if (NINT(dsoliddim(i))==1) then + soliddim(i)=.true. + else + soliddim(i)=.false. + endif + if (NINT(dburgers(i))==1) then + burgers(i)=.true. + else + burgers(i)=.false. + endif + enddo + + pi=dacos(-1.d0) + + ra=radius(nlayer+1) + rb = 3480.d3 + rc = 1221.5d3 + r0=ra + cst = 365.25d0*24.d0*3600.d0*1000.d0 + + benchmark_spada=.false. + !frequencies=1E-5/cst*2*pi + + hs(:,:)=complex(0.d0,0.d0) + ks=hs + ls=hs + + call love_numbers(frequencies,source_type,hs,ls,ks) + + hsr=dble(hs) + hsi=dimag(hs) + ksr=dble(ks) + ksi=dimag(ks) + lsr=dble(ls) + lsi=dimag(ls) + +end subroutine + +!==================== + subroutine love_numbers(frequencies,source_type,hs,ls,ks) +!==================== + + use model + use lovenb_sub + + double complex, dimension(:,:), pointer :: bc + double complex, dimension(:,:), pointer :: f + character*40, dimension(:), pointer :: sourcs + integer, dimension(:), pointer :: indx,indx2 + integer :: nbc,nfext,ifreq,ntheta2,nphi2 + double precision :: cst,prec, delta!,h1,h2,k1,k2,l1,l2,dh1,dh2,dk1,dk2,dl1,dl2 + double precision :: he,ke,le + double precision :: frequencies(nfreq), fluid_tscale + double complex :: loveh,lovel,lovek + logical :: logi + integer :: i,j,k,n, fr,m,l, source_type + double complex :: hs(nfreq,degmax+1), ls(nfreq,degmax+1), ks(nfreq,degmax+1) + + + double precision :: cpu_time1(100) + integer :: cpu_count + + pi=acos(-1.d0) + + + open(unit=101, file='lastrun_log') + cpu_count=0 + cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count)) + !write(101,*), 'Done !', cpu_time1(cpu_count), 's' + !write(101,*),'' + + + !double precision, dimension(:), pointer :: vech,vecl,veck,vecd,vec + + radbc(:)=.false. + radbc(1)=.true. + radbc(nlayer+1)=.true. + do i=2,nlayer + if ((.not.soliddim(i)).or.(.not.soliddim(i-1))) then + radbc(i)=.true. + end if + end do + + + ! Test number of interfaces and of boundary conditions + nbc = 0 + do i=1,nlayer+1 + if (radbc(i)) nbc=nbc+6 + enddo + + + ifmin=source_type + ifmax=source_type + + + if (display) then + print*,'model init' + endif + + + call model_init2(nbc) + + + if (display) then + print*, 'done' + endif + + nfext = 13 ! Number of potential excitation sources + allocate( bc(nbc,nbc), indx(nbc) ) + allocate( f(nfext,nbc), sourcs(nfext) ) + + sourcs(1)='ICB --Volumetric Potential' + sourcs(2)='ICB --Pressure' + sourcs(3)='ICB --Loading' + sourcs(4)='ICB --Tangential Traction' + sourcs(5)='CMB --Volumetric Potential' + sourcs(6)='CMB --Pressure' + sourcs(7)='CMB --Loading' + sourcs(8)='CMB --Tangential Traction' + sourcs(9)='SURF--Volumetric Potential' + sourcs(10)='SURF--Pressure' + sourcs(11)='SURF--Loading' + sourcs(12)='SURF--Tangential Traction' + + + if (display) print*, 'source_type =', source_type, ', ', sourcs(source_type) + if (source_type<9) then + print*,'Error: Internal loading not supported yet, please input source_type between 9 and 12' + print*,'Reference:' + do i=1,12 + print*,i,sourcs(i) + enddo + return + endif + + +! Print/write control + !display=.false. ! Printing on the terminal + + ! Calculation optimization (avoids underflows at high degree) + !allow_layer_del = .true. ! Do we allow deletion of the central layers in the calculation if necessary? + layerrap = 1.d0 ! Max absolute ratio between love numbers at the top and bottom of the last layer + + ! Mode search parameters + freq0 = -1.d10 + logi = .true. + cst = 365.25d0*24.d0*3600.d0*1000.d0 + prec = 1.d-8 + firstmode = 0.d0 + + + cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count)) + !write(101,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + !write(101,*),'' + if (display) then + write(*,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + write(*,*),'' + endif + +! Spherical Harmonics initialization + + + +! -- SH degree loop + + + bc(:,:)=0.d0 + + + !write(101,*), 'Calculating impulse response' + if (display) then + write(*,*), 'Calculating impulse response' + endif + if (maxval(frequencies).ne.0 .and. display) then + print*, 'Frequency ranging from ', minval(dabs(frequencies)), & + ' to ', maxval(frequencies),' s^-1' + endif + + if (frequencies(1)==0 .and. display) then + print*,'Elastic calculation' + endif + deg=degmin-1 + + + + + do l = degmin,degmax + deg=deg+1 + !write(*,*) + !write(*,*) + !write(101,fmt='(A19, 10X,A2,I3,A1,I3)',ADVANCE='NO'), repeat(char(8),19),'l=',deg,'/',degmax + if (display) then + write(*,fmt='(A23, 10X,A2,I5,A1,I5)',ADVANCE='NO'), repeat(char(8),23),'l=',deg,'/',degmax + endif + + !write(*,*) + !!write(101,*), fr, T + + + + + do fr=1,nfreq + + !!write(101,*),fr + + + !the way Fourier Transform algorithms works : they calculate for fmin=0 to fmax=1/dt + !by aliasing for this algo, f=1/(2*dt) to f=1/dt is the same as f=-1/(2*dt) to f=0 + !but physics-wise, the Earth response at high frequency is not what we are looking for + !the [0 ; 1/(2*dt)] interval and its negative counterpart are what we are looking for + !therefore in our love number calculation the frenquency set must be [0 ; 1/(2*dt) ] U [-1/(2*dt) 0[ + + !if (fr==1) then ! the variable freq is, in fact, the pulsation = 2*pi*frequency + ! fluid_tscale=-575.502 / (deg+91.1765) -0.176471 ! this formula sets fluid time scale = 1e6 kyr at degree 2, 1e5 kyr at degree 20 and 1e3 kyr at degree 90, which allows no error growth and still the fluid number as the benchmark at all degree (tested up to degmax=147) + !print*,fluid_tscale + ! freq=complex(0.d0,1.d0)*2.d0*pi*(10.d0**(fluid_tscale)/cst) ! empiric time period that changes according to harmonic degree to ensure getting fluid response without crazy error growing for high degree +! it is meant to approximate the theoretical freq=0.d0 that should be assessed for fft + !elseif (fr <= nfreq/2+1) then + ! freq=complex(0.d0,1.d0)*(dble(fr-1)/T)*2.d0*pi + !else + ! freq=complex(0.d0,1.d0)*(-dble(nfreq+1-fr)/T)*2.d0*pi + !end if + + freq=complex(0.d0,frequencies(fr)) + !print*,fr,freq + + !if (deg==degmin) write(31,*),fr,real(freq)*cst,aimag(freq)*cst,time(fr) + !freq=dble(fr)/T + !-- Elastic Love number calculation + !freq=2.d0*pi/(1e4*cst)*complex(0.d0,1.d0) + + !if (deg==2) then ! for rotationnal feedback + !ifmin=9;ifmax=9 ! Sets tidal calculation + !bc(:,:)=0.d0 + !call boundary_conditions_matrix(bc,indx,nbc) + !call external_forcing(deg,f,sourcs,nfext,nbc) + !call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) + !k2tidal(fr)=lovek-1.d0 + !h2tidal(fr)=loveh + !l2tidal(fr)=lovel + !ifmin=11;ifmax=11 ! Sets back loading calculation + !endif + !print*,(deg) + bc(:,:)=0.d0 + call boundary_conditions_matrix(bc,indx,nbc) + !print*,'bc ok', layerrap + call external_forcing(deg,f,sourcs,nfext,nbc) + !print*,'forcing ok', layerrap + call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) + + !print*,'first_sol', layerrap, epsdb + + ! Automatic reduction of the number of layers when the attenuation with depth becomes too strong + if (allow_layer_del.eqv..true.) then + + do while ((layerrap<=epsdb).and.(nbc>12).or.(isnan(layerrap)) ) + ! print*,'trying to delete layer', layerrap, epsdb, nbc + !write(101,*) + !write(101,*) 'Rapport Nombre de Love surface/profondeur faible : ', layerrap + !write(101,*) ' Changement d''interface de debut d''integration' + + if (display) then + write(*,*) + write(*,*) 'Surface/Depth Love number ratio small: ', layerrap + write(*,*) ' Changing the interface where the integration starts' + endif + + nbc = nbc-6 + n = 1 + do while (.not.radbc(n)) + n = n+1 + end do + radbc(n) =.false. + !write(101,*) ' New start interface: ', radius(n+1)/1.d3,' km' + if (display) then + write(*,*) ' New start interface: ', radius(n+1)/1.d3,' km' + endif + + deallocate( bc, indx, f, sourcs ) + allocate( bc(nbc,nbc), indx(nbc) ) + allocate( f(nfext,nbc), sourcs(nfext) ) + + + bc(:,:)=0.d0 + call boundary_conditions_matrix(bc,indx,nbc) + call external_forcing(deg,f,sourcs,nfext,nbc) + call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) + + end do + end if + + !-- Saving Love numbers + + + hs(fr,deg+1) = loveh + ks(fr,deg+1) = lovek - 1.d0 + ls(fr,deg+1) = lovel + + !print*,dble(hs(fr,deg+1)), dimag(hs(fr,deg+1)),dble(ks(fr,deg+1)), dimag(ks(fr,deg+1)) + + !if (fr==1) then ! if this is supposed to be fluid response + ! ! then cut off the imaginary part, which is inherited from approximating t=infinity to a few million years (see above the setting of freq) + ! hs(fr,deg+1)=real(hs(fr,deg+1)) + ! ks(fr,deg+1)=real(ks(fr,deg+1)) + ! ls(fr,deg+1)=real(ls(fr,deg+1)) + !endif +!~~~~~~~~~~~~~~TEST ZONE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + if (benchmark_spada) then + write(17,*), deg,fr,real(hs(fr,deg+1)), aimag(hs(fr,deg+1)),real(ks(fr,deg+1)), aimag(ks(fr,deg+1)) + endif +!~~~~~~~~~~~~~~END TEST ZONE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + !write(*,*) 'Elastic Love Numbers: ',ke,he,le + !write(ifile1,*) deg,ke,he,le + end do + + end do + + + cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count)) + !write(101,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + !write(101,*),'' + if (display) then + write(*,*), 'Love number calc done :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + write(*,*),'' + endif + + deallocate( bc, indx, f, sourcs ) +end subroutine Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90.bak =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90.bak (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90.bak (revision 22004) @@ -0,0 +1,362 @@ + program lnb_setup + + use lnb_param + use util + + integer :: source_type, fr, i + double precision :: T,cst, fluid_tscale,g + double complex, allocatable :: hs(:,:), ls(:,:), ks(:,:) + double precision, allocatable :: frequencies(:) + + degmax=180 + degmin=1 + nfreq=500 + nlayer=6 + allocate( radbc(nlayer+1), soliddim(nlayer) ) + allocate( roc(nlayer), lac(nlayer), muc(nlayer), muc2(nlayer), vic(nlayer), vic2(nlayer)) + allocate(radius(nlayer+1), burgers(nlayer), frequencies(nfreq)) + allocate( hs(nfreq,degmax+1), ls(nfreq,degmax+1), ks(nfreq,degmax+1) ) + + pi=dacos(-1.d0) + + allow_layer_del=.true. + + + radius(1) = 100.d3; radius(2) = 1221.5d3; radius(3) = 3480.d3; + radius(4) = 5701.d3; radius(5) = 5951.d3; radius(6) = 6271.d3; radius(7) = 6371.d3 + vic(1) = 0.d0; vic(2) = 0.d0; vic(3) = 2.d22; vic(4) = 5.d20; vic(5) = 5.d20; vic(6) = 1.d40 + incompressible=.false. + !ra=6378.137d0 + muc(1) = 1.d6; muc(2) = 0.d0; muc(3) = 2.20425d11 + muc(4) = 0.870257d11; muc(5) = 0.870257d11; muc(6) = 0.510268d11 + roc(1) = 10987.9d0; roc(2) = 10986.9d0; roc(3) = 4903.58d0 + roc(4) = 3628.29d0; roc(5) = 3628.29d0; roc(6) = 3113.94d0 + roc(6)=3300.d0 + + + !roc(6)=3300.d0 + muc2=muc;vic2=vic;lac=1e11 + !burgers(:)=.false. + soliddim(:) = .true. + soliddim(2) = .false. + ra=radius(nlayer+1) + !ra = 6378.137d3 + rb = 3480.d3 + rc = 1221.5d3 + + go0=10 + mu0=1e11 + r0=ra + + radbc(:)=.false. + radbc(1)=.true. + radbc(nlayer+1)=.true. + do i=2,nlayer + if ((.not.soliddim(i)).or.(.not.soliddim(i-1))) then + radbc(i)=.true. + end if + end do + + source_type=11 + cst = 365.25d0*24.d0*3600.d0*1000.d0 + T=8*cst + burgers(1:6)=.false. + print*,'ok',T, burgers(1:6), ' ', soliddim(1:nlayer), radbc(1:nlayer),ra, rb, rc + + + do i=1,nlayer + print*, radius(i+1), roc(i), muc(i), lac(i), vic(i) + print*, '' + enddo + benchmark_spada=.false. + !frequencies=1E-5/cst*2*pi + + do fr=1,nfreq + if (fr==1) then ! the variable freq is, in fact, the pulsation = 2*pi*frequency + fluid_tscale=-575.502 / (deg+91.1765) -0.176471 ! this formula sets fluid time scale = 1e6 kyr at degree 2, 1e5 kyr at degree 20 and 1e3 kyr at degree 90, which allows no error growth and still the fluid number as the benchmark at all degree (tested up to degmax=147) + !print*,fluid_tscale + frequencies(fr)=2.d0*pi*(10.d0**(fluid_tscale)/cst) +! it is meant to approximate the theoretical freq=0.d0 that should be assessed for fft + elseif (fr <= nfreq/2+1) then + frequencies(fr)=(dble(fr-1)/T)*2.d0*pi + else + frequencies(fr)=(-dble(nfreq+1-fr)/T)*2.d0*pi + end if + enddo + + display=.true. + call love_numbers(frequencies,source_type,hs,ls,ks) + !print*, hs + do i=1,degmax + do fr=1,nfreq + write(17,*), real(hs(fr,i+1)),aimag(hs(fr,i+1)), real(ls(fr,i+1)),aimag(ls(fr,i+1)), real(ks(fr,i+1)),aimag(ks(fr,i+1)) + enddo + enddo + + +end program + +!==================== + subroutine love_numbers(frequencies,source_type,hs,ls,ks) +!==================== + + use lnb_param + use model + use lovenb_sub + + double complex, dimension(:,:), pointer :: bc + double complex, dimension(:,:), pointer :: f + character*40, dimension(:), pointer :: sourcs + integer, dimension(:), pointer :: indx,indx2 + integer :: nbc,nfext,ifreq,ntheta2,nphi2 + double precision :: cst,prec, delta!,h1,h2,k1,k2,l1,l2,dh1,dh2,dk1,dk2,dl1,dl2 + double precision :: he,ke,le + double precision :: frequencies(nfreq), fluid_tscale + double complex :: loveh,lovel,lovek + logical :: logi + integer :: i,j,k,n, fr,m,l, source_type + double complex :: hs(nfreq,degmax+1), ls(nfreq,degmax+1), ks(nfreq,degmax+1) + + + double precision :: cpu_time1(100) + integer :: cpu_count + + + pi=acos(-1.d0) + + + open(unit=101, file='lastrun_log') + cpu_count=0 + cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count)) + !write(101,*), 'Done !', cpu_time1(cpu_count), 's' + !write(101,*),'' + + + !double precision, dimension(:), pointer :: vech,vecl,veck,vecd,vec + + + + ! Test number of interfaces and of boundary conditions + nbc = 0 + do i=1,nlayer+1 + if (radbc(i)) nbc=nbc+6 + enddo + + + ifmin=source_type + ifmax=source_type + + + print*,'model init' + call model_init2(nbc,nlayer) + print*, 'done' + + nfext = 13 ! Number of potential excitation sources + allocate( bc(nbc,nbc), indx(nbc) ) + allocate( f(nfext,nbc), sourcs(nfext) ) + + sourcs(1)='ICB --Volumetric Potential' + sourcs(2)='ICB --Pressure' + sourcs(3)='ICB --Loading' + sourcs(4)='ICB --Tangential Traction' + sourcs(5)='CMB --Volumetric Potential' + sourcs(6)='CMB --Pressure' + sourcs(7)='CMB --Loading' + sourcs(8)='CMB --Tangential Traction' + sourcs(9)='SURF--Volumetric Potential' + sourcs(10)='SURF--Pressure' + sourcs(11)='SURF--Loading' + sourcs(12)='SURF--Tangential Traction' + + + print*, 'source_type =', source_type, ', ', sourcs(source_type) + if (source_type<9) then + print*,'Error: Internal loading not supported yet, please input source_type between 9 and 12' + print*,'Reference:' + do i=1,12 + print*,i,sourcs(i) + enddo + return + endif + +! Print/write control + !display=.false. ! Printing on the terminal + + ! Calculation optimization (avoids underflows at high degree) + !allow_layer_del = .true. ! Do we allow deletion of the central layers in the calculation if necessary? + layerrap = 1.d0 ! Max absolute ratio between love numbers at the top and bottom of the last layer + + ! Mode search parameters + freq0 = -1.d10 + logi = .true. + cst = 365.25d0*24.d0*3600.d0*1000.d0 + prec = 1.d-8 + firstmode = 0.d0 + + + cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count)) + !write(101,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + !write(101,*),'' + if (display) then + write(*,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + write(*,*),'' + endif + +! Spherical Harmonics initialization + + + +! -- SH degree loop + + + bc(:,:)=0.d0 + + + !write(101,*), 'Calculating impulse response' + if (display) then + write(*,*), 'Calculating impulse response' + endif + + deg=degmin-1 + do l = degmin,degmax + deg=deg+1 + !write(*,*) + !write(*,*) + !write(101,fmt='(A19, 10X,A2,I3,A1,I3)',ADVANCE='NO'), repeat(char(8),19),'l=',deg,'/',degmax + if (display) then + write(*,fmt='(A23, 10X,A2,I5,A1,I5)',ADVANCE='NO'), repeat(char(8),23),'l=',deg,'/',degmax + endif + + !write(*,*) + !!write(101,*), fr, T + + + + + + + do fr=1,nfreq + + !!write(101,*),fr + + + !the way Fourier Transform algorithms works : they calculate for fmin=0 to fmax=1/dt + !by aliasing for this algo, f=1/(2*dt) to f=1/dt is the same as f=-1/(2*dt) to f=0 + !but physics-wise, the Earth response at high frequency is not what we are looking for + !the [0 ; 1/(2*dt)] interval and its negative counterpart are what we are looking for + !therefore in our love number calculation the frenquency set must be [0 ; 1/(2*dt) ] U [-1/(2*dt) 0[ + + !if (fr==1) then ! the variable freq is, in fact, the pulsation = 2*pi*frequency + ! fluid_tscale=-575.502 / (deg+91.1765) -0.176471 ! this formula sets fluid time scale = 1e6 kyr at degree 2, 1e5 kyr at degree 20 and 1e3 kyr at degree 90, which allows no error growth and still the fluid number as the benchmark at all degree (tested up to degmax=147) + !print*,fluid_tscale + ! freq=complex(0.d0,1.d0)*2.d0*pi*(10.d0**(fluid_tscale)/cst) ! empiric time period that changes according to harmonic degree to ensure getting fluid response without crazy error growing for high degree +! it is meant to approximate the theoretical freq=0.d0 that should be assessed for fft + !elseif (fr <= nfreq/2+1) then + ! freq=complex(0.d0,1.d0)*(dble(fr-1)/T)*2.d0*pi + !else + ! freq=complex(0.d0,1.d0)*(-dble(nfreq+1-fr)/T)*2.d0*pi + !end if + + freq=complex(0.d0,frequencies(fr)) + !print*,freq + + !if (deg==degmin) write(31,*),fr,real(freq)*cst,aimag(freq)*cst,time(fr) + !freq=dble(fr)/T + !-- Elastic Love number calculation + !freq=2.d0*pi/(1e4*cst)*complex(0.d0,1.d0) + + !if (deg==2) then ! for rotationnal feedback + !ifmin=9;ifmax=9 ! Sets tidal calculation + !bc(:,:)=0.d0 + !call boundary_conditions_matrix(bc,indx,nbc) + !call external_forcing(deg,f,sourcs,nfext,nbc) + !call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) + !k2tidal(fr)=lovek-1.d0 + !h2tidal(fr)=loveh + !l2tidal(fr)=lovel + !ifmin=11;ifmax=11 ! Sets back loading calculation + !endif + !print*,(deg) + bc(:,:)=0.d0 + call boundary_conditions_matrix(bc,indx,nbc) + call external_forcing(deg,f,sourcs,nfext,nbc) + call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) + + ! Automatic reduction of the number of layers when the attenuation with depth becomes too strong + if (allow_layer_del.eqv..true.) then + + do while ((layerrap<=epsdb).and.(nbc>12).or.(isnan(layerrap)) ) + + !write(101,*) + !write(101,*) 'Rapport Nombre de Love surface/profondeur faible : ', layerrap + !write(101,*) ' Changement d''interface de debut d''integration' + + if (display) then + write(*,*) + write(*,*) 'Surface/Depth Love number ratio small: ', layerrap + write(*,*) ' Changing the interface where the integration starts' + endif + + nbc = nbc-6 + n = 1 + do while (.not.radbc(n)) + n = n+1 + end do + radbc(n) =.false. + !write(101,*) ' New start interface: ', radius(n+1)/1.d3,' km' + if (display) then + write(*,*) ' New start interface: ', radius(n+1)/1.d3,' km' + endif + + deallocate( bc, indx, f, sourcs ) + allocate( bc(nbc,nbc), indx(nbc) ) + allocate( f(nfext,nbc), sourcs(nfext) ) + + + bc(:,:)=0.d0 + call boundary_conditions_matrix(bc,indx,nbc) + call external_forcing(deg,f,sourcs,nfext,nbc) + call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) + + end do + end if + + !-- Saving Love numbers + + + hs(fr,deg+1) = loveh + ks(fr,deg+1) = lovek - 1.d0 + ls(fr,deg+1) = lovel + + !if (fr==1) then ! if this is supposed to be fluid response + ! ! then cut off the imaginary part, which is inherited from approximating t=infinity to a few million years (see above the setting of freq) + ! hs(fr,deg+1)=real(hs(fr,deg+1)) + ! ks(fr,deg+1)=real(ks(fr,deg+1)) + ! ls(fr,deg+1)=real(ls(fr,deg+1)) + !endif +!~~~~~~~~~~~~~~TEST ZONE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + if (benchmark_spada) then + write(17,*), deg,fr,real(hs(fr,deg+1)), aimag(hs(fr,deg+1)),real(ks(fr,deg+1)), aimag(ks(fr,deg+1)) + endif +!~~~~~~~~~~~~~~END TEST ZONE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + + !write(*,*) 'Elastic Love Numbers: ',ke,he,le + !write(ifile1,*) deg,ke,he,le + end do + + + + + + end do + + cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count)) + !write(101,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + !write(101,*),'' + if (display) then + write(*,*), 'Love number calc done :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's' + write(*,*),'' + endif + +end subroutine Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 (revision 22004) @@ -0,0 +1,609 @@ +module lovenb_sub + + use util + use lnb_param + + implicit none + + contains + + subroutine boundary_conditions_matrix(bc,indx,nbc) +!------------------------------------- + + implicit none + double complex :: la, mu + double precision :: xmin,xmax,eps,hini,hmin + double precision :: x,ro1,ro2,g,d,one + double complex, dimension (:,:), pointer :: bc + integer, dimension(:), pointer :: indx + double complex, dimension(6) :: ystart,dydx + double complex, dimension(6,6) :: ystep + integer :: ny,i,j,k,ii,jj,kk,nok,nbad,is,nbc,ibc,ifirst,nstep + logical :: solid,vali + + !-- Initialisation of boundary condition matrix + + indx(:) = 0 + + !-- Integration + + eps= 1.d-10 + hini=1.d-10 + hmin=0.d0 + + ! Underflows risk tests in the numerical integration + ! The starting ystart are chosen to be different than 1 to limit under and overflows in the Runge-Kutta + if (-300.d0 zabs(lovel1s/lovel)) layerrap = zabs(lovel1s/lovel) + if (layerrap > zabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0))) & + layerrap = zabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0)) + endif + + loveh1s = loveh1 + lovel1s = lovel1 + lovek1s = lovek1 + rads = radius(j) + !if (radius(j)==rb) then + ! write(*,*) 'CMB ',n,loveh1,lovel1,lovek1-(rb/ra)**deg + !else if (radius(j)==rc) then + ! write(*,*) 'ICB ',n,loveh1,lovel1,lovek1-(rc/ra)**deg + !else + ! if (j<10) then + ! write(*,*) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg + ! else + ! write(*,*) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg + ! end if + !endif + end if + end do + + end do + end if + + + !301 format(a4,5x,'n=',i3,5x,'h=',f14.10,5x,'l=',f14.10,5x,'k=',f14.10) + 301 format(a4,5x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) + 305 format(a4,i1,4x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) + 306 format(a4,i2,3x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) + 302 format(i5,5x,f14.10,5x,f14.10,5x,f14.10) + 303 format(i5,5x,f14.10,5x,f14.10,5x,f14.10,5x,f14.10) + 304 format(i5,5x,d20.10,5x,d20.10,5x,d20.10) + + end subroutine solution + + + subroutine solution_nbloutput(n,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) +!------------------------------------------------- + + integer :: n,nfext,nbc,i,j,ic,ibc,ii,jj,icmin,ifirst + double precision, dimension(:,:), pointer :: bc,f + integer, dimension(:), pointer :: indx + character*40, dimension(:), pointer :: sourcs + double precision, dimension(nbc) :: sc + double precision :: loveh,lovel,lovek,loveh1,lovel1,lovek1,loveh1s,lovel1s,lovek1s + double precision :: sumh,suml,sumk,d,val,delta,sumy2,sumy4,sumy6,rads + double precision, dimension(6) :: icbh,icbl,icbk,y2srf,y4srf,y6srf + double precision, dimension(nbc,nbc) :: bcsav + logical :: ok + + + !do ic=1,6 +! if (nbc==18) then +! icbh(ic) = bc(1,ic) +! icbl(ic) = bc(3,ic) +! icbk(ic) = bc(5,ic) +!! elseif (nbc==24) then +! icbh(ic) = bc(7,ic+3) +! icbl(ic) = bc(9,ic+3) +! icbk(ic) = bc(11,ic+3) +! end if +! y2srf(ic) = bc(nbc-4,ic+nbc-9) +! y4srf(ic) = bc(nbc-2,ic+nbc-9) +! y6srf(ic) = bc(nbc,ic+nbc-9) + ! end do + + bcsav(:,:) = bc(:,:) + rads=0 + ifirst = 0 + i = 0 + do while (ifirst==0) + i = i+1 + if (radbc(i)) ifirst = i + end do + !do i = 1,nbc +! do j = 1,nbc +! if (bc(i,j)/=0.d0) print*, i,j,bc(i,j) +! enddo + !enddo + !write(*,*) + + !-- LU decomposition and determinant + ok = .true. + + call ludcmp(bc,nbc,nbc,indx,d) + do i=1,nbc + d=d*bc(i,i) + end do + write(*,*) + !write(*,*),dlog(freq),' Matrix det = ',d/dabs(d)*dlog(dabs(d)) + +!-- Resolution + + if (ok) then + do i = ifmin,ifmax + + do j=1,nbc + sc(j)=f(i,j) + end do + + write(*,*) 'Source: ',sourcs(i) + call lubksb(bc,nbc,nbc,indx,sc) + + loveh = sc(nbc-2)*ra*go_surf + lovel = sc(nbc-1)*ra*go_surf + lovek = sc(nbc)*ra*go0 + delta = 1.d0-dble(n+1)/dble(n)*(lovek-1.d0)+2.d0/dble(n)*loveh + + sumy2=0.d0 + sumy4=0.d0 + sumy6=0.d0 + + do ic=1,6 + sumy2 = sumy2+bcsav(nbc-4,ic+nbc-9)*sc(ic+nbc-9) + sumy4 = sumy4+bcsav(nbc-2,ic+nbc-9)*sc(ic+nbc-9) + sumy6 = sumy6+bcsav(nbc,ic+nbc-9)*sc(ic+nbc-9) + end do + + !write(*,*) 'SURF',n,loveh,lovel,lovek-1.d0 + + ibc = 0 + do j = nlayer,1,-1 + if (radbc(j)) then + ibc=ibc+1 + if (soliddim(j)) then + loveh1 = sc(nbc - ibc*6 -3 +1)*ra*go_surf*valini + lovel1 = sc(nbc -ibc*6 -3 +3)*ra*go_surf*valini + lovek1 = sc(nbc -ibc*6 -3 +5)*ra*go0*valini + else + sumh=0.d0 + suml=0.d0 + sumk=0.d0 + ii = nbc - (ibc+1)*6 + jj = nbc - (ibc+1)*6 -3 + if (j==ifirst) then + icmin = 4 + else + icmin = 1 + end if + do ic=icmin,6 + !print*, ii+1,jj+ic + sumh = sumh+bcsav(ii+1,jj+ic)*sc(jj+ic) + suml = suml+bcsav(ii+3,jj+ic)*sc(jj+ic) + sumk = sumk+bcsav(ii+5,jj+ic)*sc(jj+ic) + end do + loveh1 = sumh*ra*go_surf*valini + lovel1 = suml*ra*go_surf*valini + lovek1 = sumk*ra*go0*valini + end if + if (j==ifirst) then + layerrap = dabs(loveh1s/loveh) + if (layerrap > dabs(lovel1s/lovel)) layerrap = dabs(lovel1s/lovel) + if (layerrap > dabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0))) & + layerrap = dabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0)) + endif + + loveh1s = loveh1 + lovel1s = lovel1 + lovek1s = lovek1 + rads = radius(j) + !if (radius(j)==rb) then + ! write(*,301) 'CMB ',n,loveh1,lovel1,lovek1-(rb/ra)**deg + !else if (radius(j)==rc) then + ! write(*,301) 'ICB ',n,loveh1,lovel1,lovek1-(rc/ra)**deg + !else + ! if (j<10) then + ! write(*,305) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg + ! else + ! write(*,306) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg + ! end if + !endif + end if + end do + + end do + end if + + !write(*,*),dlog10(aimag(freq)*365.d0*24.d0*3600.d0*1000.d0),' Matrix det = ',d, lovek + + + + !301 format(a4,5x,'n=',i3,5x,'h=',f14.10,5x,'l=',f14.10,5x,'k=',f14.10) + 301 format(a4,5x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) + 305 format(a4,i1,4x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) + 306 format(a4,i2,3x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) + 302 format(i5,5x,f14.10,5x,f14.10,5x,f14.10) + 303 format(i5,5x,e14.10,5x,e14.10,5x,e14.10,5x,e14.10) + 304 format(i5,5x,d20.10,5x,d20.10,5x,d20.10) + + end subroutine solution_nbloutput + + + +end module lovenb_sub Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 (revision 22004) @@ -0,0 +1,422 @@ +module model + + implicit none + + contains + + subroutine model_init2(nbc) +!--------------------- + use lnb_param + + double precision :: ro,g,drad,la,mu,ro1,ro2,x + logical :: solid,test + integer :: nbc,i,n,j + + + + !- Calculus parameters + pi=dacos(-1.d0) + epsdb=epsilon(1.d0) + + !write(101,*), 'Epsilon machine = ',epsdb + !write(101,*), 'Constants:' + !write(101,*), ' pi = ',pi + !write(101,*), ' G = ',GG + !write(101,*), 'Dimensioned scales:' + !write(101,*), ' rigidity mu0 = ',mu0 + !write(101,*), ' gravity go0 = ',go0 + !write(101,*), ' distance r0 = ',r0 + + + + + !call earth_homogeneous_layers_prembased(radius,nlayer,roc,lac,muc,g,solid) !yields the Earth layers radii, density, lamé parameters, gravity and solid/liquid state boolean identifier + + !if (benchmark_spada) then + ! muc(1) = 1.d6; muc(2) = 0.d0; muc(3) = 2.28340d11 + ! muc(4) = 1.05490d11; muc(5) = 0.70363d11; muc(6) = 0.50605d11 + ! roc(1) = 10751.d0; roc(2) = 10750.d0; roc(3) = 4978.d0 + ! roc(6)=3300.d0 + !else + !muc(1) = 1.d6; muc(2) = 0.d0; muc(3) = 2.20425d11 + !muc(4) = 0.870257d11; muc(5) = 0.870257d11; muc(6) = 0.510268d11 + !roc(1) = 10987.9d0; roc(2) = 10986.9d0; roc(3) = 4903.58d0 + !roc(4) = 3628.29d0; roc(5) = 3628.29d0; roc(6) = 3113.94d0 + !endif + + !vic2=vic/vicratio + !muc2=muc/mucratio + + if (display) then + print*, ' Layer top boundary rho la mu eta' + do i = 1,nlayer + ! print*, radius(i+1),roc(i),lac(i),muc(i),vic(i) + if (burgers(i)) then + print*, radius(i+1),roc(i),lac(i),muc(i),vic(i),muc2(i),vic2(i), 'B' + else + print*, radius(i+1),roc(i),lac(i),muc(i),vic(i) + endif + end do + + endif + !write(101,*), + + + call earth_nlayers(1.d0,ro,la,mu,g,solid) + go_surf = g + ro_mean = 3.d0/4.d0*go_surf/(pi*GG*ra) + if (display) then + print*, 'Surface Gravity: ', go_surf, ' m.s^-2, mean density: ', ro_mean, ' kg.m^-3' + print*, 'Center of the Earth approximated at r=', radius(1)/1e3, ' km' + endif + + ! Test number of interfaces and of boundary conditions + n = 0 + do i=1,nlayer+1 + if (radbc(i)) n=n+1 + enddo + if (nbc/=(6*n)) then + !write(101,*), + print*, 'Error nbc should be ',6*n,' instead of ',nbc + stop + endif + + ! Test solid-liquid boundaries + test = .false. + do i = 1,nlayer + if (radbc(i)) test = .true. + if (test) then + if ((.not.soliddim(i)).and.((.not.radbc(i)).or.(.not.radbc(i+1)))) then + print*, 'Error in model definition: one fluid layer is not set with correct boundary conditions' + stop + endif + endif + enddo + + !-- Display a few quantities + !write(101,*), ' surface radius = ',ra + !write(101,*), ' CMB radius = ',rb + !write(101,*), ' ICB radius = ',rc + !write(101,*), ' go surface = ',go_surf + !write(101,*), ' mean density = ',ro_mean + + !-- Diplsay ICB conditions + !x=rc/ra + !write(101,*), + !call earth_nlayers_viscoelas(x-epsdb,ro1,complex(la,0.d0),complex(mu,0d0),g,solid) + !write(101,*), 'ICB, g = ',g + !write(101,*), 'ICB, ro(Inner Core) = ', ro1 + !call earth_nlayers_viscoelas(x,ro2,complex(la,0.d0),complex(mu,0d0),g,solid) + !write(101,*), 'ICB, g = ',g + !write(101,*), 'ICB, ro(Outer Core) = ', ro2 + !go_icb = g + + !-- Display CMB conditions + !x=rb/ra + !write(101,*), + !call earth_nlayers_viscoelas(x-epsdb,ro1,complex(la,0.d0),complex(mu,0d0),g,solid) + !write(101,*), 'CMB, g = ',g + !write(101,*), 'CMB, ro(Outer Core) = ', ro1 + !call earth_nlayers_viscoelas(x,ro2,complex(la,0.d0),complex(mu,0d0),g,solid) + !write(101,*), 'CMB, g = ',g + !write(101,*), 'CMB, ro(Mantle) = ', ro2 + !go_cmb = g + + end subroutine model_init2 + + + subroutine earth_homogeneous_layers_prembased(radin,nl,rot,lat,mut,GR,solid) + use lnb_param + + double precision :: x0,GR + logical :: solid + integer :: i,j,nl + doubleprecision :: la,mu,ro,vp,vs,RT + doubleprecision :: t1,t2,t3,t4,r1,r2 + doubleprecision, dimension(14) :: r + doubleprecision, dimension(13,4) :: d,p,s + doubleprecision, dimension(:), pointer :: radin + doubleprecision, dimension(:), pointer :: rot,lat,mut + doubleprecision, dimension(nl+1) :: rad + + + RT = ra + r(1) = 0.d0; r(2) = 1221.5d0; r(3) = 3480.d0; r(4) = 3630.d0 + r(5) = 5600.d0; r(6) = 5701.d0; r(7) = 5771.d0; r(8) = 5971.d0 + r(9) = 6151.d0; r(10) = 6291.d0; r(11) = 6346.6d0 + r(12) = 6356.d0; r(13) = 6368.d0; r(14) = ra + +! if ((r(14)*1.d3) /= ra) stop'Problem in prem: Earth radius /= ra !' +! + d(:,:) = 0.d0 + d(1,1) = 13.0885d0; d(1,3) = -8.8381d0 + d(2,1) = 12.5815d0; d(2,2) = -1.2638d0; d(2,3) = -3.6426d0; d(2,4) = -5.5281d0 + d(3,1) = 7.9565d0 ; d(3,2) = -6.4761; d(3,3) = 5.5283d0; d(3,4) = -3.0807d0 + d(4,1) = 7.9565d0 ; d(4,2) = -6.4761; d(4,3) = 5.5283d0; d(4,4) = -3.0807d0 + d(5,1) = 7.9565d0 ; d(5,2) = -6.4761; d(5,3) = 5.5283d0; d(5,4) = -3.0807d0 + d(6,1) = 5.3197d0 ; d(6,2) = -1.4836d0 + d(7,1) = 11.2494d0; d(7,2) = -8.0298d0 + d(8,1) = 7.1089d0 ; d(8,2) = -3.8045d0 + d(9,1) = 2.6910d0 ; d(9,2) = 0.6924d0 + d(10,1) = 2.6910d0; d(10,2) = 0.6924d0 + d(11,1) = 2.9d0 + d(12,1) = 2.6d0 + +! ocean + if (soliddim(13).eqv..false.) then + d(13,1) = 1.02d0 + +! continental + else + d(13,1) = d(12,1) + end if + + p(:,:) = 0.d0 + p(1,1) = 11.2622d0 ; p(1,3) = -6.3640d0 + p(2,1) = 11.0487d0 ; p(2,2) = -4.0362d0; p(2,3) = 4.8023d0; p(2,4) = -13.5732d0 + p(3,1) = 15.3891d0 ; p(3,2) = -5.3181d0; p(3,3) = 5.5242d0; p(3,4) = -2.5514d0 + p(4,1) = 24.952d0 ; p(4,2) = -40.4673d0; p(4,3) = 51.4832d0; p(4,4) = -26.6419d0 + p(5,1) = 29.2766d0 ; p(5,2) = -23.6027d0; p(5,3) = 5.5242d0; p(5,4) = -2.5514d0 + p(6,1) = 19.0957d0 ; p(6,2) = -9.8672d0 + p(7,1) = 39.7027d0 ; p(7,2) = -32.6166d0 + p(8,1) = 20.3926d0 ; p(8,2) = -12.2569d0 + p(9,1) = 4.1875d0 ; p(9,2) = 3.9382d0 + p(10,1) = 4.1875d0 ; p(10,2) = 3.9382d0 + p(11,1) = 6.8d0 + p(12,1) = 5.8d0 +! +! ocean + if (soliddim(13).eqv..false.) then + p(13,1) = 1.45d0 +! +! continental + else + p(13,1) = p(12,1) + end if +!---- +! + s(:,:) = 0.d0 +! + s(1,1) = 3.6678d0; s(1,3) = -4.4475d0 + + s(3,1) = 6.9254d0; s(3,2) = 1.4672d0; s(3,3) = -2.0834d0; s(3,4) = 0.9783d0 + s(4,1) = 11.1671d0; s(4,2) = -13.7818d0; s(4,3) = 17.4575d0; s(4,4) = -9.2777d0 + s(5,1) = 22.3459d0; s(5,2) = -17.2473d0; s(5,3) = -2.0834d0; s(5,4) = 0.9783d0 + s(6,1) = 9.9839d0; s(6,2) = -4.9324 + s(7,1) = 22.3512d0; s(7,2) = -18.5856d0 + s(8,1) = 8.9496d0; s(8,2) = -4.4597 + s(9,1) = 2.1519d0; s(9,2) = 2.3481d0 + s(10,1) = 2.1519d0; s(10,2) = 2.3481d0 + s(11,1) = 3.9d0 + s(12,1) = 3.2d0 +! +! ocean (please don't modify) + if (soliddim(13).eqv..false.) then +! +! continental + else + s(13,1) = s(12,1) + end if +! +! + r(:) = r(:)*1.d3 + + !- handling the first layer : central sphere + rad = radin + rad(1) = 0.d0 + + do j = 1,nl + + ro = 0.d0 + vp = 0.d0 + vs = 0.d0 + + do i = 1,13 + + r1 = 0.d0 + r2 = 0.d0 + if ((rad(j) > r(i)).and.(rad(j) <= r(i+1))) then + if (rad(j+1) <= r(i+1)) then + r2 = rad(j+1) + r1 = rad(j) + else + r2 = r(i+1) + r1 = rad(j) + end if + else if (rad(j) <= r(i)) then + if ((rad(j+1) > r(i)).and.(rad(j+1) <= r(i+1))) then + r2 = rad(j+1) + r1 = r(i) + else if (rad(j+1) > r(i+1)) then + r2 = r(i+1) + r1 = r(i) + end if + end if + + t1 = d(i,1)/3.d0 + t2 = d(i,2)/(Rt*4.d0) + t3 = d(i,3)/((Rt**2)*5.d0) + t4 = d(i,4)/((Rt**3)*6.d0) + ro = ro + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - & + ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) ) + + t1 = p(i,1)/3.d0 + t2 = p(i,2)/(Rt*4.d0) + t3 = p(i,3)/((Rt**2)*5.d0) + t4 = p(i,4)/((Rt**3)*6.d0) + vp = vp + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - & + ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) ) + + t1 = s(i,1)/3.d0 + t2 = s(i,2)/(Rt*4.d0) + t3 = s(i,3)/((Rt**2)*5.d0) + t4 = s(i,4)/((Rt**3)*6.d0) + vs = vs + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - & + ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) ) + + end do + ro = ro*3 / (rad(j+1)**3-rad(j)**3) + vp = vp*3 /(rad(j+1)**3-rad(j)**3) + vs = vs*3 / (rad(j+1)**3-rad(j)**3) + mu = ro*vs**2 + la = ro*vp**2 - 2.d0*mu + ro = ro*1.d3 + la = la*1.d9 + mu = mu*1.d9 + + !print*, rad(j+1), ro, la , mu + rot(j) = ro + lat(j) = la + mut(j) = mu + end do + + + end subroutine earth_homogeneous_layers_prembased + + + subroutine earth_nlayers(x0,ro,la,mu,GR,solid) + use lnb_param +!--------------------------------------- +! Gets planet properties at non-dimentionalized radius x0: +!ro: density +!la: Lame lambda constant +!mu: Lame mu constant (shear modulus) +!GR: gravity +!solid: boolean, true if solid, otherwise liquid + integer :: i,j + doubleprecision :: la,mu,CST,GR + doubleprecision :: x,x0,xx0,ro + doubleprecision ::r1,r2 + logical :: solid + double precision, dimension(nlayer+1) :: r + + + + solid = .true. + + r(:) = radius(:) + r(1) = 0.d0 + x = x0 * ra + CST = 4.d0*pi*GG/3.d0 + if ( x > r(nlayer+1) ) x = r(nlayer+1) + + GR = 0.d0 + + do i = 1,nlayer + + if (x > r(i)) then + + r2 = r(i+1) + r1 = r(i) + ro = roc(i) + + if (x <= r(i+1)) then + + GR = GR + ro*(x**3-r1**3) + if (x > epsdb) GR = GR*CST/(x**2) + + else + GR = GR + ro*(r2**3-r1**3) + end if + end if + + if ((x >= r(i)).and.(x < r(i+1))) then + la = lac(i) + mu = muc(i) + ro = roc(i) + if ( soliddim(i).eqv..false. ) solid = .false. + + end if + + end do + + if (x>=r(nlayer+1)) then + la = lac(nlayer) + mu = muc(nlayer) + ro = roc(nlayer) + if ( soliddim(nlayer).eqv..false. ) solid = .false. + + end if + + end subroutine earth_nlayers + + subroutine earth_nlayers_viscoelas(x0,ro,la,mu,GR,solid) + use lnb_param +!--------------------------------------- +! Defines the earth rheology at non-dimensionalized radius x0 and angular frequency freq (global/public variable) + + integer :: i + double complex :: la, mu, cst + double precision :: GR,la0,mu00,ka,mu1,mu2,vi2 + double precision :: x,x0,xx0,ro,vi + logical :: notfound,solid,burg + + + + call earth_nlayers(x0,ro,la0,mu00,GR,solid) + + x = x0 * ra + do i = 1,nlayer + if ((x >= radius(i)).and.(x < radius(i+1))) then + vi = vic(i) + vi2 = vic2(i) + mu2=muc2(i) + burg=burgers(i) + end if + end do + + if (burg.eqv..true.) then + mu1=mu00 + + if ((vi*freq) /= 0.d0) then + ka=la0 + 2.d0/3.d0*mu00 + !print*,mu1,mu2,vi,vi2 + mu=mu1*freq*(freq+mu2/vi2)/((freq+mu2/vi2)*(freq+mu1/vi)+mu1/vi2*freq) + !la=4*la0*mu**2/(2*mu*(3*la0+2*mu)-3*la0*mu) + la=ka-2.d0/3.d0*mu + !print*,freq,mu,la + else + mu = mu00 + la = la0 + endif + else + ka = la0 + 2.d0/3.d0*mu00 + cst = vi*freq + if (cst /= 0.d0) then + mu = mu00/(1.d0+mu00/cst) + la = (la0 + mu00*ka/cst)/(1.d0 + mu00/cst) + else + mu = mu00 + la = la0 + endif + endif + + ! + !!write(101,*), freq,cst,la0,la,mu00,mu,vi + + end subroutine earth_nlayers_viscoelas + + + + +end module model Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/util.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/util.f90 (revision 22004) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/util.f90 (revision 22004) @@ -0,0 +1,226 @@ +module util + + use model + use lnb_param + + contains + +! +!_______DERIVS_________________________________________________________________ +! + subroutine derivs(x,y,dydx) + + implicit none + double complex :: la, mu + double precision :: x,ro,g + double complex :: flm,rlm,rm0,rg0,frh,fgr,fn + logical :: solid + integer :: iy,id,n2,ny + double complex, dimension(6) :: y,dydx,f + double complex, dimension(6,6) :: dy + + call earth_nlayers_viscoelas(x,ro,la,mu,g,solid) + !print*,la,mu +! +!_______Expressions +! + flm=(la+2.d0*mu) + rlm=(3.d0*la+2.d0*mu)/(la+2.d0*mu) + rm0=mu/mu0 + rg0=g/go0 + frh=ro*g*ra/mu0 + fgr=4.d0*Pi*GG*ro*ra/go0 + fn=dble(deg*(deg+1)) + +! +!_______Systeme differentiel yi dans un solide +! + if(solid) then + + ny = 6 + + dy(1,1)= (-2.d0*la/flm)/x + dy(1,2)= mu0/flm + dy(1,3)= (la*fn/flm)/x + dy(1,4)= 0.d0 + dy(1,5)= 0.d0 + dy(1,6)= 0.d0 + + dy(2,1)= 4.d0*(-frh+rm0*rlm/x)/x + ro*freq**2*ra/mu0 + dy(2,2)=(-4.d0*mu/flm)/x + dy(2,3)= fn*(frh-2.d0*rm0*rlm/x)/x + dy(2,4)= fn/x + dy(2,5)= 0.d0 + dy(2,6)= -frh/rg0 + + dy(3,1)= -1.d0/x + dy(3,2)= 0.d0 + dy(3,3)= 1.d0/x + dy(3,4)= 1.d0/rm0 + dy(3,5)= 0.d0 + dy(3,6)= 0.d0 + + dy(4,1)= (frh-2.d0*rm0*rlm/x)/x + dy(4,2)= ( -la/flm)/x + dy(4,3)= (2.d0*rm0*(la*(2.d0*fn-1.d0)+2.d0*mu*(fn-1.d0))/flm)/(x*x) + ro*freq**2*ra/mu0 + dy(4,4)= -3.d0/x + dy(4,5)= -(frh/rg0)/x + dy(4,6)= 0.d0 + + dy(5,1)= fgr + dy(5,2)= 0.d0 + dy(5,3)= 0.d0 + dy(5,4)= 0.d0 + dy(5,5)= 0.d0 + dy(5,6)= 1.d0 + + dy(6,1)= 0.d0 + dy(6,2)= 0.d0 + dy(6,3)= -(fgr*fn)/x + dy(6,4)= 0.d0 + dy(6,5)= fn/(x*x) + dy(6,6)= -2.d0/x + + else + + ny = 2 + + dy(:,:)= 0.d0 + + dy(1,1)= fgr/rg0 + dy(1,2)= 1.d0 + dy(2,1)= (-4.d0*(fgr/rg0)+fn/x)/x + dy(2,2)= -2.d0/x-fgr/rg0 + + end if +! +!_______Derivees des yi +! + dydx(:)=0.d0 + do id=1,ny + do iy=1,ny + dydx(id)=dydx(id)+dy(id,iy)*y(iy) + end do + end do + + end subroutine derivs + + + subroutine propagFDS(ystart,xmin,xmax,nstep) + implicit none + double precision :: xmin,xmax,dr,x + integer :: nstep,i + double complex, dimension(6) :: ystart,dydx,y,y1 + + dr = (xmax -xmin)/nstep + y = ystart + x = xmin + + do i = 1,nstep + + call derivs(x,y,dydx) + y1 = y + dr*dydx + x = x + dr + y = y1 + + end do + + ystart = y + + end subroutine propagFDS + +!_____LUDCMP___________________________________________________________________ +! + SUBROUTINE ludcmp(a,n,np,indx,d) + IMPLICIT DOUBLE PRECISION (a-h,o-z) + INTEGER n,np,indx(n),NMAX + REAL*8 d,a(np,np),TINY + PARAMETER (NMAX=500,TINY=1.0D-30) + INTEGER i,imax,j,k + REAL*8 aamax,dum,sum,vv(NMAX) + d=1d0 + do 12 i=1,n + aamax=0d0 + do 11 j=1,n + if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j)) +11 continue + if (aamax.eq.0d0) stop 'singular matrix in ludcmp' + vv(i)=1./aamax +12 continue + do 19 j=1,n + do 14 i=1,j-1 + sum=a(i,j) + do 13 k=1,i-1 + sum=sum-a(i,k)*a(k,j) +13 continue + a(i,j)=sum +14 continue + aamax=0d0 + do 16 i=j,n + sum=a(i,j) + do 15 k=1,j-1 + sum=sum-a(i,k)*a(k,j) +15 continue + a(i,j)=sum + dum=vv(i)*dabs(sum) + if (dum.ge.aamax) then + imax=i + aamax=dum + endif +16 continue + if (j.ne.imax)then + do 17 k=1,n + dum=a(imax,k) + a(imax,k)=a(j,k) + a(j,k)=dum +17 continue + d=-d + vv(imax)=vv(j) + endif + indx(j)=imax + if(a(j,j).eq.0d0)a(j,j)=TINY + if(j.ne.n)then + dum=1d0/a(j,j) + do 18 i=j+1,n + a(i,j)=a(i,j)*dum +18 continue + endif +19 continue + return + END SUBROUTINE ludcmp +! +!_____LUBKSB___________________________________________________________________ +! + SUBROUTINE lubksb(a,n,np,indx,b) + IMPLICIT DOUBLE PRECISION (a-h,o-z) + INTEGER n,np,indx(n) + REAL*8 a(np,np),b(np) + INTEGER i,ii,j,ll + REAL*8 sum + ii=0 + do 12 i=1,n + ll=indx(i) + sum=b(ll) + b(ll)=b(i) + if (ii.ne.0)then + do 11 j=ii,i-1 + sum=sum-a(i,j)*b(j) +11 continue + else if (sum.ne.0d0) then + ii=i + endif + b(i)=sum +12 continue + do 14 i=n,1,-1 + sum=b(i) + do 13 j=i+1,n + sum=sum-a(i,j)*b(j) +13 continue + b(i)=sum/a(i,i) +14 continue + return + END SUBROUTINE lubksb + + + +end module util Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22004) @@ -8,5 +8,5 @@ #include "./ModelProcessorx.h" -void CreateElementsVerticesAndMaterials(Elements* elements,Vertices* vertices,Materials* materials,IoModel* iomodel,const int nummodels){ +void CreateElementsVerticesAndMaterials(Elements* elements,Vertices* vertices,Materials* materials,IoModel* iomodel,const int nummodels,int solution_type){ /*Intermediary*/ @@ -15,4 +15,6 @@ bool control_analysis; bool dakota_analysis; + int nnat,dummy; + int* nature=NULL; /*Fetch parameters: */ @@ -23,4 +25,7 @@ /*Did we already create the elements? : */ _assert_(elements->Size()==0); + + /*Setup matpar counter in iomodel before we call on element constructors: */ + iomodel->matparcounter=iomodel->numberofelements+1; /*Create elements*/ @@ -112,4 +117,100 @@ } break; + case MaterialsEnum: + + //we have several types of materials. Retrieve this info first: + iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature"); + + //make sure materials that are not tied to elements come last: for now, only Matlitho qualifies. + for(int i=0;iFetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); + for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); + switch(iomodel->domaindim){ + case 2: + elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + break; + case 3: + break; + default: + _error_("Mesh not supported yet"); + } + break; + case MatlithoEnum: + iomodel->FetchData(9,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.isburgers","md.materials.issolid","md.materials.density"); + materials->AddObject(new Matlitho(materials->Size()+1,iomodel)); + iomodel->DeleteData(9,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.isburgers","md.materials.issolid","md.materials.density"); + break; + + case MatenhancediceEnum: + iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum); + for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); + switch(iomodel->domaindim){ + case 2: + elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + elements->InputDuplicate(MaterialsRheologyEEnum,MaterialsRheologyEbarEnum); + break; + case 3: + break; + default: + _error_("Mesh not supported yet"); + } + break; + case MatdamageiceEnum: + iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,"md.damage.D",DamageDEnum); + for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); + switch(iomodel->domaindim){ + case 2: + elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + elements->InputDuplicate(DamageDEnum,DamageDbarEnum); + break; + case 3: + break; + default: + _error_("Mesh not supported yet"); + } + break; + case MatestarEnum: + iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum); + for(i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matestar(i+1,i,iomodel)); + switch(iomodel->domaindim){ + case 2: + elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + elements->InputDuplicate(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum); + elements->InputDuplicate(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum); + break; + case 3: + break; + default: + _error_("Mesh not supported yet"); + } + break; + + default: + _error_("Materials "<(nature); + break; + default: _error_("Materials "<AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials + materials->AddObject(new Matpar(iomodel));//put it at the end of the materials /*Create vertices: */ @@ -129,5 +230,5 @@ if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r"); - CreateNumberNodeToElementConnectivity(iomodel); + CreateNumberNodeToElementConnectivity(iomodel,solution_type); for(i=0;inumberofvertices;i++){ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp (revision 22004) @@ -14,5 +14,5 @@ #include "./ModelProcessorx.h" -void CreateNumberNodeToElementConnectivity(IoModel* iomodel){ +void CreateNumberNodeToElementConnectivity(IoModel* iomodel,int solution_type){ /*Intermediary*/ @@ -28,25 +28,32 @@ /*Some checks if debugging*/ - _assert_(iomodel->numberofvertices); - _assert_(iomodel->numberofelements); - _assert_(iomodel->elements); + if (solution_type==LoveSolutionEnum){ + /*do nothing, we don't have a mesh. Just initialize to NULL*/ + } + else{ - /*Allocate ouput*/ - connectivity=xNewZeroInit(iomodel->numberofvertices); + /*Some checks if debugging*/ + _assert_(iomodel->numberofvertices); + _assert_(iomodel->numberofelements); + _assert_(iomodel->elements); - /*Get element width*/ - switch(iomodel->meshelementtype){ - case TriaEnum: elementswidth=3; break; - case TetraEnum: elementswidth=4; break; - case PentaEnum: elementswidth=6; break; - default: _error_("mesh not supported yet"); - } + /*Allocate ouput*/ + connectivity=xNewZeroInit(iomodel->numberofvertices); - /*Create connectivity table*/ - for (i=0;inumberofelements;i++){ - for (j=0;jelements[elementswidth*i+j]; - _assert_(vertexid>0 && vertexid-1numberofvertices); - connectivity[vertexid-1]+=1; + /*Get element width*/ + switch(iomodel->meshelementtype){ + case TriaEnum: elementswidth=3; break; + case TetraEnum: elementswidth=4; break; + case PentaEnum: elementswidth=6; break; + default: _error_("mesh not supported yet"); + } + + /*Create connectivity table*/ + for (i=0;inumberofelements;i++){ + for (j=0;jelements[elementswidth*i+j]; + _assert_(vertexid>0 && vertexid-1numberofvertices); + connectivity[vertexid-1]+=1; + } } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22004) @@ -23,4 +23,9 @@ char* fieldname = NULL; IssmDouble time; + + /*love numbers: */ + IssmDouble* frequencies = NULL; + int nfreq=0; + int dummy; /*parameters for mass flux:*/ @@ -61,4 +66,24 @@ parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); + + /*gia: */ + parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum)); + + /*love numbers: */ + iomodel->FetchData(&nfreq,"md.love.nfreq"); + parameters->AddObject(new IntParam(LoveNfreqEnum,nfreq)); + + iomodel->FetchData(&frequencies,&dummy,&dummy,"md.love.frequencies"); + parameters->AddObject(new DoubleVecParam(LoveFrequenciesEnum,frequencies,nfreq)); + xDelete(frequencies); + + parameters->AddObject(iomodel->CopyConstantObject("md.love.sh_nmax",LoveShNmaxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.sh_nmin",LoveShNminEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.g0",LoveG0Enum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.r0",LoveR0Enum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.mu0",LoveMu0Enum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.allow_layer_deletion",LoveAllowLayerDeletionEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.forcing_type",LoveForcingTypeEnum)); + {/*This is specific to ice...*/ parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum)); @@ -80,5 +105,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.transient.iscoupler",TransientIscouplerEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum)); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 22004) @@ -36,5 +36,5 @@ /*Create elements, vertices and materials, independent of analysis_enum: */ - CreateElementsVerticesAndMaterials(elements,vertices,materials,iomodel,nummodels); + CreateElementsVerticesAndMaterials(elements,vertices,materials,iomodel,nummodels,solution_enum); /*Create Parameters*/ @@ -68,6 +68,8 @@ /*Solution specific updates*/ + if(VerboseMProcessor()) _printf0_(" updating elements and materials for control parameters" << "\n"); UpdateElementsAndMaterialsControl(elements,materials,iomodel); #ifdef _HAVE_DAKOTA_ + if(VerboseMProcessor()) _printf0_(" updating elements and materials for uncertainty quantification" << "\n"); UpdateElementsAndMaterialsDakota(elements,materials,iomodel); #endif @@ -77,4 +79,5 @@ /*Output definitions dataset: */ + if(VerboseMProcessor()) _printf0_(" creating output definitions" << "\n"); CreateOutputDefinitions(elements,parameters,iomodel); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 22003) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 22004) @@ -12,5 +12,5 @@ /*Creation of fem datasets: general drivers*/ -void CreateElementsVerticesAndMaterials(Elements* elements,Vertices* vertices,Materials* materials, IoModel* iomodel,const int nummodels); +void CreateElementsVerticesAndMaterials(Elements* elements,Vertices* vertices,Materials* materials, IoModel* iomodel,const int nummodels,int solution_type); void CreateParameters(Parameters*parameters,IoModel* iomodel,char* rootpath,FILE* toolkitfile,const int solution_type); void CreateParametersAutodiff(Parameters* parameters,IoModel* iomodel); @@ -37,4 +37,4 @@ /*Connectivity*/ void CreateSingleNodeToElementConnectivity(IoModel* iomodel); -void CreateNumberNodeToElementConnectivity(IoModel* iomodel); +void CreateNumberNodeToElementConnectivity(IoModel* iomodel,int solution_type); #endif Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 22004) @@ -70,5 +70,5 @@ /*Initialize vector with number of vertices*/ - int numvertices=femmodel->vertices->NumberOfVertices(); + int numvertices=femmodel->vertices->NumberOfVertices(); if (numvertices==0)return; Vector* vec_mask_ice=new Vector(numvertices); //vertices that have ice at next time step /*Fill vector with values: */ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 22003) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 22004) @@ -30,4 +30,5 @@ #include "./GetVectorFromControlInputsx/GetVectorFromControlInputsx.h" #include "./GiaDeflectionCorex/GiaDeflectionCorex.h" +#include "./FourierLoveCorex/FourierLoveCorex.h" #include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h" #include "./SetActiveNodesLSMx/SetActiveNodesLSMx.h" Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22003) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22004) @@ -118,4 +118,6 @@ case 3: return MaticeEnum; case 4: return MatenhancediceEnum; + case 5: return MatlithoEnum; + case 6: return MaterialsEnum; default: _error_("Marshalled materials code \""<',0); + md = checkfield(md,'fieldname','love.frequencies','NaN',1,'Inf',1,'numel',md.love.nfreq); + md = checkfield(md,'fieldname','love.sh_nmax','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.sh_nmin','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.g0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.r0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.mu0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.allow_layer_deletion','values',[0 1]); + md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',1,'>',0, '<=', 12); + if md.love.sh_nmin<=1 & md.love.forcing_type==9 + error('Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.') + end + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'fieldname','sh_nmax','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','sh_nmin','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','g0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','r0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','mu0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','allow_layer_deletion','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','forcing_type','format','Integer'); + + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + error('not implemented yet!'); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/fourierlove.py =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 22004) +++ /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 22004) @@ -0,0 +1,104 @@ +from fielddisplay import fielddisplay +from checkfield import checkfield +from WriteData import WriteData + +class fourierlove(object): + """ + Fourier Love Number class definition + + Usage: + fourierlove=fourierlove(); + """ + + def __init__(self): # {{{ + self.nfreq = float('NaN'); + self.frequencies = float('NaN'); + self.sh_nmax = float('NaN'); + self.sh_nmin = float('NaN'); + self.g0 = float('NaN'); + self.r0 = float('NaN'); + self.mu0 = float('NaN'); + self.allow_layer_deletion = float('NaN'); + self.forcing_type = float('NaN'); + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + + string=' Fourier Love class:' + + string="%s\n%s"%(string,fielddisplay(self,'nfreq','number of frequencies sampled (default 1, elastic) [Hz]')) + string="%s\n%s"%(string,fielddisplay(self,'frequencies','frequencies sampled (convention defaults to 0 for the elastic case) [Hz]')) + string="%s\n%s"%(string,fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)')) + string="%s\n%s"%(string,fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default 1)')) + string="%s\n%s"%(string,fielddisplay(self,'g0','adimensioning constant for gravity (default 10) [m/s^2]')) + string="%s\n%s"%(string,fielddisplay(self,'r0','adimensioning constant for radius (default 6378*10^3) [m]')) + string="%s\n%s"%(string,fielddisplay(self,'mu0','adimensioning constant for stress (default 10^11) [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)')) + string="%s\n%s"%(string,fielddisplay(self,'forcing_type','integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :')) + string="%s\n%s"%(string,' 1: Inner core boundary -- Volumic Potential') + string="%s\n%s"%(string,' 2: Inner core boundary -- Pressure') + string="%s\n%s"%(string,' 3: Inner core boundary -- Loading') + string="%s\n%s"%(string,' 4: Inner core boundary -- Tangential traction') + string="%s\n%s"%(string,' 5: Core mantle boundary -- Volumic Potential') + string="%s\n%s"%(string,' 6: Core mantle boundary -- Pressure') + string="%s\n%s"%(string,' 7: Core mantle boundary -- Loading') + string="%s\n%s"%(string,' 8: Core mantle boundary -- Tangential traction') + string="%s\n%s"%(string,' 9: Surface -- Volumic Potential') + string="%s\n%s"%(string,' 10: Surface -- Pressure') + string="%s\n%s"%(string,' 11: Surface -- Loading') + string="%s\n%s"%(string,' 12: Surface -- Tangential traction ') + + return string; + + + #}}} + def extrude(self,md): # {{{ + return self + #}}} + def setdefaultparameters(self): # {{{ + + #we setup an elastic love number computation by default. + self.nfreq=1 + self.frequencies=[0]; #Hz + self.sh_nmax=256 # .35 degree, 40 km at the equator. + self.sh_nmin=1 + self.g0=10 # m/s^2 + self.r0=6378*1e3 #m + self.mu0=10^11 # Pa + self.allow_layer_deletion=1 + self.forcing_type = 11 + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + md = checkfield(md,'fieldname','love.nfreq','NaN',1,'Inf',1,'numel',[1],'>',0); + md = checkfield(md,'fieldname','love.frequencies','NaN',1,'Inf',1,'numel',[md.love.nfreq]); + md = checkfield(md,'fieldname','love.sh_nmax','NaN',1,'Inf',1,'numel',[1],'>',0); + md = checkfield(md,'fieldname','love.sh_nmin','NaN',1,'Inf',1,'numel',[1],'>',0); + md = checkfield(md,'fieldname','love.g0','NaN',1,'Inf',1,'numel',[1],'>',0); + md = checkfield(md,'fieldname','love.r0','NaN',1,'Inf',1,'numel',[1],'>',0); + md = checkfield(md,'fieldname','love.mu0','NaN',1,'Inf',1,'numel',[1],'>',0); + md = checkfield(md,'fieldname','love.allow_layer_deletion','values',[0,1]); + md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',[1],'>',0, '<=', 12); + if md.love.sh_nmin<=1 and md.love.forcing_type==9: + raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.") + + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + + WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'fieldname','sh_nmax','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','sh_nmin','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','g0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','r0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','mu0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','allow_layer_deletion','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','forcing_type','format','Integer'); + + # }}} Index: /issm/trunk-jpl/src/m/classes/maskpsl.m =================================================================== --- /issm/trunk-jpl/src/m/classes/maskpsl.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/maskpsl.m (revision 22004) @@ -42,4 +42,6 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ + + if strcmpi(solution,'LoveSolution'), return; end; md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices 1]); Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 22004) @@ -113,16 +113,16 @@ self.numlayers=2; - %surface, then the lab (lithosphere/asthenosphere boundary) then the center of the earth + %center of the earth (approximation, must not be 0), then the lab (lithosphere/asthenosphere boundary) then the surface %(with 1d3 to avoid numerical singularities) - self.radius=[6378*1e3;6278*1e3;1e3]; - - self.viscosity=[1e40;1e21]; %lithosphere and mantle viscosity (respectively) [Pa.s] - self.lame_mu=[6.7*10^10;1.45*1e11]; % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa] - self.lame_lambda=[6.7*10^10;1.45*1e11]; % (Pa) %lithosphere and mantle lamba parameter (respectively) [Pa] + self.radius=[1e3;6278*1e3;6378*1e3]; + + self.viscosity=[1e21;1e40]; %mantle and lithosphere viscosity (respectively) [Pa.s] + self.lame_mu=[1.45*1e11;6.7*10^10]; % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa] + self.lame_lambda=self.lame_mu; % (Pa) %mantle and lithosphere lamba parameter (respectively) [Pa] self.burgers_viscosity=[NaN;NaN]; self.burgers_mu=[NaN;NaN]; - self.isburgers=[false;false]; - self.density=[3.32*1e3;3.34*1e3]; % (Pa) %lithosphere and mantle density [kg/m^3] - self.issolid=[true;true]; % is layer solid or liquid. + self.isburgers=[0;0]; + self.density=[5.51*1e3;5.50*1e3]; % (Pa) %mantle and lithosphere density [kg/m^3] + self.issolid=[1;1]; % is layer solid or liquid. otherwise @@ -188,37 +188,30 @@ md = checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',1); md = checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1 1],'>',0); - md = checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); - md = checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); - md = checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0,'<',2); + md = checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<',2); md = checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); - md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); - md = checkfield(md,'fieldname','materials.isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0,'<',2); - md = checkfield(md,'fieldname','materials.burgers_viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); - md = checkfield(md,'fieldname','materials.burgers_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); - - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); - end - end - - end % }}} - function intnat = naturetointeger(strnat) % {{{ - intnat=zeros(length(strnat),1); - for i=1:length(strnat), - switch strnat{i}, - case 'damageice' - intnat(i)=1; - case 'estar' - intnat(i)=2; - case 'ice' - intnat(i)=3; - case 'enhancedice' - intnat(i)=4; - case 'litho' - intnat(i)=5; - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); - end - end + md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',1); + md = checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers 1],'>=',0); + + for i=1:md.materials.numlayers, + if md.materials.isburgers(i) & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))), + error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice'); + end + end + if md.materials.issolid(1)==0 | md.materials.lame_mu(1)==0 + error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.'); + end + ind=find(md.materials.issolid==0); + if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0 + error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.']) + end + otherwise + error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); + end + end + end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -226,5 +219,5 @@ %1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum WriteData(fid,prefix,'name','md.materials.type','data',6,'format','Integer'); - WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',1); + WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3); for i=1:length(self.nature), @@ -248,14 +241,14 @@ WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String'); case 'litho' - WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda ','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','BooleanMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','BooleanMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Integer'); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3); otherwise error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); @@ -316,4 +309,21 @@ end end - - + function intnat = naturetointeger(strnat) % {{{ + intnat=zeros(length(strnat),1); + for i=1:length(strnat), + switch strnat{i}, + case 'damageice' + intnat(i)=1; + case 'estar' + intnat(i)=2; + case 'ice' + intnat(i)=3; + case 'enhancedice' + intnat(i)=4; + case 'litho' + intnat(i)=5; + otherwise + error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); + end + end + end % }}} Index: /issm/trunk-jpl/src/m/classes/materials.py =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.py (revision 22004) +++ /issm/trunk-jpl/src/m/classes/materials.py (revision 22004) @@ -0,0 +1,280 @@ +import numpy as np +from fielddisplay import fielddisplay +from project3d import project3d +from checkfield import checkfield +from WriteData import WriteData + +def naturetointeger(strnat): #{{{ + + intnat=np.array(len(strnat)) + for i in range(len(intnat)): + if strnat[i]=='damageice': + intnat[i]=1 + elif strnat[i]=='estar': + intnat[i]=2 + elif strnat[i]=='ice': + intnat[i]=3 + elif strnat[i]=='enhancedice': + intnat[i]=4 + elif strnat[i]=='litho': + intnat[i]=5 + else: + raise RuntimeError("materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')"); + + return intnat +# }}} +class materials(object): + """ + MATERIALS class definition + + Usage: + materials=materials(); + """ + + def __init__(self,*args): # {{{ + + self.nature = [] + + if not len(args): + self.nature=['ice'] + else: + self.nature=args + + for i in range(len(self.nature)): + if not(self.nature[i] == 'litho' or self.nature[i]=='ice'): + raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho')") + + #start filling in the dynamic fields: + for i in range(len(self.nature)): + nat=self.nature[i]; + if nat=='ice': + setattr(self,'rho_ice',0) + setattr(self,'rho_ice',0); + setattr(self,'rho_water',0); + setattr(self,'rho_freshwater',0); + setattr(self,'mu_water',0); + setattr(self,'heatcapacity',0); + setattr(self,'latentheat',0); + setattr(self,'thermalconductivity',0); + setattr(self,'temperateiceconductivity',0); + setattr(self,'meltingpoint',0); + setattr(self,'beta',0); + setattr(self,'mixed_layer_capacity',0); + setattr(self,'thermal_exchange_velocity',0); + setattr(self,'rheology_B',0); + setattr(self,'rheology_n',0); + setattr(self,'rheology_law',0); + elif nat=='litho': + setattr(self,'numlayers',0); + setattr(self,'radius',0); + setattr(self,'viscosity',0); + setattr(self,'lame_lambda',0); + setattr(self,'lame_mu',0); + setattr(self,'burgers_viscosity',0); + setattr(self,'burgers_mu',0); + setattr(self,'isburgers',0); + setattr(self,'density',0); + setattr(self,'issolid',0); + else: + raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho')"); + #set default parameters: + self.setdefaultparameters() + #}}} + def __repr__(self): # {{{ + string=" Materials:" + for i in range(len(self.nature)): + nat=self.nature[i]; + if nat=='ice': + string="%s\n%s"%(string,'Ice:'); + string="%s\n%s"%(string,fielddisplay(self,"rho_ice","ice density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"rho_water","water density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"rho_freshwater","fresh water density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"mu_water","water viscosity [N s/m^2]")) + string="%s\n%s"%(string,fielddisplay(self,"heatcapacity","heat capacity [J/kg/K]")) + string="%s\n%s"%(string,fielddisplay(self,"thermalconductivity","ice thermal conductivity [W/m/K]")) + string="%s\n%s"%(string,fielddisplay(self,"temperateiceconductivity","temperate ice thermal conductivity [W/m/K]")) + string="%s\n%s"%(string,fielddisplay(self,"meltingpoint","melting point of ice at 1atm in K")) + string="%s\n%s"%(string,fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"beta","rate of change of melting point with pressure [K/Pa]")) + string="%s\n%s"%(string,fielddisplay(self,"mixed_layer_capacity","mixed layer capacity [W/kg/K]")) + string="%s\n%s"%(string,fielddisplay(self,"thermal_exchange_velocity","thermal exchange velocity [m/s]")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_B","flow law parameter [Pa/s^(1/n)]")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_n","Glen's flow law exponent")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_law","law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) + elif nat=='litho': + string="%s\n%s"%(string,'Litho:'); + string="%s\n%s"%(string,fielddisplay(self,'numlayers','number of layers (default 2)')) + string="%s\n%s"%(string,fielddisplay(self,'radius','array describing the radius for each interface (numlayers+1) [m]')) + string="%s\n%s"%(string,fielddisplay(self,'viscosity','array describing each layer''s viscosity (numlayers) [Pa.s]')) + string="%s\n%s"%(string,fielddisplay(self,'lame_lambda','array describing the lame lambda parameter (numlayers) [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'lame_mu','array describing the shear modulus for each layers (numlayers) [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'burgers_viscosity','array describing each layer''s transient viscosity, only for Burgers rheologies (numlayers) [Pa.s]')) + string="%s\n%s"%(string,fielddisplay(self,'burgers_mu','array describing each layer''s transient shear modulus, only for Burgers rheologies (numlayers) [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'isburgers','array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)')) + string="%s\n%s"%(string,fielddisplay(self,'density','array describing each layer''s density (numlayers) [kg/m^3]')) + string="%s\n%s"%(string,fielddisplay(self,'issolid','array describing whether the layer is solid or liquid (default 1) (numlayers)')) + + else: + raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho')"); + + return string + #}}} + def extrude(self,md): # {{{ + for i in range(len(self.nature)): + nat=self.nature[i]; + if nat=='ice': + self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node') + self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element') + return self + #}}} + def setdefaultparameters(self): # {{{ + for i in range(len(self.nature)): + nat=self.nature[i]; + if nat=='ice': + #ice density (kg/m^3) + self.rho_ice=917. + + #ocean water density (kg/m^3) + self.rho_water=1023. + + #fresh water density (kg/m^3) + self.rho_freshwater=1000. + + #water viscosity (N.s/m^2) + self.mu_water=0.001787 + + #ice heat capacity cp (J/kg/K) + self.heatcapacity=2093. + + #ice latent heat of fusion L (J/kg) + self.latentheat=3.34*10^5 + + #ice thermal conductivity (W/m/K) + self.thermalconductivity=2.4 + + #wet ice thermal conductivity (W/m/K) + self.temperateiceconductivity=.24 + + #the melting point of ice at 1 atmosphere of pressure in K + self.meltingpoint=273.15 + + #rate of change of melting point with pressure (K/Pa) + self.beta=9.8*10^-8 + + #mixed layer (ice-water interface) heat capacity (J/kg/K) + self.mixed_layer_capacity=3974. + + #thermal exchange velocity (ice-water interface) (m/s) + self.thermal_exchange_velocity=1.00*10^-4 + + #Rheology law: what is the temperature dependence of B with T + #available: none, paterson and arrhenius + self.rheology_law='Paterson' + + elif nat=='litho': + #we default to a configuration that enables running GIA solutions using giacaron and/or giaivins. + self.numlayers=2 + + #center of the earth (approximation, must not be 0), then the lab (lithosphere/asthenosphere boundary) then the surface + #(with 1d3 to avoid numerical singularities) + self.radius=[1e3,6278*1e3,6378*1e3] + + self.viscosity=[1e21,1e40] #mantle and lithosphere viscosity (respectively) [Pa.s] + self.lame_mu=[1.45*1e11,6.7*10^10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa] + self.lame_lambda=self.lame_mu # (Pa) #mantle and lithosphere lamba parameter (respectively) [Pa] + self.burgers_viscosity=[NaN,NaN] + self.burgers_mu=[NaN,NaN] + self.isburgers=[0,0] + self.density=[5.51*1e3,5.50*1e3] # (Pa) #mantle and lithosphere density [kg/m^3] + self.issolid=[1,1] # is layer solid or liquid. + + else: + raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho')"); + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + for i in range(len(self.nature)): + nat=self.nature[i]; + if nat=='ice': + md = checkfield(md,'fieldname','materials.rho_ice','>',0) + md = checkfield(md,'fieldname','materials.rho_water','>',0) + md = checkfield(md,'fieldname','materials.rho_freshwater','>',0) + md = checkfield(md,'fieldname','materials.mu_water','>',0) + md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements]) + md = checkfield(md,'fieldname','materials.rheology_law','values',['None','BuddJacka','Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval']) + elif nat=='litho': + if 'LoveAnalysis' not in analyses: + return md + + md = checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',[1]) + md = checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1,1],'>',0) + md = checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers,1],'>=',0) + md = checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers,1],'>=',0) + md = checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers,1],'>=',0,'<',2) + md = checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers,1],'>',0) + md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers,1],'>=',0) + md = checkfield(md,'fieldname','materials.isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers,1],'>=',0,'<=',1) + md = checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers,1],'>=',0) + md = checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers,1],'>=',0) + + for i in range(md.materials.numlayers): + if md.materials.isburgers[i] and (np.isnan(md.materials.burgers_viscosity[i] or np.isnan(md.materials.burgers_mu[i]))): + raise RuntimeError("materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice") + + if md.materials.issolid[0]==0 or md.materials.lame_mu[0]==0: + raise RuntimeError('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.') + + for i in range(md.materials.numlayers-1): + if md.materials.issolid[i] and md.materials.issolid[i+1]: #if there are at least two consecutive indices that contain issolid = 0 + raise RuntimeError("%s%i%s"%('2 or more adjacent fluid layers detected starting at layer',i,'This is not supported yet. Consider merging them.')) + + else: + raise RuntimeError("materials checkconsistency error message: nature of the material not supported yet! ('ice' or 'litho')"); + + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + + #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum + WriteData(fid,prefix,'name','md.materials.type','data',6,'format','Integer') + WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3) + + for i in range(len(self.nature)): + nat=self.nature[i]; + if nat=='ice': + + WriteData(fid,prefix,'name','md.materials.type','data',3,'format','Integer') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String') + + elif nat=='litho': + WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Integer') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3) + + else: + raise RuntimeError("materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')") + + # }}} Index: /issm/trunk-jpl/src/m/classes/mesh2d.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh2d.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/mesh2d.m (revision 22004) @@ -80,4 +80,6 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ + + if strcmpi(solution,'LoveSolution'), return; end md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); Index: /issm/trunk-jpl/src/m/classes/mesh2dvertical.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh2dvertical.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/mesh2dvertical.m (revision 22004) @@ -79,4 +79,6 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ + + if strcmpi(solution,'LoveSolution'), return; end md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); Index: /issm/trunk-jpl/src/m/classes/mesh3dprisms.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dprisms.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/mesh3dprisms.m (revision 22004) @@ -90,4 +90,6 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ + + if strcmpi(solution,'LoveSolution'), return; end md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dsurface.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.m (revision 22004) @@ -80,4 +80,6 @@ end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ + + if strcmpi(solution,'LoveSolution'), return; end md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); Index: /issm/trunk-jpl/src/m/classes/mesh3dtetras.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dtetras.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/mesh3dtetras.m (revision 22004) @@ -90,4 +90,6 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ + + if strcmpi(solution,'LoveSolution'), return; end md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 22003) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 22004) @@ -40,5 +40,6 @@ levelset = 0; calving = 0; - gia = 0; + love = 0; + gia = 0; esa = 0; @@ -1150,5 +1151,6 @@ md.levelset = levelset(); md.calving = calving(); - md.gia = giaivins(); + md.gia = giaivins(); + md.love = fourierlove(); md.esa = esa(); md.autodiff = autodiff(); @@ -1324,4 +1326,5 @@ disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving')); disp(sprintf('%19s: %-22s -- %s','gia' ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution')); + disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution')); disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters')); Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 22003) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 22004) @@ -18,4 +18,5 @@ from levelset import levelset from calving import calving +from fourierlove import fourierlove from calvinglevermann import calvinglevermann #from calvingpi import calvingpi @@ -109,4 +110,5 @@ self.levelset = levelset() self.calving = calving() + self.love = fourierlove() self.gia = giaivins() @@ -153,5 +155,6 @@ 'levelset',\ 'calving',\ - 'gia',\ + 'gia',\ + 'love',\ 'autodiff',\ 'inversion',\ @@ -195,4 +198,5 @@ string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("levelset","[%s,%s]" % ("1x1",obj.levelset.__class__.__name__),"parameters for moving boundaries (level-set method)")) string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("calving","[%s,%s]" % ("1x1",obj.calving.__class__.__name__),"parameters for calving")) + string="%s\n%s" % (string,'%19s: %-22s -- %s' % ("love","[%s,%s]" % ("1x1",obj.love.__class__.__name__),"parameters for love solution")) string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("autodiff","[%s,%s]" % ("1x1",obj.autodiff.__class__.__name__),"automatic differentiation parameters")) string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("inversion","[%s,%s]" % ("1x1",obj.inversion.__class__.__name__),"parameters for inverse methods")) Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 22003) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 22004) @@ -66,4 +66,6 @@ elseif strcmp(solutiontype,'GiaSolution') analyses={'GiaIvinsAnalysis'}; + elseif strcmp(solutiontype,'LoveSolution') + analyses={'LoveAnalysis'}; elseif strcmp(solutiontype,'EsaSolution') analyses={'EsaAnalysis'}; Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 22003) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 22004) @@ -33,4 +33,7 @@ elif solutiontype == 'GiaSolution': analyses=['GiaIvinsAnalysis'] + + elif solutiontype == 'LoveSolution': + analyses=['LoveAnalysis'] elif solutiontype == 'TransientSolution': Index: /issm/trunk-jpl/src/m/solve/marshall.py =================================================================== --- /issm/trunk-jpl/src/m/solve/marshall.py (revision 22003) +++ /issm/trunk-jpl/src/m/solve/marshall.py (revision 22004) @@ -20,8 +20,4 @@ raise IOError("marshall error message: could not open '%s.bin' file for binary writing." % md.miscellaneous.name) - #Go through all model fields: check that it is a class and call checkconsistency - fields=vars(md) - -# for field in fields.iterkeys(): for field in md.properties(): Index: /issm/trunk-jpl/src/m/solve/solve.m =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.m (revision 22003) +++ /issm/trunk-jpl/src/m/solve/solve.m (revision 22004) @@ -21,4 +21,5 @@ % - 'Esa' or 'esa' % - 'Sealevelrise' or 'slr' +% - 'Love' or 'lv' % % extra options: @@ -64,4 +65,6 @@ elseif strcmpi(solutionstring,'gia') || strcmpi(solutionstring,'Gia') solutionstring = 'GiaSolution'; +elseif strcmpi(solutionstring,'lv') || strcmpi(solutionstring,'Love') + solutionstring = 'LoveSolution'; elseif strcmpi(solutionstring,'esa') || strcmpi(solutionstring,'Esa') solutionstring = 'EsaSolution'; Index: /issm/trunk-jpl/src/m/solve/solve.py =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.py (revision 22003) +++ /issm/trunk-jpl/src/m/solve/solve.py (revision 22004) @@ -30,4 +30,5 @@ - 'Gia' or 'gia' - 'Sealevelrise' or 'slr' + - 'Love' or 'lv' extra options: @@ -66,4 +67,6 @@ elif solutionstring.lower() == 'gia' or solutionstring.lower() == 'gia': solutionstring = 'GiaSolution'; + elif solutionstring.lower() == 'lv' or solutionstring.lower() == 'love': + solutionstring = 'LoveSolution'; elif solutionstring.lower() == 'slr' or solutionstring.lower() == 'sealevelrise': solutionstring = 'SealevelriseSolution';