1 | /*!\file ThicknessAbsMisfitx
|
---|
2 | * \brief: compute misfit between observations and model
|
---|
3 | */
|
---|
4 |
|
---|
5 | #include "./ThicknessAbsMisfitx.h"
|
---|
6 |
|
---|
7 | #include "../../shared/shared.h"
|
---|
8 | #include "../../toolkits/toolkits.h"
|
---|
9 |
|
---|
10 | void ThicknessAbsMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
|
---|
11 |
|
---|
12 | /*output: */
|
---|
13 | IssmDouble J=0.;
|
---|
14 | IssmDouble J_sum;
|
---|
15 |
|
---|
16 | /*Compute Misfit: */
|
---|
17 | for(int i=0;i<elements->Size();i++){
|
---|
18 | Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
|
---|
19 | J+=ThicknessAbsMisfit(element);
|
---|
20 | }
|
---|
21 |
|
---|
22 | /*Sum all J from all cpus of the cluster:*/
|
---|
23 | ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
|
---|
24 | ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
|
---|
25 | J=J_sum;
|
---|
26 |
|
---|
27 | /*Assign output pointers: */
|
---|
28 | *pJ=J;
|
---|
29 | }
|
---|
30 |
|
---|
31 | IssmDouble ThicknessAbsMisfit(Element* element){
|
---|
32 |
|
---|
33 | IssmDouble thickness,thicknessobs,weight;
|
---|
34 | IssmDouble Jelem=0.;
|
---|
35 | IssmDouble misfit,Jdet;
|
---|
36 | IssmDouble* xyz_list = NULL;
|
---|
37 |
|
---|
38 | /*If on water, return 0: */
|
---|
39 | if(!element->IsIceInElement()) return 0.;
|
---|
40 |
|
---|
41 | /* Get node coordinates*/
|
---|
42 | element->GetVerticesCoordinates(&xyz_list);
|
---|
43 |
|
---|
44 | /*Retrieve all inputs we will be needing: */
|
---|
45 | Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
|
---|
46 | Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input);
|
---|
47 | Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input);
|
---|
48 |
|
---|
49 | /* Start looping on the number of gaussian points: */
|
---|
50 | Gauss* gauss=element->NewGauss(2);
|
---|
51 | for(int ig=gauss->begin();ig<gauss->end();ig++){
|
---|
52 |
|
---|
53 | gauss->GaussPoint(ig);
|
---|
54 |
|
---|
55 | /* Get Jacobian determinant: */
|
---|
56 | element->JacobianDeterminant(&Jdet,xyz_list,gauss);
|
---|
57 |
|
---|
58 | /*Get all parameters at gaussian point*/
|
---|
59 | weights_input->GetInputValue(&weight,gauss,ThicknessAbsMisfitEnum);
|
---|
60 | thickness_input->GetInputValue(&thickness,gauss);
|
---|
61 | thicknessobs_input->GetInputValue(&thicknessobs,gauss);
|
---|
62 |
|
---|
63 | /*Compute ThicknessAbsMisfitEnum*/
|
---|
64 | Jelem+=0.5*(thickness-thicknessobs)*(thickness-thicknessobs)*weight*Jdet*gauss->weight;
|
---|
65 | }
|
---|
66 |
|
---|
67 | /*clean up and Return: */
|
---|
68 | xDelete<IssmDouble>(xyz_list);
|
---|
69 | delete gauss;
|
---|
70 | return Jelem;
|
---|
71 | }
|
---|