source: issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp@ 16604

#include "./EnthalpyAnalysis.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

/*Model processing*/
int  EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
        return 1;
}/*}}}*/
void EnthalpyAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/

        int     numoutputs;
        char**  requestedoutputs = NULL;

        parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));

        iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
        parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
        if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
        iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum);
}/*}}}*/
void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

        bool dakota_analysis;
        bool isenthalpy;

        /*Now, is the model 3d? otherwise, do nothing: */
        if(iomodel->meshtype==Mesh2DhorizontalEnum)return;

        /*Is enthalpy requested?*/
        iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
        if(!isenthalpy) return;

        /*Fetch data needed: */
        iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);

        /*Update elements: */
        int counter=0;
        for(int i=0;i<iomodel->numberofelements;i++){
                if(iomodel->my_elements[i]){
                        Element* element=(Element*)elements->GetObjectByOffset(counter);
                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
                        counter++;
                }
        }

        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);

        iomodel->FetchDataToInput(elements,ThicknessEnum);
        iomodel->FetchDataToInput(elements,SurfaceEnum);
        iomodel->FetchDataToInput(elements,BedEnum);
        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
        iomodel->FetchDataToInput(elements,FrictionPEnum);
        iomodel->FetchDataToInput(elements,FrictionQEnum);
        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
        iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
        iomodel->FetchDataToInput(elements,PressureEnum);
        iomodel->FetchDataToInput(elements,TemperatureEnum);
        iomodel->FetchDataToInput(elements,WaterfractionEnum);
        iomodel->FetchDataToInput(elements,EnthalpyEnum);
        iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
        iomodel->FetchDataToInput(elements,WatercolumnEnum);
        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
        iomodel->FetchDataToInput(elements,VxEnum);
        iomodel->FetchDataToInput(elements,VyEnum);
        iomodel->FetchDataToInput(elements,VzEnum);
        InputUpdateFromConstantx(elements,0.,VxMeshEnum);
        InputUpdateFromConstantx(elements,0.,VyMeshEnum);
        InputUpdateFromConstantx(elements,0.,VzMeshEnum);
        if(dakota_analysis){
                elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
                elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
                elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
                elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
                elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
        }

        /*Free data: */
        iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
}/*}}}*/
void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/

        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
        ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
}/*}}}*/
void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/

        /*Intermediary*/
        int        count;
        int        M,N;
        bool       spcpresent = false;
        IssmDouble heatcapacity;
        IssmDouble referencetemperature;

        /*Output*/
        IssmDouble *spcvector  = NULL;
        IssmDouble* times=NULL;
        IssmDouble* values=NULL;

        /*Fetch parameters: */
        iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
        iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);

        /*return if 2d mesh*/
        if(iomodel->meshtype==Mesh2DhorizontalEnum) return;

        /*Fetch data: */
        iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);

        //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS 
        /*Transient or static?:*/
        if(M==iomodel->numberofvertices){
                /*static: just create Constraints objects*/
                count=0;

                for(int i=0;i<iomodel->numberofvertices;i++){
                        /*keep only this partition's nodes:*/
                        if((iomodel->my_vertices[i])){

                                if (!xIsNan<IssmDouble>(spcvector[i])){

                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
                                        count++;

                                }
                        }
                }
        }
        else if (M==(iomodel->numberofvertices+1)){
                /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve 
                 * various times and values to initialize an SpcTransient object: */
                count=0;

                /*figure out times: */
                times=xNew<IssmDouble>(N);
                for(int j=0;j<N;j++){
                        times[j]=spcvector[(M-1)*N+j];
                }

                /*Create constraints from x,y,z: */
                for(int i=0;i<iomodel->numberofvertices;i++){

                        /*keep only this partition's nodes:*/
                        if((iomodel->my_vertices[i])){

                                /*figure out times and values: */
                                values=xNew<IssmDouble>(N);
                                spcpresent=false;
                                for(int j=0;j<N;j++){
                                        values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
                                        if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
                                }

                                if(spcpresent){
                                        constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum));
                                        count++;
                                }
                                xDelete<IssmDouble>(values);
                        }
                }
        }
        else{
                _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
        }

        /*Free ressources:*/
        iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum);
        xDelete<IssmDouble>(times);
        xDelete<IssmDouble>(values);
}/*}}}*/
void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/

        /*No loads */
}/*}}}*/