Changeset 13314

Timestamp:

09/10/12 23:06:01 (13 years ago)

Author:

utke

Message:

setup for fos_forward

File:

• issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp

@@ -12 +12 @@
 
 void AutodiffDriversx(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,Results* results){
-
-        int  i                  ,j;
         bool isautodiff       = false;
-        int  num_dependents;
-        int  num_independents;
-        int  numberofvertices;
-
-        /*outputs: */
-        IssmDouble  *matJ = NULL;
-        IssmDouble  *axp  = NULL;
-        double     **J    = NULL;
-        double      *xp   = NULL;
-                
-        /*diverse: */
-        int configuration_type;
-        int solveSize;
-        int edf_n,edf_m;
-
         /*AD mode on?: */
         parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
-
         if(isautodiff){
-                #ifdef _HAVE_ADOLC_
-                parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
-                parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
-                parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
-
-                if(!(num_dependents*num_independents)) return;
-                
+                #ifdef _HAVE_ADOLC_
+                int  num_dependents;
+                parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
+                int  num_independents;
+                parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
+                if(!(num_dependents*num_independents)) return;
+                int  numberofvertices;
+                parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
                 /*retrieve state variable: */
+                IssmDouble  *axp  = NULL;
                 parameters->FindParam(&axp,&num_independents,AutodiffXpEnum);
-
-                /*allocate driver results: */
-                J=xNew<double*>(num_independents);
-                xp=xNew<double>(num_independents); 
-
-                for(i=0;i<num_independents;i++){
-                        J[i]=xNew<double>(num_dependents);
+                /* driver argument */
+                double *xp=xNew<double>(num_independents);
+                for(int i=0;i<num_independents;i++){
                         xp[i]=reCast<double,IssmDouble>(axp[i]);
                 }
-
                 /* get the dimension for the solverx arguments*/
+                int configuration_type;
                 parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-                solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum);
-                edf_n=solveSize*(solveSize+1), edf_m=solveSize;
+                int solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum);
+                int edf_n=solveSize*(solveSize+1), edf_m=solveSize;
 
                 /*get the EDF pointer:*/
@@ -62 +42 @@
                 
                 /*set the forward method function pointers: */
+                anEDF_for_solverx_p->zos_forward=EDF_for_solverx;
                 anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
                 anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx;
-
+                anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
+                // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
                 /*allocate the space for the parameters to invoke the forward methods:*/
                 anEDF_for_solverx_p->dp_x=xNew<double>(edf_n);
@@ -72 +54 @@
                 anEDF_for_solverx_p->dp_Y=xNew<double>(edf_m);
                 anEDF_for_solverx_p->dpp_Y=xNew<double>(edf_m, num_independents);
+                anEDF_for_solverx_p->dp_U=xNew<double>(edf_m);
+                anEDF_for_solverx_p->dpp_U=xNew<double>(num_dependents,edf_m);
+                anEDF_for_solverx_p->dp_Z=xNew<double>(edf_n);
+                anEDF_for_solverx_p->dpp_Z=xNew<double>(num_dependents,edf_n);
+                /* Call AD driver:*/
+                // single direction:
+                double *tangentDir=xNewZeroInit<double>(num_independents);
+                unsigned int anIndepNum=55; // <---------- make this selectable via config
+                tangentDir[anIndepNum]=1.0;
+                double *theJacVecProduct=xNew<double>(num_dependents);
+                double *theOutput=xNew<double>(num_dependents);
+                if (fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, theJacVecProduct ))
+                  _error_("fos_forward returned non-zero error code");
+                IssmDouble *theJacVecProductA=xNew<IssmDouble>(num_dependents);
+                for (int i=0;i<num_dependents;++i)
+                  theJacVecProductA[i]=theJacVecProduct[i];
+                results->AddObject(new DoubleVecExternalResult(results->Size()+1,AutodiffJacobianEnum,theJacVecProductA,num_dependents,1,0.0));
 
-                /* Call AD driver:*/
+#if 0
+                double **J=xNew<double>(num_dependents,num_independents);
                 jacobian(1,num_dependents,num_independents,xp,J);
+                IssmDouble *matJ=xNew<IssmDouble>(num_dependents*num_independents);
+                for (int i=0;i<num_dependents*num_independents;++i)
+                  matJ[i]=J[i];
+                results->AddObject(new DoubleMatExternalResult(results->Size()+1,AutodiffJacobianEnum,matJ,num_independents,num_dependents,1,0.0));
+                xDelete(matJ);
+                xDelete(J);
+#endif
 
-                /* delete the allocated space for the parameters:*/
-                xDelete(anEDF_for_solverx_p->dp_x);
-                xDelete(anEDF_for_solverx_p->dp_X);
-                xDelete(anEDF_for_solverx_p->dpp_X,edf_n);
-                xDelete(anEDF_for_solverx_p->dp_y);
-                xDelete(anEDF_for_solverx_p->dp_Y);
-                xDelete(anEDF_for_solverx_p->dpp_Y,edf_m);
-                
-                /*Add jacobian matrix to results: */
-                matJ=xNew<IssmDouble>(num_dependents*num_independents);
-                for (i=0;i<num_independents;++i){
-                        double* rowi=J[i];
-                        for (j=0;j<num_dependents;++j){
-                                *(matJ+num_dependents*i+j)=rowi[j];
-                        }
-                }
-                results->AddObject(new DoubleMatExternalResult(results->Size()+1,AutodiffJacobianEnum,matJ,num_independents,num_dependents,1,0.0));
-
+                /* delete the allocated space for the parameters:*/
+                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);
                 /*Free ressources: */
-                xDelete<double*>(J);
-                xDelete<double>(xp);
-                xDelete<IssmDouble>(matJ);
-                xDelete<IssmDouble>(axp);
+                xDelete(xp);
+                xDelete(axp); // did we allocate this?
                 #else
                 _error_("Should not be requesting AD drivers when an AD library is not available!");