Index: /issm/trunk-jpl/m4/analyses.m4 =================================================================== --- /issm/trunk-jpl/m4/analyses.m4 (revision 23935) +++ /issm/trunk-jpl/m4/analyses.m4 (revision 23936) @@ -304,4 +304,18 @@ AC_MSG_RESULT($HAVE_HYDROLOGYSHREVE) dnl }}} +dnl with-HydrologyGlaDS{{{ +AC_ARG_WITH([HydrologyGlaDS], + AS_HELP_STRING([--with-HydrologyGlaDS = YES], [compile with HydrologyGlaDS capabilities (default is yes)]), + [HYDROLOGYGLADS=$withval],[HYDROLOGYGLADS=yes]) +AC_MSG_CHECKING(for HydrologyGlaDS capability compilation) + +HAVE_HYDROLOGYGLADS=no +if test "x$HYDROLOGYGLADS" = "xyes"; then + HAVE_HYDROLOGYGLADS=yes + AC_DEFINE([_HAVE_HYDROLOGYGLADS_],[1],[with HydrologyGlaDS capability]) +fi +AM_CONDITIONAL([HYDROLOGYGLADS], [test x$HAVE_HYDROLOGYGLADS = xyes]) +AC_MSG_RESULT($HAVE_HYDROLOGYGLADS) +dnl }}} dnl with-HydrologyShakti{{{ AC_ARG_WITH([HydrologyShakti], Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 23935) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 23936) @@ -436,4 +436,7 @@ issm_sources += ./analyses/HydrologyPismAnalysis.cpp endif +if HYDROLOGYGLADS +issm_sources += ./analyses/HydrologyGlaDSAnalysis.cpp +endif if HYDROLOGYDCINEFFICIENT issm_sources += ./analyses/HydrologyDCInefficientAnalysis.cpp Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 23935) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 23936) @@ -77,4 +77,7 @@ case HydrologyShreveAnalysisEnum : return new HydrologyShreveAnalysis(); #endif + #ifdef _HAVE_HYDROLOGYGLADS_ + case HydrologyGlaDSAnalysisEnum : return new HydrologyGlaDSAnalysis(); + #endif #ifdef _HAVE_HYDROLOGYSHAKTI_ case HydrologyShaktiAnalysisEnum : return new HydrologyShaktiAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23936) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23936) @@ -0,0 +1,400 @@ +#include "./HydrologyGlaDSAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +void HydrologyGlaDSAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + if(hydrology_model!=HydrologyGlaDSEnum) return; + + _error_("Not implemented"); + //IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologyGlaDSAnalysisEnum,P1Enum); + +}/*}}}*/ +void HydrologyGlaDSAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int hydrology_model; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want GlaDS?*/ + if(hydrology_model!=HydrologyGlaDSEnum) return; + + _error_("Not implemented"); + + /*Create discrete loads for Moulins*/ + CreateSingleNodeToElementConnectivity(iomodel); + for(int i=0;inumberofvertices;i++){ + if (iomodel->domaintype!=Domain3DEnum){ + /*keep only this partition's nodes:*/ + if(iomodel->my_vertices[i]){ + loads->AddObject(new Moulin(i+1,i,iomodel,HydrologyGlaDSAnalysisEnum)); + } + } + else if(reCast(iomodel->Data("md.mesh.vertexonbase")[i])){ + if(iomodel->my_vertices[i]){ + loads->AddObject(new Moulin(i+1,i,iomodel,HydrologyGlaDSAnalysisEnum)); + } + } + } + iomodel->DeleteData(1,"md.mesh.vertexonbase"); + + /*Deal with Neumann BC*/ + int M,N; + int *segments = NULL; + iomodel->FetchData(&segments,&M,&N,"md.mesh.segments"); + + /*Check that the size seem right*/ + _assert_(N==3); _assert_(M>=3); + for(int i=0;imy_elements[segments[i*3+2]-1]){ + loads->AddObject(new Neumannflux(i+1,i,iomodel,segments,HydrologyGlaDSAnalysisEnum)); + } + } + + xDelete(segments); + +}/*}}}*/ +void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/ + + /*Fetch parameters: */ + int hydrology_model; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want GlaDS?*/ + if(hydrology_model!=HydrologyGlaDSEnum) return; + + _error_("Not implemented"); + + if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); + ::CreateNodes(nodes,iomodel,HydrologyGlaDSAnalysisEnum,P1Enum); + iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); +}/*}}}*/ +int HydrologyGlaDSAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void HydrologyGlaDSAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Fetch data needed: */ + int hydrology_model,frictionlaw; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want GlaDS?*/ + if(hydrology_model!=HydrologyGlaDSEnum) return; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum); + iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum); + if(iomodel->domaintype!=Domain2DhorizontalEnum){ + iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); + iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); + } + iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); + + _error_("not impelemented yet"); +}/*}}}*/ +void HydrologyGlaDSAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + int numoutputs; + char** requestedoutputs = NULL; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want GlaDS?*/ + if(hydrology_model!=HydrologyGlaDSEnum) return; + + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + //parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum)); + _error_("not implemented"); + + /*Requested outputs*/ + iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs"); + parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs"); +}/*}}}*/ + +/*Finite Element Analysis*/ +void HydrologyGlaDSAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* HydrologyGlaDSAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* HydrologyGlaDSAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ +ElementMatrix* HydrologyGlaDSAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + _error_("Not implemented"); + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* dbasis = xNew(2*numnodes); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Get englacial storage coefficient*/ + IssmDouble storage,dt; + element->FindParam(&storage,HydrologyStorageEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(1); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + for(int i=0;ivalues[i*numnodes+j] += gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]) + + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]; + } + } + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(dbasis); + delete gauss; + return Ke; +}/*}}}*/ +ElementVector* HydrologyGlaDSAnalysis::CreatePVector(Element* element){/*{{{*/ + + /*Skip if water or ice shelf element*/ + if(element->IsFloating()) return NULL; + + _error_("not implemented"); + + /*Intermediaries */ + IssmDouble Jdet,meltrate,G,dh[2],B,A,n; + IssmDouble gap,bed,thickness,head,ieb,head_old; + IssmDouble lr,br,vx,vy,beta; + IssmDouble alpha2,frictionheat; + IssmDouble PMPheat,dpressure_water[2],dbed[2]; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum); + IssmDouble g = element->FindParam(ConstantsGEnum); + IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); + Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); + Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); + Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); + Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); + + + /*Get englacial storage coefficient*/ + IssmDouble storage,dt; + element->FindParam(&storage,HydrologyStorageEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(element,2); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + geothermalflux_input->GetInputValue(&G,gauss); + base_input->GetInputValue(&bed,gauss); + base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss); + thickness_input->GetInputValue(&thickness,gauss); + gap_input->GetInputValue(&gap,gauss); + head_input->GetInputValue(&head,gauss); + head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); + englacial_input->GetInputValue(&ieb,gauss); + lr_input->GetInputValue(&lr,gauss); + br_input->GetInputValue(&br,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + headold_input->GetInputValue(&head_old,gauss); + + /*Get ice A parameter*/ + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + A=pow(B,-n); + + /*Compute beta term*/ + if(gapGetAlpha2(&alpha2,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + frictionheat=alpha2*(vx*vx+vy*vy); + + /*Get water and ice pressures*/ + IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); + IssmDouble pressure_water = rho_water*g*(head-bed); + if(pressure_water>pressure_ice) pressure_water = pressure_ice; + + /*Get water pressure from previous time step to use in lagged creep term*/ + IssmDouble pressure_water_old = rho_water*g*(head_old-bed); + if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice; + + /*Compute change in sensible heat due to changes in pressure melting point*/ + dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); + dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); + PMPheat=0.; + + meltrate = 1/latentheat*(G+frictionheat+rho_water*g*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); + _assert_(meltrate>0.); + + for(int i=0;ivalues[i]+=Jdet*gauss->weight* + ( + meltrate*(1/rho_water-1/rho_ice) + +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap + -beta*sqrt(vx*vx+vy*vy) + +ieb + +storage*head_old/dt + )*basis[i]; + } + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete friction; + delete gauss; + return pe; +}/*}}}*/ +void HydrologyGlaDSAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented"); + element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum); +}/*}}}*/ +void HydrologyGlaDSAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + + _error_("not implemented"); + + /*Intermediary*/ + IssmDouble dh[3]; + int* doflist = NULL; + IssmDouble* xyz_list = NULL; + + /*Get gravity from parameters*/ + IssmDouble g = element->FindParam(ConstantsGEnum); + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Fetch dof list and allocate solution vector*/ + element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numnodes); + + /*Get thickness and base on nodes to apply cap on water head*/ + IssmDouble* eff_pressure = xNew(numnodes); + IssmDouble* thickness = xNew(numnodes); + IssmDouble* bed = xNew(numnodes); + IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + element->GetInputListOnNodes(&thickness[0],ThicknessEnum); + element->GetInputListOnNodes(&bed[0],BaseEnum); + + /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/ + IssmDouble* head_old = xNew(numnodes); + element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum); + IssmDouble relaxation; + element->FindParam(&relaxation,HydrologyRelaxationEnum); + + /*Use the dof list to index into the solution vector: */ + for(int i=0;i hrho_ice*thickness[i]/rho_water+bed[i]){ + values[i] = rho_ice*thickness[i]/rho_water+bed[i]; + } + + /*Make sure that negative pressure is not allowed*/ + // if(values[i](values[i])) _error_("NaN found in solution vector"); + if(xIsInf(values[i])) _error_("Inf found in solution vector"); + } + + /*Add input to the element: */ + element->AddInput(HydrologyHeadEnum,values,element->GetElementType()); + element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum); + + /*Update reynolds number according to new solution*/ + element->GetVerticesCoordinates(&xyz_list); + Input* head_input = element->GetInput(HydrologyHeadEnum);_assert_(head_input); + head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list); + + /*Free resources:*/ + xDelete(values); + xDelete(thickness); + xDelete(bed); + xDelete(xyz_list); + xDelete(doflist); + xDelete(eff_pressure); + xDelete(head_old); +}/*}}}*/ +void HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 23936) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 23936) @@ -0,0 +1,33 @@ +/*! \file HydrologyGlaDSAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _HydrologyGlaDSAnalysis_ +#define _HydrologyGlaDSAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class HydrologyGlaDSAnalysis: public Analysis{ + + public: + /*Model processing*/ + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + void CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr=false); + 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 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/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 23935) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 23936) @@ -30,4 +30,5 @@ #include "./HydrologyDCInefficientAnalysis.h" #include "./HydrologyShreveAnalysis.h" +#include "./HydrologyGlaDSAnalysis.h" #include "./HydrologyShaktiAnalysis.h" #include "./HydrologyPismAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23935) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23936) @@ -823,4 +823,5 @@ if(ishydrology){ analyses_temp[numanalyses++]=HydrologyShreveAnalysisEnum; + analyses_temp[numanalyses++]=HydrologyGlaDSAnalysisEnum; analyses_temp[numanalyses++]=HydrologyShaktiAnalysisEnum; analyses_temp[numanalyses++]=HydrologyPismAnalysisEnum; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23935) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23936) @@ -424,4 +424,5 @@ iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); if(hydrology_model==HydrologydcEnum){ + /*FIXME: this cshould go to Analysis!!!*/ parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_porosity",HydrologydcSedimentPorosityEnum)); @@ -450,4 +451,7 @@ /*Nothing to add*/ } + else if(hydrology_model==HydrologyGlaDSEnum){ + /*Nothing to add*/ + } else{ _error_("Hydrology model "<=',0); + md = checkfield(md,'fieldname','hydrology.relaxation','>=',0); md = checkfield(md,'fieldname','hydrology.storage','>=',0); - md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1); + md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1); end % }}} function disp(self) % {{{ @@ -75,5 +75,5 @@ fielddisplay(self,'relaxation','under-relaxation coefficient for nonlinear iteration'); fielddisplay(self,'storage','englacial storage coefficient (void ratio)'); - fielddisplay(self,'requested_outputs','additional outputs requested'); + fielddisplay(self,'requested_outputs','additional outputs requested'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -91,13 +91,13 @@ WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','neumannflux','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double'); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','storage','format','Double'); - outputs = self.requested_outputs; - pos = find(ismember(outputs,'default')); - if ~isempty(pos), - outputs(pos) = []; %remove 'default' from outputs - outputs = [outputs defaultoutputs(self,md)]; %add defaults - end - WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); end % }}} end