Ignore:
Timestamp:
04/22/14 10:39:19 (11 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 17804

Location:
issm/trunk
Files:
18 edited
2 copied

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        1616.deps
        1717.dirstamp
         18.libs

  • issm/trunk/src/c/cores

    • Property svn:ignore

      •  

        old new  
         1*.obj
        12.deps
        23.dirstamp

  • issm/trunk/src/c/cores/AnalysisConfiguration.cpp

    r16533 r17806  
    2626
    2727                case StressbalanceSolutionEnum:
    28                         numanalyses=4;
     28                        numanalyses=6;
    2929                        analyses=xNew<int>(numanalyses);
    3030                        analyses[0]=StressbalanceAnalysisEnum;
    … …  
    3232                        analyses[2]=StressbalanceSIAAnalysisEnum;
    3333                        analyses[3]=L2ProjectionBaseAnalysisEnum;
     34                        analyses[4]=ExtrudeFromBaseAnalysisEnum;
     35                        analyses[5]=DepthAverageAnalysisEnum;
    3436                        break;
    3537
    … …  
    5557
    5658                case HydrologySolutionEnum:
    57                         numanalyses=4;
     59                        numanalyses=5;
    5860                        analyses=xNew<int>(numanalyses);
    5961                        analyses[0]=HydrologyShreveAnalysisEnum;
    … …  
    6163                        analyses[2]=HydrologyDCEfficientAnalysisEnum;
    6264                        analyses[3]=L2ProjectionBaseAnalysisEnum;
     65                        analyses[4]=L2ProjectionEPLAnalysisEnum;
    6366                        break;
    6467
    … …  
    114117
    115118                case TransientSolutionEnum:
    116                         numanalyses=12;
     119                        numanalyses=20;
    117120                        analyses=xNew<int>(numanalyses);
    118121                        analyses[ 0]=StressbalanceAnalysisEnum;
    … …  
    128131                        analyses[10]=ExtrudeFromBaseAnalysisEnum;
    129132                        analyses[11]=ExtrudeFromTopAnalysisEnum;
     133                        analyses[12]=LevelsetAnalysisEnum;
     134                        analyses[13]=ExtrapolationAnalysisEnum;
     135                        analyses[14]=LsfReinitializationAnalysisEnum;
     136                        analyses[15]=DamageEvolutionAnalysisEnum;
     137                        analyses[16]=HydrologyShreveAnalysisEnum;
     138                        analyses[17]=HydrologyDCInefficientAnalysisEnum;
     139                        analyses[18]=HydrologyDCEfficientAnalysisEnum;
     140                        analyses[19]=L2ProjectionEPLAnalysisEnum;
    130141                        break;
    131142

  • issm/trunk/src/c/cores/CorePointerFromSolutionEnum.cpp

    r16518 r17806  
    2424
    2525                case StressbalanceSolutionEnum:
    26                         #ifdef _HAVE_STRESSBALANCE_
    2726                        solutioncore=&stressbalance_core;
    28                         #else
    29                         _error_("ISSM was not compiled with stressbalance capabilities. Exiting");
    30                         #endif
    3127                        break;
    3228                case SteadystateSolutionEnum:
    33                         #ifdef _HAVE_STEADYSTATE_
    3429                        solutioncore=&steadystate_core;
    35                         #else
    36                         _error_("ISSM was not compiled with steady state capabilities. Exiting");
    37                         #endif
    3830                        break;
    3931                case ThermalSolutionEnum:
    40                         #ifdef _HAVE_THERMAL_
    4132                        solutioncore=&thermal_core;
    42                         #else
    43                         _error_("ISSM was not compiled with thermal capabilities. Exiting");
    44                         #endif
    4533                        break;
    4634                case BalancethicknessSolutionEnum:
    47                         #ifdef _HAVE_BALANCED_
    4835                        solutioncore=&balancethickness_core;
    49                         #else
    50                         _error_("ISSM was not compiled with balanced capabilities. Exiting");
    51                         #endif
    5236                        break;
    5337                case BalancethicknessSoftSolutionEnum:
    54                         #ifdef _HAVE_BALANCED_
    5538                        solutioncore=&dummy_core;
    56                         #else
    57                         _error_("ISSM was not compiled with balanced capabilities. Exiting");
    58                         #endif
    5939                        break;
    6040                case BalancevelocitySolutionEnum:
    61                         #ifdef _HAVE_BALANCED_
    6241                        solutioncore=&balancevelocity_core;
    63                         #else
    64                         _error_("ISSM was not compiled with balanced capabilities. Exiting");
    65                         #endif
    6642                        break;
    6743                case HydrologySolutionEnum:
    68                         #ifdef _HAVE_HYDROLOGY_
    6944                        solutioncore=&hydrology_core;
    70                         #else
    71                         _error_("ISSM was not compiled with hydrology capabilities. Exiting");
    72                         #endif
    7345                        break;
    7446                case SurfaceSlopeSolutionEnum:
    75                         #ifdef _HAVE_SLOPE_
    7647                        solutioncore=&surfaceslope_core;
    77                         #else
    78                         _error_("ISSM was not compiled with slope capabilities. Exiting");
    79                         #endif
    8048                        break;
    8149                case BedSlopeSolutionEnum:
    82                         #ifdef _HAVE_SLOPE_
    8350                        solutioncore=&bedslope_core;
    84                         #else
    85                         _error_("ISSM was not compiled with slope capabilities. Exiting");
    86                         #endif
    8751                        break;
    8852                case TransientSolutionEnum:
    89                         #ifdef _HAVE_TRANSIENT_
    9053                        solutioncore=&transient_core;
    91                         #else
    92                         _error_("ISSM was not compiled with transient capabilities. Exiting");
    93                         #endif
    9454                        break;
    9555                case MasstransportSolutionEnum:
    96                         #ifdef _HAVE_MASSTRANSPORT_
    9756                        solutioncore=&masstransport_core;
    98                         #else
    99                         _error_("ISSM was not compiled with masstransport capabilities. Exiting");
    100                         #endif
    10157                        break;
    102 
    10358                case GiaSolutionEnum:
    104                         #ifdef _HAVE_GIA_
     59                        #if _HAVE_GIA_
    10560                        solutioncore=&gia_core;
    10661                        #else
    107                         _error_("ISSM was not compiled with gia capabilities. Exiting");
     62                        _error_("ISSM not compiled with Gia capability");
    10863                        #endif
    10964                        break;
    11065                case DamageEvolutionSolutionEnum:
    111                         #ifdef _HAVE_DAMAGE_
    11266                        solutioncore=&damage_core;
    113                         #else
    114                         _error_("ISSM was not compiled with damage evolution capabilities. Exiting");
    115                         #endif
    11667                        break;
    117 
    11868                default:
    11969                        _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");
    … …  
    12474        _assert_(psolutioncore);
    12575        *psolutioncore=solutioncore;
    126 
    12776}

  • issm/trunk/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp

    r16518 r17806  
    4242        }
    4343        else if(control_analysis){
    44                 #ifdef _HAVE_CONTROL_
    4544                if(tao_analysis) solutioncore=controltao_core;
    4645                else solutioncore=control_core;
    47                 #else
    48                 _error_("ISSM was not compiled with control support, cannot carry out control analysis!");
    49                 #endif
    5046        }
    5147        else CorePointerFromSolutionEnum(&solutioncore,parameters,solutiontype);  /*This means we retrieve a core solution that is not a wrapper*/

  • issm/trunk/src/c/cores/bedslope_core.cpp

    r16529 r17806  
    1414        /*parameters: */
    1515        bool save_results;
    16         int  meshtype;
     16        int  domaintype;
    1717
    1818        /*Recover some parameters: */
    1919        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    20         femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
     20        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    2121
    2222        if(VerboseSolution()) _printf0_("   computing slope\n");
    … …  
    2828        solutionsequence_linear(femmodel);
    2929
    30         if(meshtype!=Mesh2DverticalEnum){
     30        if(domaintype!=Domain2DverticalEnum){
    3131                femmodel->parameters->SetParam(BedSlopeYEnum,InputToL2ProjectEnum);
    3232                solutionsequence_linear(femmodel);
    … …  
    3535        if(save_results){
    3636                if(VerboseSolution()) _printf0_("   saving results\n");
    37                 if(meshtype!=Mesh2DverticalEnum){
     37                if(domaintype!=Domain2DverticalEnum){
    3838                        int outputs[2] = {BedSlopeXEnum,BedSlopeYEnum};
    3939                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);

  • issm/trunk/src/c/cores/cores.h

    r16534 r17806  
    2424void thermal_core(FemModel* femmodel);
    2525void surfaceslope_core(FemModel* femmodel);
     26void levelsetfunctionslope_core(FemModel* femmodel);
    2627void bedslope_core(FemModel* femmodel);
    2728void meshdeformation_core(FemModel* femmodel);
    … …  
    2930void controltao_core(FemModel* femmodel);
    3031void masstransport_core(FemModel* femmodel);
     32void depthaverage_core(FemModel* femmodel);
    3133void extrudefrombase_core(FemModel* femmodel);
    3234void extrudefromtop_core(FemModel* femmodel);

  • issm/trunk/src/c/cores/damage_core.cpp

    r16518 r17806  
    1414        /*intermediary*/
    1515        bool   save_results;
    16         bool   dakota_analysis  = false;
     16        bool   dakota_analysis     = false;
    1717        int    solution_type;
     18        int    numoutputs          = 0;
     19        char   **requested_outputs = NULL;
    1820
    19         if(VerboseSolution()) _printf0_("   computing damage\n");
    2021       
    2122        //first recover parameters common to all solutions
    … …  
    2324        femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
    2425        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     26        femmodel->parameters->FindParam(&numoutputs,DamageEvolutionNumRequestedOutputsEnum);
     27        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,DamageEvolutionRequestedOutputsEnum);
    2528
    2629        if(dakota_analysis && solution_type!=TransientSolutionEnum){
    … …  
    2932        }
    3033
     34        if(VerboseSolution()) _printf0_("   computing damage\n");
    3135        femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
    3236        solutionsequence_linear(femmodel);
    … …  
    3438        if(save_results){
    3539                if(VerboseSolution()) _printf0_("   saving results\n");
    36                 int outputs = DamageDEnum;
    37                 femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
     40                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
     41        }
     42       
     43        /*Free resources:*/     
     44        if(numoutputs){
     45                for(int i=0;i<numoutputs;i++){
     46                        xDelete<char>(requested_outputs[i]);
     47                }
     48                xDelete<char*>(requested_outputs);
    3849        }
    3950}

  • issm/trunk/src/c/cores/gia_core.cpp

    r16518 r17806  
    5050        if(save_results){
    5151                if(VerboseSolution()) _printf0_("   saving results\n");
    52                 const int outputs[2] = {GiaWEnum,GiadWdtEnum};
     52                int outputs[2] = {GiaWEnum,GiadWdtEnum};
    5353                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
    5454        }

  • issm/trunk/src/c/cores/hydrology_core.cpp

    r16518 r17806  
    1212void hydrology_core(FemModel* femmodel){
    1313
    14         int i;
    15 
    1614        /*intermediary*/
    17         int        step,nsteps;
    18         int        output_frequency,hydrology_model;
     15        int        hydrology_model;
    1916        bool       save_results;
    2017        bool       modify_loads=true;
    2118        bool       isefficientlayer;
    22         IssmDouble starttime,final_time;
    23         IssmDouble time,dt;
    2419
    2520        /*first recover parameters common to all solutions*/
    26         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    27         femmodel->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
    28         femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    2921        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    30         femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
    3122        femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
    3223
    … …  
    3728        }
    3829
    39         /*Compute number of time steps: */
    40         if((dt==0)|| (final_time==0)){
    41                 dt=0;
    42                 nsteps=1;
     30        /*Using the Shreve based Model*/
     31        if (hydrology_model==HydrologyshreveEnum){
     32                if(VerboseSolution()) _printf0_("   computing water column\n");
     33                femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
     34                solutionsequence_nonlinear(femmodel,modify_loads);
     35               
     36                /*transfer water column thickness to old water column thickness: */
     37               
     38                InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
     39               
     40                if(save_results){
     41                        if(VerboseSolution()) _printf0_("   saving results \n");
     42                        int outputs[3] = {WatercolumnEnum,HydrologyWaterVxEnum,HydrologyWaterVyEnum};
     43                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
     44                       
     45                        /*unload results*/
     46                        if(VerboseSolution()) _printf0_("   saving temporary results\n");
     47                        OutputResultsx(femmodel);
     48                }
    4349        }
    44         else nsteps=reCast<int,IssmDouble>((final_time-starttime)/dt);
    4550
    46         /*initialize: */
    47         step=0;
    48         time=starttime;
    49 
    50         /*Loop through time: */
    51         for(i=0;i<nsteps;i++){
    52 
    53                 if(nsteps)if(VerboseSolution()) _printf0_("time step:" << i+1 << "/" << nsteps << "\n");
    54                 time+=dt;
    55                 step+=1;
    56                 femmodel->parameters->SetParam(time,TimeEnum);
    57                 femmodel->parameters->SetParam(step,StepEnum);
    58 
    59                 if (hydrology_model==HydrologyshreveEnum){
    60                         if(VerboseSolution()) _printf0_("   computing water column\n");
    61                         femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
    62                         solutionsequence_nonlinear(femmodel,modify_loads);
    63 
    64                         /*transfer water column thickness to old water column thickness: */
    65 
    66                         InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
    67 
    68                         if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
    69                                 if(VerboseSolution()) _printf0_("   saving results \n");
    70                                 int outputs[3] = {WatercolumnEnum,HydrologyWaterVxEnum,HydrologyWaterVyEnum};
     51        /*Using the double continuum model*/
     52        else if (hydrology_model==HydrologydcEnum){
     53                InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
     54                femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     55                if (isefficientlayer){
     56                        InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
     57                        InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
     58                }
     59               
     60                /*Proceed now to heads computations*/
     61                if(VerboseSolution()) _printf0_("   computing water head\n");
     62                solutionsequence_hydro_nonlinear(femmodel);
     63                if(save_results){
     64                        if(VerboseSolution()) _printf0_("   saving results \n");
     65                        if(isefficientlayer){
     66                                int outputs[9] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,HydrologydcMaskEplactiveEltEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum,EffectivePressureEnum};
     67                                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],9);
     68                        }
     69                        else{
     70                                int outputs[3] = {SedimentHeadEnum,SedimentHeadResidualEnum,EffectivePressureEnum};
    7171                                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
    72 
    73                                 /*unload results*/
    74                                 if(VerboseSolution()) _printf0_("   saving temporary results\n");
    75                                 OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
    7672                        }
     73                        /*unload results*/
     74                        if(VerboseSolution()) _printf0_("   saving temporary results\n");
     75                        OutputResultsx(femmodel);
    7776                }
    78 
    79                 else if (hydrology_model==HydrologydcEnum){
    80                         InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
    81                         femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    82                         if (isefficientlayer){
    83                                 InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
    84                         }
    85 
    86                         if(VerboseSolution()) _printf0_("   computing water head\n");
    87                         solutionsequence_hydro_nonlinear(femmodel);
    88                         if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
    89                                 if(VerboseSolution()) _printf0_("   saving results \n");
    90                                 if(isefficientlayer){
    91                                         int outputs[4] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum};
    92                                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
    93                                 }
    94                                 else{
    95                                         int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum};
    96                                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
    97                                 }
    98                                 /*unload results*/
    99                                 if(VerboseSolution()) _printf0_("   saving temporary results\n");
    100                                 OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
    101                         }
    102                 }
    103 
    10477        }
    10578}
     79

  • issm/trunk/src/c/cores/masstransport_core.cpp

    r16518 r17806  
    1414        /*parameters: */
    1515        int    i;
    16         int    numoutputs,meshtype;
     16        int    numoutputs,domaintype;
    1717        bool   save_results;
    1818        bool   issmbgradients,ispdd,isdelta18o,isFS,isfreesurface,dakota_analysis;
    … …  
    2525        /*recover parameters: */
    2626        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    27         femmodel->parameters->FindParam(&issmbgradients,SurfaceforcingsIssmbgradientsEnum);
    28         femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum);
    29         femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
    3027        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    3128        femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum);
    3229        femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
    3330        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    34         femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
     31        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    3532        femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
    3633        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
    … …  
    4037                InputDuplicatex(femmodel,QmuSurfaceEnum,SurfaceEnum);
    4138                InputDuplicatex(femmodel,QmuThicknessEnum,ThicknessEnum);
    42                 InputDuplicatex(femmodel,QmuBedEnum,BedEnum);
     39                InputDuplicatex(femmodel,QmuBaseEnum,BaseEnum);
    4340        }
    4441
    45         if(issmbgradients){
    46           if(VerboseSolution())_printf_("       call smb gradients module\n\n");
    47           SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    48         }
    49         if(ispdd){
    50                 if(isdelta18o){
    51                         if(VerboseSolution()) _printf0_("   call Delta18oParametrization module\n");
    52                         Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    53                 }
    54                 if(VerboseSolution()) _printf0_("   call positive degree day module\n");
    55                 PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    56         }
     42        /*Calculate new Surface Mass Balance (SMB)*/
     43        SurfaceMassBalancex(femmodel);
    5744
     45        /*Transport mass or free surface*/
    5846        if(isFS && isfreesurface){
    5947                if(VerboseSolution()) _printf0_("   call free surface computational core\n");
    6048                femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum);
    6149                solutionsequence_linear(femmodel);
    62                 if(meshtype==Mesh2DverticalEnum){
    63                         femmodel->parameters->SetParam(BedEnum,InputToExtrudeEnum);
     50                if(domaintype!=Domain2DhorizontalEnum){
     51                        femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
    6452                        extrudefrombase_core(femmodel);
    6553                }
    6654                femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum);
    6755                solutionsequence_linear(femmodel);
    68                 if(meshtype==Mesh2DverticalEnum){
     56                if(domaintype!=Domain2DhorizontalEnum){
    6957                        femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
    7058                        extrudefromtop_core(femmodel);
    … …  
    7462                if(VerboseSolution()) _printf0_("   call computational core\n");
    7563                solutionsequence_linear(femmodel);
     64                if(domaintype==Domain2DverticalEnum){
     65                        femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
     66                        extrudefrombase_core(femmodel);
     67                        femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
     68                        extrudefrombase_core(femmodel);
     69                        femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
     70                        extrudefrombase_core(femmodel);
     71                }
    7672        }
    7773
    … …  
    8480
    8581        /*Free ressources:*/
    86         if(numoutputs){for(int i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);}
     82        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    8783}

  • issm/trunk/src/c/cores/steadystate_core.cpp

    r16518 r17806  
    5050
    5151                if(VerboseSolution()) _printf0_("   computing temperature and velocity for step: " << step << "\n");
    52                 #ifdef _HAVE_THERMAL_
    5352                thermal_core(femmodel);
    5453                if(!isenthalpy)femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);/*Could be MeltingAnalysis...*/
    5554                GetSolutionFromInputsx(&tg,femmodel);
    56                 #else
    57                 _error_("ISSM was not compiled with thermal capabilities. Exiting");
    58                 #endif
    5955
    6056                if(VerboseSolution()) _printf0_("   computing new velocity\n");
    … …  
    8783        delete tg;
    8884        delete ug;     
    89         if(numoutputs){ for (i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);}
     85        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    9086}
    9187bool steadystateconvergence(Vector<IssmDouble>* tg,Vector<IssmDouble>* tg_old,Vector<IssmDouble>* ug,Vector<IssmDouble>* ug_old,IssmDouble reltol){

  • issm/trunk/src/c/cores/stressbalance_core.cpp

    r16518 r17806  
    66#include "../toolkits/toolkits.h"
    77#include "../classes/classes.h"
     8#include "../analyses/analyses.h"
    89#include "../shared/shared.h"
    910#include "../modules/modules.h"
    … …  
    1314
    1415        /*parameters: */
    15         bool  dakota_analysis;
    16         int   meshtype;
    17         bool  isSIA,isSSA,isL1L2,isHO,isFS;
    18         bool  conserve_loads    = true;
    19         bool  save_results;
    20         int   newton;
    21         int   solution_type;
    22         int   numoutputs        = 0;
    23         char** requested_outputs = NULL;
    24         int    i;
    25 
     16        bool       dakota_analysis;
     17        int        domaintype;
     18        bool       isSIA,isSSA,isL1L2,isHO,isFS;
     19        bool       save_results;
     20        int        solution_type;
     21        int        numoutputs        = 0;
     22        char     **requested_outputs = NULL;
     23        Analysis  *analysis          = NULL;
    2624
    2725        /* recover parameters:*/
    28         femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
     26        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    2927        femmodel->parameters->FindParam(&isSIA,FlowequationIsSIAEnum);
    3028        femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
    … …  
    3230        femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
    3331        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    34         femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
    3532        femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
    3633        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    … …  
    4340                InputDuplicatex(femmodel,QmuVxEnum,VxEnum);
    4441                InputDuplicatex(femmodel,QmuVyEnum,VyEnum);
    45                 InputDuplicatex(femmodel,QmuVzEnum,VzEnum);
     42                if(domaintype==Domain3DEnum) InputDuplicatex(femmodel,QmuVzEnum,VzEnum);
    4643                if(isFS) InputDuplicatex(femmodel,QmuPressureEnum,PressureEnum);
    4744        }
    4845
    49         /*Compute slopes: */
    50         if(isSIA) surfaceslope_core(femmodel);
     46        /*Compute slopes if necessary */
     47        if(isSIA || (isFS && domaintype==Domain2DverticalEnum)) surfaceslope_core(femmodel);
    5148        if(isFS){
    5249                bedslope_core(femmodel);
    5350                femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
    54                 ResetCoordinateSystemx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     51                ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    5552        }
    5653
     54        /*Compute SIA velocities*/
    5755        if(isSIA){
    58                 if(VerboseSolution()) _printf0_("   computing SIA velocities\n");
    5956
    6057                /*Take the last velocity into account so that the velocity on the SSA domain is not zero*/
    6158                if(isSSA || isL1L2 || isHO ) ResetBoundaryConditions(femmodel,StressbalanceSIAAnalysisEnum);
    62                 femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum);
    63                 solutionsequence_linear(femmodel);
     59
     60                analysis = new StressbalanceSIAAnalysis();
     61                analysis->Core(femmodel);
     62                delete analysis;
     63
     64                /*Reset velocities for other ice flow models*/
    6465                if(isSSA || isL1L2 || isHO) ResetBoundaryConditions(femmodel,StressbalanceAnalysisEnum);
    6566        }
    6667
    67         if ((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor
    68                 if(VerboseSolution()) _printf0_("   computing velocities\n");
    69 
    70                 femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
    71                 if(newton>0)
    72                  solutionsequence_newton(femmodel);
    73                 else
    74                  solutionsequence_nonlinear(femmodel,conserve_loads);
     68        /*Compute stressbalance for SSA L1L2 HO and FS*/
     69        if(isSSA || isL1L2 || isHO || isFS){
     70                analysis = new StressbalanceAnalysis();
     71                analysis->Core(femmodel);
     72                delete analysis;
    7573        }
    7674
    77         if ((isSSA || isL1L2 || isHO) && isFS){
    78                 if(VerboseSolution()) _printf0_("   computing coupling betweem lower order models and full-FS\n");
    79                 solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads);
     75        /*Compute vertical velocities*/
     76        if (domaintype==Domain3DEnum && (isSIA || isSSA || isL1L2 || isHO)){
     77                analysis = new StressbalanceVerticalAnalysis();
     78                analysis->Core(femmodel);
     79                delete analysis;
    8080        }
    8181
    82         if (meshtype==Mesh3DEnum && (isSIA || isSSA || isL1L2 || isHO)){
    83                 if(VerboseSolution()) _printf0_("   computing vertical velocities\n");
    84                 femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
    85                 solutionsequence_linear(femmodel);
    86         }
    8782
    8883        if(save_results){
    … …  
    9489
    9590        /*Free ressources:*/   
    96         if(numoutputs){ for (i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);}
     91        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    9792}

  • issm/trunk/src/c/cores/surfaceslope_core.cpp

    r16529 r17806  
    1414        /*parameters: */
    1515        bool save_results;
    16         int  meshtype;
     16        int  domaintype;
    1717
    1818        /*Recover some parameters: */
    1919        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    20         femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
     20        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    2121
    2222        if(VerboseSolution()) _printf0_("computing slope...\n");
    … …  
    2828        solutionsequence_linear(femmodel);
    2929
    30         if(meshtype!=Mesh2DverticalEnum){
     30        if(domaintype!=Domain2DverticalEnum){
    3131                femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToL2ProjectEnum);
    3232                solutionsequence_linear(femmodel);
     33        }
     34        if(domaintype==Domain2DverticalEnum){
     35                femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum);
     36                extrudefrombase_core(femmodel);
    3337        }
    3438
    3539        if(save_results){
    3640                if(VerboseSolution()) _printf0_("saving results:\n");
    37                 if(meshtype!=Mesh2DverticalEnum){
     41                if(domaintype!=Domain2DverticalEnum){
    3842                        int outputs[2] = {SurfaceSlopeXEnum,SurfaceSlopeYEnum};
    3943                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);

  • issm/trunk/src/c/cores/thermal_core.cpp

    r16518 r17806  
    66#include "../toolkits/toolkits.h"
    77#include "../classes/classes.h"
     8#include "../analyses/analyses.h"
    89#include "../shared/shared.h"
    910#include "../modules/modules.h"
    … …  
    1718        int    solution_type,numoutputs;
    1819        char** requested_outputs = NULL;
     20        EnthalpyAnalysis * enthalpy_analysis = NULL;
    1921
    2022        /*first recover parameters common to all solutions*/
    … …  
    4648                InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum);
    4749
    48                 /*Post process*/
    49                 if(solution_type!=SteadystateSolutionEnum) PostprocessingEnthalpyx(femmodel);
     50                if(solution_type!=SteadystateSolutionEnum){
     51                        /*Post process*/
     52                        enthalpy_analysis = new EnthalpyAnalysis();
     53                        enthalpy_analysis->PostProcessing(femmodel);
     54                        delete enthalpy_analysis;
     55                }
    5056        }
    5157        else{
    … …  
    6369                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    6470        }
     71
     72        /*Free ressources:*/   
     73        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
     74
    6575}

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

    r16518 r17806  
    2222        int    i;
    2323        IssmDouble starttime,finaltime,dt,yts;
    24         bool   isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isdelta18o,isgia;
     24        bool   isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology;
    2525        bool   save_results,dakota_analysis;
    2626        bool   time_adapt=false;
    2727        int    output_frequency;
    28         int    meshtype,groundingline_migration;
     28        int    domaintype,groundingline_migration,smb_model;
    2929        int    numoutputs         = 0;
     30  Analysis *analysis = NULL;
    3031        char** requested_outputs = NULL;
    3132
    … …  
    3637
    3738        //first recover parameters common to all solutions
    38         femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
     39        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    3940        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    4041        femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
    … …  
    4950        femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
    5051        femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
     52        femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
     53        femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
     54        femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
    5155        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    5256        if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
    5357        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
    5458        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
    55         femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
    5659
    5760        /*initialize: */
    … …  
    6467                        InputDuplicatex(femmodel,QmuVxEnum,VxEnum);
    6568                        InputDuplicatex(femmodel,QmuVyEnum,VyEnum);
    66                         if(meshtype==Mesh3DEnum){
     69                        if(domaintype==Domain3DEnum){
    6770                                InputDuplicatex(femmodel,QmuVzEnum,VzEnum);
    6871                                if(isFS)InputDuplicatex(femmodel,QmuPressureEnum,PressureEnum);
    … …  
    7275                        InputDuplicatex(femmodel,QmuThicknessEnum,ThicknessEnum);
    7376                        InputDuplicatex(femmodel,QmuSurfaceEnum,SurfaceEnum);
    74                         InputDuplicatex(femmodel,QmuBedEnum,BedEnum);
     77                        InputDuplicatex(femmodel,QmuBaseEnum,BaseEnum);
    7578                        InputDuplicatex(femmodel,QmuMaskIceLevelsetEnum,MaskIceLevelsetEnum);
    7679                }
    7780                if(isgroundingline) InputDuplicatex(femmodel,QmuMaskGroundediceLevelsetEnum,MaskGroundediceLevelsetEnum);
    78                 if(isthermal && meshtype==Mesh3DEnum){
     81                if(isthermal && domaintype==Domain3DEnum){
    7982                        //Update Vertex Position after updating Thickness and Bed
    8083                        femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum);
    … …  
    9396
    9497        while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.
    95 
    9698                /*Increment*/
    9799                if(time_adapt){
    … …  
    106108
    107109                if(VerboseSolution()) _printf0_("iteration " << step << "/" << floor((finaltime-time)/dt)+step << "  time [yr]: " << time/yts << " (time step: " << dt/yts << ")\n");
    108                 if(step%output_frequency==0 || time==finaltime)
     110                if(step%output_frequency==0 || (time >= finaltime - (yts*DBL_EPSILON)) || step==1)
    109111                 save_results=true;
    110112                else
    … …  
    112114                femmodel->parameters->SetParam(save_results,SaveResultsEnum);
    113115
    114                 if(isthermal && meshtype==Mesh3DEnum){
    115                         if(VerboseSolution()) _printf0_("   computing temperatures\n");
    116                         #ifdef _HAVE_THERMAL_
     116                if(isthermal && domaintype==Domain3DEnum){
     117                        if(VerboseSolution()) _printf0_("   computing thermal regime\n");
    117118                        thermal_core(femmodel);
    118                         #else
    119                         _error_("ISSM was not compiled with thermal capabilities. Exiting");
    120                         #endif
     119                }
     120
     121                if(ishydrology){
     122                        if(VerboseSolution()) _printf0_("   computing water heads\n");
     123                        hydrology_core(femmodel);
    121124                }
    122125
    123126                if(isstressbalance){
    124127                        if(VerboseSolution()) _printf0_("   computing new velocity\n");
    125                         #ifdef _HAVE_STRESSBALANCE_
    126128                        stressbalance_core(femmodel);
    127                         #else
    128                         _error_("ISSM was not compiled with stressbalance capabilities. Exiting");
    129                         #endif
     129                }
     130
     131                if(isdamageevolution){
     132                        if(VerboseSolution()) _printf0_("   computing damage\n");
     133                        damage_core(femmodel);
     134                }
     135
     136                if(islevelset){
     137                        if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
     138                        /* smoothen slope of lsf for computation of normal on ice domain*/
     139                        levelsetfunctionslope_core(femmodel);
     140
     141                        /* extrapolate velocities onto domain with no ice */
     142                        Analysis* extanalysis = new ExtrapolationAnalysis();
     143                        const int nvars=2;
     144                        int vars[nvars] = {VxEnum, VyEnum};
     145                        for(int iv=0;iv<nvars;iv++){
     146                                femmodel->parameters->SetParam(vars[iv],ExtrapolationVariableEnum);
     147                                extanalysis->Core(femmodel);
     148                        }
     149                        delete extanalysis;     
     150
     151                        /* solve level set equation */
     152                        analysis = new LevelsetAnalysis();
     153                        analysis->Core(femmodel);
     154                        delete analysis;
     155
     156                        /* update vertices included for next calculation */
     157                        GetMaskOfIceVerticesLSMx(femmodel);
     158
     159                        /* add computation domain mask to outputs */
     160                        int outputs[1] = {IceMaskNodeActivationEnum};
     161                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
    130162                }
    131163
    … …  
    136168                        femmodel->UpdateVertexPositionsx();
    137169                }
    138 
     170               
    139171                if(isgroundingline){
    140172                        if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
    141                         #ifdef _HAVE_GROUNDINGLINE_
    142173                        GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    143                         #else
    144                         _error_("ISSM was not compiled with grounding line migration capabilities. Exiting");
    145                         #endif
     174
     175                        if(groundingline_migration==ContactEnum){
     176                                femmodel->parameters->SetParam(MaskGroundediceLevelsetEnum,InputToExtrudeEnum);
     177                                extrudefrombase_core(femmodel);
     178                                femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
     179                                extrudefrombase_core(femmodel);
     180                        }
    146181                        if(save_results){
    147                                 int outputs[3] = {SurfaceEnum,BedEnum,MaskGroundediceLevelsetEnum};
     182                                int outputs[3] = {SurfaceEnum,BaseEnum,MaskGroundediceLevelsetEnum};
    148183                                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
    149184                        }
    … …  
    156191                        _error_("ISSM was not compiled with gia capabilities. Exiting");
    157192                        #endif
    158 
    159193                }
    160194
    161195                /*unload results*/
     196                if(VerboseSolution()) _printf0_("   computing requested outputs\n");
     197                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs,save_results);
     198                if(isgroundingline && (groundingline_migration==SubelementMigrationEnum || groundingline_migration==SubelementMigration2Enum)){
     199                        int outputs[1] = {MaskGroundediceLevelsetEnum};
     200                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1,save_results);
     201                }
     202
    162203                if(save_results){
    163                         if(VerboseSolution()) _printf0_("   saving transient results\n");
    164                         femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    165                         if(isdelta18o){
    166                                 int outputs[2] = {SurfaceforcingsMonthlytemperaturesEnum,SurfaceforcingsPrecipitationEnum};
    167                                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
    168                         }
    169                         if(isgroundingline && (groundingline_migration==SubelementMigrationEnum || groundingline_migration==SubelementMigration2Enum)){
    170                                 int outputs[1] = {MaskGroundediceLevelsetEnum};
    171                                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
    172                         }
    173204                        if(VerboseSolution()) _printf0_("   saving temporary results\n");
    174                         OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
     205                        OutputResultsx(femmodel);
    175206                }
    176207        }
    … …  
    179210
    180211        /*Free ressources:*/   
    181         if(numoutputs){ for (i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);}
     212        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    182213}
Note: See TracChangeset for help on using the changeset viewer.