Changeset 26747


Ignore:
Timestamp:
12/22/21 12:12:27 (3 years ago)
Author:
bulthuis
Message:

CGH:add spatially-varying tau and phi for sampling analysis

Location:
issm/trunk-jpl
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp

    r26047 r26747  
    7676  iomodel->FetchDataToInput(inputs,elements,"md.sampling.kappa",SamplingKappaEnum);
    7777        iomodel->FetchDataToInput(inputs,elements,"md.sampling.beta",SamplingBetaEnum,0.);
     78        iomodel->FetchDataToInput(inputs,elements,"md.sampling.tau",SamplingTauEnum);
    7879        iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleEnum,0.);
     80        iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleNoiseEnum,0.);
     81        if(iomodel->solution_enum==TransientSolutionEnum) iomodel->FetchDataToInput(inputs,elements,"md.sampling.phi",SamplingPhiEnum,0.);
    7982
    8083}/*}}}*/
    … …  
    8487        char**  requestedoutputs = NULL;
    8588
    86         parameters->AddObject(iomodel->CopyConstantObject("md.sampling.tau",SamplingTauEnum));
    8789        parameters->AddObject(iomodel->CopyConstantObject("md.sampling.alpha",SamplingAlphaEnum));
    8890        parameters->AddObject(iomodel->CopyConstantObject("md.sampling.robin",SamplingRobinEnum));
    89         if(solution_enum==TransientSolutionEnum) parameters->AddObject(iomodel->CopyConstantObject("md.sampling.phi",SamplingPhiEnum));
    9091        parameters->AddObject(iomodel->CopyConstantObject("md.sampling.seed",SamplingSeedEnum));
    9192
    … …  
    273274    element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
    274275    IssmDouble* newsample = xNew<IssmDouble>(numnodes);
     276                IssmDouble* tau = xNew<IssmDouble>(numnodes);
     277                element->GetInputListOnNodes(&tau[0],SamplingTauEnum);
    275278
    276279    /*Use the dof list to index into the solution vector: */
    277280    for(int i=0;i<numnodes;i++){
    278                         newsample[i]=solution[doflist[i]];
     281                        newsample[i]=solution[doflist[i]] / tau[i]; // new
     282
    279283                        /*Check solution*/
    280284                        if(xIsNan<IssmDouble>(newsample[i])) _error_("NaN found in solution vector");
    … …  
    287291    /*Free ressources:*/
    288292    xDelete<IssmDouble>(newsample);
     293                xDelete<IssmDouble>(tau);
    289294    xDelete<int>(doflist);
    290295 }/*}}}*/
    … …  
    418423        *pMff=Mff;
    419424}/*}}}*/
     425void                                            SamplingAnalysis::UpdateTransientSample(FemModel *      femmodel){
     426
     427        for(int j=0;j<femmodel->elements->Size();j++){
     428                        Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
     429                        UpdateTransientSample(element);
     430         }
     431
     432}
     433void                                            SamplingAnalysis::UpdateTransientSample(Element *       element){
     434
     435        /*Intermediaries */
     436        IssmDouble phi, sample, noise;
     437
     438        /*Fetch number vertices for this element*/
     439        int numvertices = element->GetNumberOfVertices();
     440
     441        /*Initialize new sample*/
     442        IssmDouble* sample_new = xNew<IssmDouble>(numvertices);
     443
     444        /*Retrieve all inputs and parameters*/
     445        Input*  sample_input=element->GetInput(SampleOldEnum); _assert_(sample_input);
     446        Input*  noise_input=element->GetInput(SampleNoiseEnum); _assert_(noise_input);
     447        Input*  phi_input=element->GetInput(SamplingPhiEnum); _assert_(phi_input);
     448
     449        /* Start  looping on the number of gaussian points: */
     450        Gauss* gauss=element->NewGauss();
     451  for(int iv=0;iv<numvertices;iv++){
     452        gauss->GaussVertex(iv);
     453
     454                /*Get input values at gauss points*/
     455                sample_input->GetInputValue(&sample,gauss);
     456                noise_input->GetInputValue(&noise,gauss);
     457                phi_input->GetInputValue(&phi,gauss);
     458
     459                /*Get new sample*/
     460                sample_new[iv] = phi*sample + noise;
     461
     462  }
     463
     464  element->AddInput(SampleEnum,sample_new,element->GetElementType());
     465
     466  /*Clean up and return*/
     467  xDelete<IssmDouble>(sample_new);
     468  delete gauss;
     469
     470}

  • issm/trunk-jpl/src/c/analyses/SamplingAnalysis.h

    r26047 r26747  
    4040                void           MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel);
    4141
     42                void                                    UpdateTransientSample(FemModel* femmodel);
     43                void                                    UpdateTransientSample(Element *         element);
     44
    4245};
    4346#endif

  • issm/trunk-jpl/src/c/cores/sampling_core.cpp

    r26004 r26747  
    55#include "../classes/classes.h"
    66#include "../solutionsequences/solutionsequences.h"
     7#include "../analyses/analyses.h" // new
     8#include "../modules/modules.h"
    79
    810void sampling_core(FemModel* femmodel){
    … …  
    3032        /*Generate random sample*/
    3133        if(VerboseSolution()) _printf0_("   call computational core\n");
     34        SamplingAnalysis* analysis = new SamplingAnalysis();
    3235        femmodel->SetCurrentConfiguration(SamplingAnalysisEnum);
     36
     37        if(solution_type==TransientSolutionEnum){
     38                InputDuplicatex(femmodel,SampleEnum,SampleOldEnum);
     39
     40    int seed;
     41                femmodel->parameters->FindParam(&seed,SamplingSeedEnum);
     42                if(seed>=0){
     43                        int step;
     44                        femmodel->parameters->FindParam(&step,StepEnum);
     45                        seed = seed + 13923272*step; // change default seed for transient simulations (by considering an arbitrary shift based on the step number)
     46                        femmodel->parameters->SetParam(seed,SamplingSeedEnum);
     47                }
     48
     49        }
     50
    3351        solutionsequence_sampling(femmodel);
     52
     53        if(solution_type==TransientSolutionEnum){
     54
     55                InputDuplicatex(femmodel,SampleEnum,SampleNoiseEnum);
     56
     57                analysis->UpdateTransientSample(femmodel);
     58
     59        }
     60
     61        delete analysis;
    3462
    3563        if(save_results){
    … …  
    4674        /*profiler*/
    4775        femmodel->profiler->Stop(SAMPLINGCORE);
     76
     77
    4878}

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r26556 r26747  
    150150        femmodel->parameters->FindParam(&issampling,TransientIssamplingEnum);
    151151        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
    152         femmodel->parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); 
     152        femmodel->parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
    153153
    154154        #if defined(_HAVE_OCEAN_ )

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r26676 r26747  
    352352syn keyword cConstant ModelnameEnum
    353353syn keyword cConstant SamplingAlphaEnum
    354 syn keyword cConstant SamplingPhiEnum
    355354syn keyword cConstant SamplingNumRequestedOutputsEnum
    356355syn keyword cConstant SamplingRequestedOutputsEnum
    357356syn keyword cConstant SamplingRobinEnum
    358357syn keyword cConstant SamplingSeedEnum
    359 syn keyword cConstant SamplingTauEnum
    360358syn keyword cConstant SaveResultsEnum
    361359syn keyword cConstant SolidearthPartitionIceEnum
    … …  
    803801syn keyword cConstant RheologyBbarAbsGradientEnum
    804802syn keyword cConstant SampleEnum
     803syn keyword cConstant SampleOldEnum
     804syn keyword cConstant SampleNoiseEnum
    805805syn keyword cConstant SamplingBetaEnum
    806806syn keyword cConstant SamplingKappaEnum
     807syn keyword cConstant SamplingPhiEnum
     808syn keyword cConstant SamplingTauEnum
    807809syn keyword cConstant SealevelEnum
    808810syn keyword cConstant SealevelGRDEnum
    … …  
    15801582syn keyword cType Cfsurfacesquare
    15811583syn keyword cType Channel
    1582 syn keyword cType classes
    15831584syn keyword cType Constraint
    15841585syn keyword cType Constraints
    … …  
    15871588syn keyword cType ControlInput
    15881589syn keyword cType Covertree
     1590syn keyword cType DataSetParam
    15891591syn keyword cType DatasetInput
    1590 syn keyword cType DataSetParam
    15911592syn keyword cType Definition
    15921593syn keyword cType DependentObject
    … …  
    16011602syn keyword cType ElementInput
    16021603syn keyword cType ElementMatrix
     1604syn keyword cType ElementVector
    16031605syn keyword cType Elements
    1604 syn keyword cType ElementVector
    16051606syn keyword cType ExponentialVariogram
    16061607syn keyword cType ExternalResult
    … …  
    16091610syn keyword cType Friction
    16101611syn keyword cType Gauss
    1611 syn keyword cType GaussianVariogram
    1612 syn keyword cType gaussobjects
    16131612syn keyword cType GaussPenta
    16141613syn keyword cType GaussSeg
    16151614syn keyword cType GaussTetra
    16161615syn keyword cType GaussTria
     1616syn keyword cType GaussianVariogram
    16171617syn keyword cType GenericExternalResult
    16181618syn keyword cType GenericOption
    … …  
    16301630syn keyword cType IssmDirectApplicInterface
    16311631syn keyword cType IssmParallelDirectApplicInterface
    1632 syn keyword cType krigingobjects
    16331632syn keyword cType Load
    16341633syn keyword cType Loads
    … …  
    16411640syn keyword cType Matice
    16421641syn keyword cType Matlitho
    1643 syn keyword cType matrixobjects
    16441642syn keyword cType MatrixParam
    16451643syn keyword cType Misfit
    … …  
    16541652syn keyword cType Observations
    16551653syn keyword cType Option
     1654syn keyword cType OptionUtilities
    16561655syn keyword cType Options
    1657 syn keyword cType OptionUtilities
    16581656syn keyword cType Param
    16591657syn keyword cType Parameters
    … …  
    16691667syn keyword cType Regionaloutput
    16701668syn keyword cType Results
     1669syn keyword cType RiftStruct
    16711670syn keyword cType Riftfront
    1672 syn keyword cType RiftStruct
    16731671syn keyword cType SealevelGeometry
    16741672syn keyword cType Seg
    16751673syn keyword cType SegInput
     1674syn keyword cType SegRef
    16761675syn keyword cType Segment
    1677 syn keyword cType SegRef
    16781676syn keyword cType SpcDynamic
    16791677syn keyword cType SpcStatic
    … …  
    16941692syn keyword cType Vertex
    16951693syn keyword cType Vertices
     1694syn keyword cType classes
     1695syn keyword cType gaussobjects
     1696syn keyword cType krigingobjects
     1697syn keyword cType matrixobjects
    16961698syn keyword cType AdjointBalancethickness2Analysis
    16971699syn keyword cType AdjointBalancethicknessAnalysis

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r26676 r26747  
    346346        ModelnameEnum,
    347347        SamplingAlphaEnum,
    348         SamplingPhiEnum,
    349348        SamplingNumRequestedOutputsEnum,
    350349        SamplingRequestedOutputsEnum,
    351350        SamplingRobinEnum,
    352351        SamplingSeedEnum,
    353         SamplingTauEnum,
    354352        SaveResultsEnum,
    355353        SolidearthPartitionIceEnum,
    … …  
    534532        TransientIsgroundinglineEnum,
    535533        TransientIshydrologyEnum,
    536         TransientIsmasstransportEnum, 
     534        TransientIsmasstransportEnum,
    537535        TransientIsoceantransportEnum,
    538536        TransientIsmovingfrontEnum,
    … …  
    799797        RheologyBbarAbsGradientEnum,
    800798        SampleEnum,
     799        SampleOldEnum,
     800        SampleNoiseEnum,
    801801        SamplingBetaEnum,
    802802        SamplingKappaEnum,
     803        SamplingPhiEnum,
     804        SamplingTauEnum,
    803805        SealevelEnum,
    804806        SealevelGRDEnum,
    … …  
    10361038        TemperatureSEMICEnum,
    10371039   ThermalforcingAutoregressionNoiseEnum,
    1038         ThermalforcingValuesAutoregressionEnum, 
     1040        ThermalforcingValuesAutoregressionEnum,
    10391041        ThermalSpctemperatureEnum,
    10401042        ThicknessAbsGradientEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r26676 r26747  
    354354                case ModelnameEnum : return "Modelname";
    355355                case SamplingAlphaEnum : return "SamplingAlpha";
    356                 case SamplingPhiEnum : return "SamplingPhi";
    357356                case SamplingNumRequestedOutputsEnum : return "SamplingNumRequestedOutputs";
    358357                case SamplingRequestedOutputsEnum : return "SamplingRequestedOutputs";
    359358                case SamplingRobinEnum : return "SamplingRobin";
    360359                case SamplingSeedEnum : return "SamplingSeed";
    361                 case SamplingTauEnum : return "SamplingTau";
    362360                case SaveResultsEnum : return "SaveResults";
    363361                case SolidearthPartitionIceEnum : return "SolidearthPartitionIce";
    … …  
    805803                case RheologyBbarAbsGradientEnum : return "RheologyBbarAbsGradient";
    806804                case SampleEnum : return "Sample";
     805                case SampleOldEnum : return "SampleOld";
     806                case SampleNoiseEnum : return "SampleNoise";
    807807                case SamplingBetaEnum : return "SamplingBeta";
    808808                case SamplingKappaEnum : return "SamplingKappa";
     809                case SamplingPhiEnum : return "SamplingPhi";
     810                case SamplingTauEnum : return "SamplingTau";
    809811                case SealevelEnum : return "Sealevel";
    810812                case SealevelGRDEnum : return "SealevelGRD";

  • issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

    r26676 r26747  
    345345syn keyword juliaConstC ModelnameEnum
    346346syn keyword juliaConstC SamplingAlphaEnum
    347 syn keyword juliaConstC SamplingPhiEnum
    348347syn keyword juliaConstC SamplingNumRequestedOutputsEnum
    349348syn keyword juliaConstC SamplingRequestedOutputsEnum
    350349syn keyword juliaConstC SamplingRobinEnum
    351350syn keyword juliaConstC SamplingSeedEnum
    352 syn keyword juliaConstC SamplingTauEnum
    353351syn keyword juliaConstC SaveResultsEnum
    354352syn keyword juliaConstC SolidearthPartitionIceEnum
    … …  
    796794syn keyword juliaConstC RheologyBbarAbsGradientEnum
    797795syn keyword juliaConstC SampleEnum
     796syn keyword juliaConstC SampleOldEnum
     797syn keyword juliaConstC SampleNoiseEnum
    798798syn keyword juliaConstC SamplingBetaEnum
    799799syn keyword juliaConstC SamplingKappaEnum
     800syn keyword juliaConstC SamplingPhiEnum
     801syn keyword juliaConstC SamplingTauEnum
    800802syn keyword juliaConstC SealevelEnum
    801803syn keyword juliaConstC SealevelGRDEnum

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r26676 r26747  
    360360              else if (strcmp(name,"Modelname")==0) return ModelnameEnum;
    361361              else if (strcmp(name,"SamplingAlpha")==0) return SamplingAlphaEnum;
    362               else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum;
    363362              else if (strcmp(name,"SamplingNumRequestedOutputs")==0) return SamplingNumRequestedOutputsEnum;
    364363              else if (strcmp(name,"SamplingRequestedOutputs")==0) return SamplingRequestedOutputsEnum;
    365364              else if (strcmp(name,"SamplingRobin")==0) return SamplingRobinEnum;
    366365              else if (strcmp(name,"SamplingSeed")==0) return SamplingSeedEnum;
    367               else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
    368366              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    369367              else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
    … …  
    383381              else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
    384382              else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum;
     383              else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
     384              else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
    389               else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
    390               else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum;
     388              if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum;
    391389              else if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum;
    392390              else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum;
    … …  
    506504              else if (strcmp(name,"Step")==0) return StepEnum;
    507505              else if (strcmp(name,"Steps")==0) return StepsEnum;
     506              else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
     507              else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
    512               else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
    513               else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
     511              if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
    514512              else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
    515513              else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
    … …  
    629627              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    630628              else if (strcmp(name,"BottomPressure")==0) return BottomPressureEnum;
     629              else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
     630              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
    635               else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    636               else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
     634              if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    637635              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    638636              else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
    … …  
    752750              else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
    753751              else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
     752              else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
     753              else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
    758               else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
    759               else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
     757              if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
    760758              else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
    761759              else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
    … …  
    823821              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
    824822              else if (strcmp(name,"Sample")==0) return SampleEnum;
     823              else if (strcmp(name,"SampleOld")==0) return SampleOldEnum;
     824              else if (strcmp(name,"SampleNoise")==0) return SampleNoiseEnum;
    825825              else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum;
    826826              else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum;
     827              else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum;
     828              else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
    827829              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    828830              else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
    … …  
    873875              else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
    874876              else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
    875               else if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum;
    876               else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum;
     880              if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum;
     881              else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum;
     882              else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum;
    881883              else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
    882884              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
    … …  
    996998              else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
    997999              else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum;
    998               else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
    999               else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
     1003              if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
     1004              else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
     1005              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
    10041006              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    10051007              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    … …  
    11191121              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    11201122              else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    1121               else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    1122               else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
     1126              if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
     1127              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
     1128              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    11271129              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    11281130              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    … …  
    12421244              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    12431245              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    1244               else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    1245               else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
     1249              if (strcmp(name,"Boundary")==0) return BoundaryEnum;
     1250              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
     1251              else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    12501252              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    12511253              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    … …  
    13651367              else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
    13661368              else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    1367               else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    1368               else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"ElementInput")==0) return ElementInputEnum;
     1372              if (strcmp(name,"Indexed")==0) return IndexedEnum;
     1373              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
     1374              else if (strcmp(name,"ElementInput")==0) return ElementInputEnum;
    13731375              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
    13741376              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
    … …  
    14881490              else if (strcmp(name,"Regular")==0) return RegularEnum;
    14891491              else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
    1490               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    1491               else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
     1495              if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
     1496              else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum;
     1497              else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
    14961498              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    14971499              else if (strcmp(name,"SMBautoregression")==0) return SMBautoregressionEnum;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp

    r26612 r26747  
    5151  Vector<IssmDouble>*  ug  = NULL;
    5252  Vector<IssmDouble>*  uf  = NULL;
    53   Vector<IssmDouble>*  old_uf  = NULL;     // previous solution vector (for transient)
    54   Vector<IssmDouble>*  pf  = NULL;
     53        Vector<IssmDouble>*  pf  = NULL;
    5554  Vector<IssmDouble>*  df  = NULL;
    5655  Vector<IssmDouble>*  ys=NULL;
    … …  
    6261
    6362  /*parameters:*/
    64   int solution_type, alpha, step, seed, nsize;
    65   IssmDouble phi, tau;
     63  int alpha, seed, nsize;
    6664
    6765  /*Recover parameters: */
    68   femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    6966  femmodel->parameters->FindParam(&seed,SamplingSeedEnum);
    7067  femmodel->parameters->FindParam(&alpha,SamplingAlphaEnum);
    71   femmodel->parameters->FindParam(&tau,SamplingTauEnum);
    72 
    73   /*Recover parameters for transient simulation: */
    74   if(solution_type==TransientSolutionEnum)
    75   {
    76     femmodel->parameters->FindParam(&step,StepEnum);
    77     femmodel->parameters->FindParam(&phi,SamplingPhiEnum);
    78                 if(seed>=0) seed = seed + 13923272*step; // change default seed for transient simulations (by considering an arbitrary shif based on the step number)
    79 
    80     GetSolutionFromInputsx(&ug,femmodel);
    81     Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
    82     old_uf=uf->Duplicate();
    83   }
    8468
    8569  /*CreateAnalysis*/
    … …  
    10185  Reduceloadx(pf, Kfs, ys);
    10286  delete Kfs;
    103 
    104   /*Copy old solution for transient run */
    105   if(solution_type==TransientSolutionEnum) uf->Copy(old_uf);
    10687
    10788  /* Generate random RHS */
    … …  
    130111  }
    131112
    132   /* Divide results by tau */
    133   uf->Scale(1.0/tau);
    134 
    135   /* Update solution x_{t+1} = phi x_{t} + noise for transient */
    136   if(solution_type==TransientSolutionEnum) uf->AXPY(old_uf,phi);
    137 
    138113  /* Update input */
    139114  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);
    … …  
    141116
    142117  /*clean-up*/
    143   delete Kff; delete pf; delete df; delete uf; delete ys; delete old_uf;
     118  delete Kff; delete pf; delete df; delete uf; delete ys;
    144119  delete Ml; delete Mscale;
    145120  delete analysis;

  • issm/trunk-jpl/src/m/classes/sampling.m

    r26301 r26747  
    99                tau               = 0;
    1010                beta              = NaN;
    11                 phi               = 0;
     11                phi               = NaN;
    1212                alpha             = 0;
    1313                robin             = 0;
    … …  
    3030                        disp(sprintf('\n      %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
    3131                        fielddisplay(self,'kappa','coefficient of the identity operator');
    32                         fielddisplay(self,'tau','scaling coefficient of the solution (default: 1.0)');
     32                        fielddisplay(self,'tau','scaling coefficient of the solution');
    3333                        fielddisplay(self,'alpha','exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)');
    3434
    … …  
    4747                function self = setdefaultparameters(self) % {{{
    4848
    49                         %Scaling coefficient
    50                         self.tau=1;
    51 
    5249                        %Apply Robin boundary conditions
    5350                        self.robin=0;
    54 
    55                         %Temporal correlation factor
    56                         self.phi=0;
    5751
    5852                        %Exponent in fraction SPDE (default: 2, biLaplacian covariance
    … …  
    7771
    7872                        md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
    79                         md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
     73                        md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
    8074                        md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
    8175                        if(md.sampling.robin)
    8276                                md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
    83                         end
    84                         md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
     77            end
    8578                        md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
    8679                        md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
    … …  
    9184
    9285                        WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
    93                         WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
     86                        WriteData(fid,prefix,'object',self,'fieldname','tau','format','DoubleMat','mattype',1);
    9487                        WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
    95                         WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
     88                        WriteData(fid,prefix,'object',self,'fieldname','phi','format','DoubleMat','mattype',1);
    9689                        WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
    9790                        WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
    … …  
    110103
    111104                        nu = self.alpha-1;
    112                         KAPPA = sqrt(8*nu)/lc;
    113                         TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
    114                         md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
    115                         md.sampling.tau = TAU;
     105                        KAPPA = sqrt(8*nu)./lc;
     106                        TAU = sqrt(gamma(nu)./(gamma(self.alpha)*(4*pi)*KAPPA.^(2*nu).*sigma.^2));
     107                        md.sampling.kappa = KAPPA.*ones(md.mesh.numberofvertices,1);
     108                        md.sampling.tau = TAU.*ones(md.mesh.numberofvertices,1);
    116109
    117110                end % }}}
    … …  
    121114                        writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
    122115                        writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
    123                         writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
     116                        writejsdouble(fid,[modelname '.sampling.phi'],self.phi);
    124117                        writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
    125118                        writejsdouble(fid,[modelname '.sampling.robin'],self.robin);

  • issm/trunk-jpl/test/NightlyRun/test134.m

    r26520 r26747  
    55md = md.sampling.setparameters(md,2e5,1);
    66md.sampling.seed = 100;
     7md.sampling.phi = zeros(md.mesh.numberofvertices,1);
    78md.cluster=generic('name',oshostname(),'np',1);
     9md.cluster.np=1;
    810md=solve(md,'smp');
    911
Note: See TracChangeset for help on using the changeset viewer.