1 | #include "./OceantransportAnalysis.h"
|
---|
2 | #include <math.h>
|
---|
3 | #include "../toolkits/toolkits.h"
|
---|
4 | #include "../classes/classes.h"
|
---|
5 | #include "../classes/Inputs/TransientInput.h"
|
---|
6 | #include "../classes/Inputs/TriaInput.h"
|
---|
7 | #include "../classes/gauss/Gauss.h"
|
---|
8 | #include "../shared/shared.h"
|
---|
9 | #include "../modules/modules.h"
|
---|
10 |
|
---|
11 | /*Model processing*/
|
---|
12 | void OceantransportAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
|
---|
13 | /*No constraints*/
|
---|
14 | }/*}}}*/
|
---|
15 | void OceantransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
|
---|
16 | /*No loads*/
|
---|
17 | }/*}}}*/
|
---|
18 | void OceantransportAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
|
---|
19 | ::CreateNodes(nodes,iomodel,OceantransportAnalysisEnum,P1Enum);
|
---|
20 | }/*}}}*/
|
---|
21 | int OceantransportAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
|
---|
22 | return 3;
|
---|
23 | }/*}}}*/
|
---|
24 | void OceantransportAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
|
---|
25 |
|
---|
26 | int nature=0;
|
---|
27 | bool isdakota=0;
|
---|
28 |
|
---|
29 | /*Update elements: */
|
---|
30 | int counter=0;
|
---|
31 | for(int i=0;i<iomodel->numberofelements;i++){
|
---|
32 | if(iomodel->my_elements[i]){
|
---|
33 | Element* element=(Element*)elements->GetObjectByOffset(counter);
|
---|
34 | element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
|
---|
35 | counter++;
|
---|
36 | }
|
---|
37 | }
|
---|
38 |
|
---|
39 | /*Plug inputs into element:*/
|
---|
40 | iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_at_sea_floor", OceantransportSpcbottompressureEnum);
|
---|
41 | iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_above_geoid", OceantransportSpcdslEnum);
|
---|
42 | iomodel->FetchDataToInput(inputs,elements,"md.dsl.global_average_thermosteric_sea_level",OceantransportSpcstrEnum);
|
---|
43 |
|
---|
44 | /*Initialize sea level cumulated sea level loads :*/
|
---|
45 | iomodel->ConstantToInput(inputs,elements,0.,AccumulatedDeltaBottomPressureEnum,P1Enum);
|
---|
46 | iomodel->ConstantToInput(inputs,elements,0.,OldAccumulatedDeltaBottomPressureEnum,P1Enum);
|
---|
47 | iomodel->FetchDataToInput(inputs,elements,"md.initialization.bottompressure",BottomPressureEnum);
|
---|
48 | iomodel->FetchDataToInput(inputs,elements,"md.initialization.dsl",DslEnum);
|
---|
49 | iomodel->FetchDataToInput(inputs,elements,"md.initialization.str",StrEnum);
|
---|
50 |
|
---|
51 | }/*}}}*/
|
---|
52 | void OceantransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
|
---|
53 |
|
---|
54 | int dslmodel=0;
|
---|
55 | int numoutputs;
|
---|
56 | char** requestedoutputs = NULL;
|
---|
57 |
|
---|
58 | /*Deal with dsl multi-model ensembles: {{{*/
|
---|
59 | iomodel->FetchData(&dslmodel,"md.dsl.model");
|
---|
60 | if(dslmodel==2){
|
---|
61 | IssmDouble modelid;
|
---|
62 | int nummodels;
|
---|
63 |
|
---|
64 | /*create double param, not int param, because Dakota will be updating it as
|
---|
65 | * a double potentially: */
|
---|
66 | iomodel->FetchData(&modelid,"md.dsl.modelid");
|
---|
67 | parameters->AddObject(new DoubleParam(DslModelidEnum,modelid));
|
---|
68 | parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
|
---|
69 | iomodel->FetchData(&nummodels,"md.dsl.nummodels");
|
---|
70 |
|
---|
71 | /*quick checks: */
|
---|
72 | if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!");
|
---|
73 | if(modelid<=0 || modelid>nummodels)_error_("modelid field in dslmme object of md.dsl field should be between 1 and the number of ensemble runs!");
|
---|
74 | } /*}}}*/
|
---|
75 | /*Requested outputs {{{*/
|
---|
76 | iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.solidearth.requested_outputs");
|
---|
77 | if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelchangeRequestedOutputsEnum,requestedoutputs,numoutputs));
|
---|
78 | iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs");
|
---|
79 | /*}}}*/
|
---|
80 |
|
---|
81 | }/*}}}*/
|
---|
82 |
|
---|
83 | /*Finite Element Analysis*/
|
---|
84 | void OceantransportAnalysis::Core(FemModel* femmodel){/*{{{*/
|
---|
85 | _error_("not implemented");
|
---|
86 | }/*}}}*/
|
---|
87 | void OceantransportAnalysis::PreCore(FemModel* femmodel){/*{{{*/
|
---|
88 | _error_("not implemented");
|
---|
89 | }/*}}}*/
|
---|
90 | ElementVector* OceantransportAnalysis::CreateDVector(Element* element){/*{{{*/
|
---|
91 | /*Default, return NULL*/
|
---|
92 | return NULL;
|
---|
93 | }/*}}}*/
|
---|
94 | ElementMatrix* OceantransportAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
|
---|
95 | _error_("Not implemented");
|
---|
96 | }/*}}}*/
|
---|
97 | ElementMatrix* OceantransportAnalysis::CreateKMatrix(Element* element){/*{{{*/
|
---|
98 | _error_("not implemented yet");
|
---|
99 | }/*}}}*/
|
---|
100 | ElementVector* OceantransportAnalysis::CreatePVector(Element* element){/*{{{*/
|
---|
101 | _error_("not implemented yet");
|
---|
102 | }/*}}}*/
|
---|
103 | void OceantransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
|
---|
104 |
|
---|
105 | /*do nothing:*/
|
---|
106 | return;
|
---|
107 |
|
---|
108 | }/*}}}*/
|
---|
109 | void OceantransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
|
---|
110 | _error_("Not implemented yet");
|
---|
111 | }/*}}}*/
|
---|
112 | void OceantransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
|
---|
113 |
|
---|
114 | /*Fetch number of nodes and dof for this finite element*/
|
---|
115 | IssmDouble time;
|
---|
116 | IssmDouble bp[3];
|
---|
117 | IssmDouble dsl[3];
|
---|
118 | IssmDouble str;
|
---|
119 | int numnodes = element->GetNumberOfNodes();
|
---|
120 |
|
---|
121 | element->parameters->FindParam(&time,TimeEnum);
|
---|
122 |
|
---|
123 | TriaInput* str_input=xDynamicCast<TriaInput*>(element->GetInput(OceantransportSpcstrEnum,time)); _assert_(str_input);
|
---|
124 | TriaInput* bp_input=xDynamicCast<TriaInput*>(element->GetInput(OceantransportSpcbottompressureEnum,time)); _assert_(bp_input);
|
---|
125 | TriaInput* dsl_input=xDynamicCast<TriaInput*>(element->GetInput(OceantransportSpcdslEnum,time)); _assert_(dsl_input);
|
---|
126 |
|
---|
127 | Gauss* gauss=element->NewGauss();
|
---|
128 | for(int iv=0;iv<3;iv++){
|
---|
129 | gauss->GaussVertex(iv);
|
---|
130 | bp_input->GetInputValue(&bp[iv],gauss);
|
---|
131 | dsl_input->GetInputValue(&dsl[iv],gauss);
|
---|
132 | }
|
---|
133 | str_input->GetInputAverage(&str);
|
---|
134 | delete gauss;
|
---|
135 |
|
---|
136 | /*Add str bp and dsl levelsets as inputs to the tria element: */
|
---|
137 | element->AddInput(BottomPressureEnum,bp,P1Enum);
|
---|
138 | element->AddInput(DslEnum,dsl,P1Enum);
|
---|
139 | element->AddInput(StrEnum,&str,P0Enum);
|
---|
140 |
|
---|
141 | }/*}}}*/
|
---|
142 | void OceantransportAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
|
---|
143 | /*Default, do nothing*/
|
---|
144 | return;
|
---|
145 | }/*}}}*/
|
---|