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

../trunk-jpl/src/c/Makefile.am

@@ -345 +345 @@
                                         ./cores/WrapperCorePointerFromSolutionEnum.cpp\
                                         ./cores/CorePointerFromSolutionEnum.cpp\
                                         ./cores/ad_core.cpp\
+                                        ./cores/adgradient_core.cpp\
                                         ./main/EnvironmentInit.cpp\
                                         ./main/EnvironmentFinalize.cpp\
                                         ./analyses/EnumToAnalysis.h\

../trunk-jpl/src/c/cores/adgradient_core.cpp

@@ +1 @@
+/*!\file adgradient_core
+ * \brief: compute gradient for all scalar depenendents, then sum them up as output. This relies mainly on the fos_reverse 
+ * driver, hence the incapacity to merge this with ad_core.cpp.
+ */
+
+/*Includes: {{{*/
+#ifdef HAVE_CONFIG_H
+        #include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <set>
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+/*}}}*/
+
+void adgradient_core(FemModel* femmodel){
+
+        /*diverse: */
+        int     i;
+        int     dummy;
+        int     num_dependents=0;
+        int     num_dependents_old=0;
+        int     num_independents=0;
+        bool    isautodiff       = false;
+        int     aDepIndex=0;
+        int     my_rank=IssmComm::GetRank();
+
+        /*state variables: */
+        IssmDouble *axp = NULL;
+        IssmPDouble     *xp  = NULL;
+
+        /*intermediary: */
+        IssmPDouble *aWeightVector=NULL;
+        IssmPDouble *weightVectorTimesJac=NULL;
+        
+        /*output: */
+        IssmPDouble *totalgradient=NULL;
+
+        /*AD mode on?: */
+        femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
+
+        if(isautodiff){
+
+                #ifdef _HAVE_ADOLC_
+                        if(VerboseAutodiff())_printf0_("   start ad core\n"); 
+
+                        /*First, stop tracing: */
+                        trace_off();
+                
+                        /*retrieve parameters: */
+                        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
+                        femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
+                        
+                        /*if no dependents, no point in running a driver: */
+                        if(!(num_dependents*num_independents)) return;
+
+                        /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
+                        num_dependents_old=num_dependents;
+                        if (my_rank!=0){
+                                num_dependents=0; num_independents=0;
+                        }
+                        
+                        /*retrieve state variable: */
+                        femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
+                        
+                        /* driver argument */
+                        xp=xNew<double>(num_independents);
+                        for(i=0;i<num_independents;i++){
+                                xp[i]=reCast<double,IssmDouble>(axp[i]);
+                        }
+
+                        /*get the EDF pointer:*/
+                        ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
+
+                        /* these are always needed regardless of the interpreter */
+                        anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
+                        anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
+
+                        /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so 
+                         * as to generate num_dependents gradients: */
+
+                        /*Initialize outputs: */
+                        totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
+                        
+                        for(aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
+
+                                /*initialize direction index in the weights vector: */
+                                aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
+                                if (my_rank==0) aWeightVector[aDepIndex]=1.0;
+                                
+                                /*initialize output gradient: */
+                                weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
+
+                                /*set the forward method function pointer: */
+                                #ifdef _HAVE_GSL_
+                                anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
+                                #endif
+                                #ifdef _HAVE_MUMPS_
+                                anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
+                                #endif
+
+                                anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
+                                anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
+
+                                /*call driver: */
+                                fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
+
+                                /*Add to totalgradient: */
+                                if(my_rank==0)for(i=0;i<num_independents;i++)totalgradient[i]+=weightVectorTimesJac[i];
+
+                                /*free resources :*/
+                                xDelete(weightVectorTimesJac);
+                                xDelete(aWeightVector);
+                        }
+                
+                        /*add totalgradient to results*/
+                        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0));
+
+                        if(VerboseAutodiff())_printf0_("   end ad core\n");
+                        
+                        /* delete the allocated space for the parameters and free ressources:{{{*/
+                        xDelete(anEDF_for_solverx_p->dp_x);
+                        xDelete(anEDF_for_solverx_p->dp_X);
+                        xDelete(anEDF_for_solverx_p->dpp_X);
+                        xDelete(anEDF_for_solverx_p->dp_y);
+                        xDelete(anEDF_for_solverx_p->dp_Y);
+                        xDelete(anEDF_for_solverx_p->dpp_Y);
+                        xDelete(anEDF_for_solverx_p->dp_U);
+                        xDelete(anEDF_for_solverx_p->dpp_U);
+                        xDelete(anEDF_for_solverx_p->dp_Z);
+                        xDelete(anEDF_for_solverx_p->dpp_Z);
+                        xDelete(xp);
+                        xDelete(totalgradient);
+                        xDelete(axp); /*}}}*/
+                #else
+                        _error_("Should not be requesting AD drivers when an AD library is not available!");
+                #endif
+        }
+}

../trunk-jpl/src/c/cores/cores.h

@@ -42 +42 @@
 void transient_core(FemModel* femmodel);
 void dakota_core(FemModel* femmodel);
 void ad_core(FemModel* femmodel);
+void adgradient_core(FemModel* femmodel);
 void dummy_core(FemModel* femmodel);
 void gia_core(FemModel* femmodel);
 void damage_core(FemModel* femmodel);

../trunk-jpl/src/c/cores/controlad_core.cpp

@@ -215 +215 @@
         _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
         
         /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
-        ad_core(femmodel); 
+        adgradient_core(femmodel); 
 
         if(IssmComm::GetRank()==0){
                 IssmPDouble* G_temp=NULL;