source: issm/oecreview/Archive/18296-19100/ISSM-18895-18896.diff

Last change on this file was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 6.5 KB
RevLine 
[19102]1Index: ../trunk-jpl/src/c/Makefile.am
2===================================================================
3--- ../trunk-jpl/src/c/Makefile.am (revision 18895)
4+++ ../trunk-jpl/src/c/Makefile.am (revision 18896)
5@@ -345,6 +345,7 @@
6 ./cores/WrapperCorePointerFromSolutionEnum.cpp\
7 ./cores/CorePointerFromSolutionEnum.cpp\
8 ./cores/ad_core.cpp\
9+ ./cores/adgradient_core.cpp\
10 ./main/EnvironmentInit.cpp\
11 ./main/EnvironmentFinalize.cpp\
12 ./analyses/EnumToAnalysis.h\
13Index: ../trunk-jpl/src/c/cores/adgradient_core.cpp
14===================================================================
15--- ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 0)
16+++ ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 18896)
17@@ -0,0 +1,145 @@
18+/*!\file adgradient_core
19+ * \brief: compute gradient for all scalar depenendents, then sum them up as output. This relies mainly on the fos_reverse
20+ * driver, hence the incapacity to merge this with ad_core.cpp.
21+ */
22+
23+/*Includes: {{{*/
24+#ifdef HAVE_CONFIG_H
25+ #include <config.h>
26+#else
27+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
28+#endif
29+
30+#include <set>
31+#include "./cores.h"
32+#include "../toolkits/toolkits.h"
33+#include "../classes/classes.h"
34+#include "../shared/shared.h"
35+#include "../modules/modules.h"
36+#include "../solutionsequences/solutionsequences.h"
37+/*}}}*/
38+
39+void adgradient_core(FemModel* femmodel){
40+
41+ /*diverse: */
42+ int i;
43+ int dummy;
44+ int num_dependents=0;
45+ int num_dependents_old=0;
46+ int num_independents=0;
47+ bool isautodiff = false;
48+ int aDepIndex=0;
49+ int my_rank=IssmComm::GetRank();
50+
51+ /*state variables: */
52+ IssmDouble *axp = NULL;
53+ IssmPDouble *xp = NULL;
54+
55+ /*intermediary: */
56+ IssmPDouble *aWeightVector=NULL;
57+ IssmPDouble *weightVectorTimesJac=NULL;
58+
59+ /*output: */
60+ IssmPDouble *totalgradient=NULL;
61+
62+ /*AD mode on?: */
63+ femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
64+
65+ if(isautodiff){
66+
67+ #ifdef _HAVE_ADOLC_
68+ if(VerboseAutodiff())_printf0_(" start ad core\n");
69+
70+ /*First, stop tracing: */
71+ trace_off();
72+
73+ /*retrieve parameters: */
74+ femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
75+ femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
76+
77+ /*if no dependents, no point in running a driver: */
78+ if(!(num_dependents*num_independents)) return;
79+
80+ /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
81+ num_dependents_old=num_dependents;
82+ if (my_rank!=0){
83+ num_dependents=0; num_independents=0;
84+ }
85+
86+ /*retrieve state variable: */
87+ femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
88+
89+ /* driver argument */
90+ xp=xNew<double>(num_independents);
91+ for(i=0;i<num_independents;i++){
92+ xp[i]=reCast<double,IssmDouble>(axp[i]);
93+ }
94+
95+ /*get the EDF pointer:*/
96+ ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
97+
98+ /* these are always needed regardless of the interpreter */
99+ anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
100+ anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
101+
102+ /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
103+ * as to generate num_dependents gradients: */
104+
105+ /*Initialize outputs: */
106+ totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
107+
108+ for(aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
109+
110+ /*initialize direction index in the weights vector: */
111+ aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
112+ if (my_rank==0) aWeightVector[aDepIndex]=1.0;
113+
114+ /*initialize output gradient: */
115+ weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
116+
117+ /*set the forward method function pointer: */
118+ #ifdef _HAVE_GSL_
119+ anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
120+ #endif
121+ #ifdef _HAVE_MUMPS_
122+ anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
123+ #endif
124+
125+ anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
126+ anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
127+
128+ /*call driver: */
129+ fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
130+
131+ /*Add to totalgradient: */
132+ if(my_rank==0)for(i=0;i<num_independents;i++)totalgradient[i]+=weightVectorTimesJac[i];
133+
134+ /*free resources :*/
135+ xDelete(weightVectorTimesJac);
136+ xDelete(aWeightVector);
137+ }
138+
139+ /*add totalgradient to results*/
140+ femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0));
141+
142+ if(VerboseAutodiff())_printf0_(" end ad core\n");
143+
144+ /* delete the allocated space for the parameters and free ressources:{{{*/
145+ xDelete(anEDF_for_solverx_p->dp_x);
146+ xDelete(anEDF_for_solverx_p->dp_X);
147+ xDelete(anEDF_for_solverx_p->dpp_X);
148+ xDelete(anEDF_for_solverx_p->dp_y);
149+ xDelete(anEDF_for_solverx_p->dp_Y);
150+ xDelete(anEDF_for_solverx_p->dpp_Y);
151+ xDelete(anEDF_for_solverx_p->dp_U);
152+ xDelete(anEDF_for_solverx_p->dpp_U);
153+ xDelete(anEDF_for_solverx_p->dp_Z);
154+ xDelete(anEDF_for_solverx_p->dpp_Z);
155+ xDelete(xp);
156+ xDelete(totalgradient);
157+ xDelete(axp); /*}}}*/
158+ #else
159+ _error_("Should not be requesting AD drivers when an AD library is not available!");
160+ #endif
161+ }
162+}
163Index: ../trunk-jpl/src/c/cores/cores.h
164===================================================================
165--- ../trunk-jpl/src/c/cores/cores.h (revision 18895)
166+++ ../trunk-jpl/src/c/cores/cores.h (revision 18896)
167@@ -42,6 +42,7 @@
168 void transient_core(FemModel* femmodel);
169 void dakota_core(FemModel* femmodel);
170 void ad_core(FemModel* femmodel);
171+void adgradient_core(FemModel* femmodel);
172 void dummy_core(FemModel* femmodel);
173 void gia_core(FemModel* femmodel);
174 void damage_core(FemModel* femmodel);
175Index: ../trunk-jpl/src/c/cores/controlad_core.cpp
176===================================================================
177--- ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 18895)
178+++ ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 18896)
179@@ -215,7 +215,7 @@
180 _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
181
182 /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
183- ad_core(femmodel);
184+ adgradient_core(femmodel);
185
186 if(IssmComm::GetRank()==0){
187 IssmPDouble* G_temp=NULL;
Note: See TracBrowser for help on using the repository browser.