![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 13333
- Timestamp:
- 09/11/12 16:02:12 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp
r13328 r13333 17 17 parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); 18 18 if(isautodiff){ 19 20 21 22 23 24 25 26 27 /*retrieve state variable: */28 29 parameters->FindParam(&axp,&num_independents,AutodiffXpEnum);30 /* driver argument */31 double *xp=xNew<double>(num_independents);32 for(int i=0;i<num_independents;i++){33 xp[i]=reCast<double,IssmDouble>(axp[i]);34 }35 /* get the dimension for the solverx arguments*/36 37 parameters->FindParam(&configuration_type,ConfigurationTypeEnum);38 int solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum);39 int edf_n=solveSize*(solveSize+1), edf_m=solveSize;19 #ifdef _HAVE_ADOLC_ 20 int num_dependents; 21 parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); 22 int num_independents; 23 parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); 24 if(!(num_dependents*num_independents)) return; 25 int numberofvertices; 26 parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum); 27 /*retrieve state variable: */ 28 IssmDouble *axp = NULL; 29 parameters->FindParam(&axp,&num_independents,AutodiffXpEnum); 30 /* driver argument */ 31 double *xp=xNew<double>(num_independents); 32 for(int i=0;i<num_independents;i++){ 33 xp[i]=reCast<double,IssmDouble>(axp[i]); 34 } 35 /* get the dimension for the solverx arguments*/ 36 int configuration_type; 37 parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 38 int solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum); 39 int edf_n=solveSize*(solveSize+1), edf_m=solveSize; 40 40 41 /*get the EDF pointer:*/42 ext_diff_fct *anEDF_for_solverx_p=dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;43 44 /*set the forward method function pointers: */45 anEDF_for_solverx_p->zos_forward=EDF_for_solverx;46 anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;47 anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx;48 49 50 /*allocate the space for the parameters to invoke the forward methods:*/51 anEDF_for_solverx_p->dp_x=xNew<double>(edf_n);52 anEDF_for_solverx_p->dp_X=xNew<double>(edf_n);53 anEDF_for_solverx_p->dpp_X=xNew<double>(edf_n, num_independents);54 anEDF_for_solverx_p->dp_y=xNew<double>(edf_m);55 anEDF_for_solverx_p->dp_Y=xNew<double>(edf_m);56 anEDF_for_solverx_p->dpp_Y=xNew<double>(edf_m, num_independents);57 anEDF_for_solverx_p->dp_U=xNew<double>(edf_m);58 anEDF_for_solverx_p->dpp_U=xNew<double>(num_dependents,edf_m);59 60 61 /* Call AD driver:*/62 63 64 65 66 67 68 69 70 71 72 73 41 /*get the EDF pointer:*/ 42 ext_diff_fct *anEDF_for_solverx_p=dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; 43 44 /*set the forward method function pointers: */ 45 anEDF_for_solverx_p->zos_forward=EDF_for_solverx; 46 anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx; 47 anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx; 48 anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; 49 // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; 50 /*allocate the space for the parameters to invoke the forward methods:*/ 51 anEDF_for_solverx_p->dp_x=xNew<double>(edf_n); 52 anEDF_for_solverx_p->dp_X=xNew<double>(edf_n); 53 anEDF_for_solverx_p->dpp_X=xNew<double>(edf_n, num_independents); 54 anEDF_for_solverx_p->dp_y=xNew<double>(edf_m); 55 anEDF_for_solverx_p->dp_Y=xNew<double>(edf_m); 56 anEDF_for_solverx_p->dpp_Y=xNew<double>(edf_m, num_independents); 57 anEDF_for_solverx_p->dp_U=xNew<double>(edf_m); 58 anEDF_for_solverx_p->dpp_U=xNew<double>(num_dependents,edf_m); 59 anEDF_for_solverx_p->dp_Z=xNew<double>(edf_n); 60 anEDF_for_solverx_p->dpp_Z=xNew<double>(num_dependents,edf_n); 61 /* Call AD driver:*/ 62 // single direction: 63 double *tangentDir=xNewZeroInit<double>(num_independents); 64 unsigned int anIndepNum=55; // <---------- make this selectable via config 65 tangentDir[anIndepNum]=1.0; 66 double *theJacVecProduct=xNew<double>(num_dependents); 67 double *theOutput=xNew<double>(num_dependents); 68 for (int i=0; i<num_independents; ++i) { 69 std::cout << "fos_forward xp[" << i << "]=" << xp[i] << std::endl; 70 } 71 if (fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, theJacVecProduct )) 72 _error_("fos_forward returned non-zero error code"); 73 results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,theJacVecProduct,num_dependents,1,1,0.0)); 74 74 75 75 #if 0 76 77 jacobian(1,num_dependents,num_independents,xp,J);78 79 76 double **J=xNew<double>(num_dependents,num_independents); 77 jacobian(1,num_dependents,num_independents,xp,J); 78 results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,*J,num_independents,num_dependents,1,1,0.0)); 79 xDelete(J); 80 80 #endif 81 81 82 83 84 85 86 87 88 89 90 91 92 93 /*Free ressources: */94 xDelete(xp);95 xDelete(axp); // did we allocate this?82 /* delete the allocated space for the parameters:*/ 83 xDelete(anEDF_for_solverx_p->dp_x); 84 xDelete(anEDF_for_solverx_p->dp_X); 85 xDelete(anEDF_for_solverx_p->dpp_X); 86 xDelete(anEDF_for_solverx_p->dp_y); 87 xDelete(anEDF_for_solverx_p->dp_Y); 88 xDelete(anEDF_for_solverx_p->dpp_Y); 89 xDelete(anEDF_for_solverx_p->dp_U); 90 xDelete(anEDF_for_solverx_p->dpp_U); 91 xDelete(anEDF_for_solverx_p->dp_Z); 92 xDelete(anEDF_for_solverx_p->dpp_Z); 93 /*Free ressources: */ 94 xDelete(xp); 95 xDelete(axp); // did we allocate this? 96 96 #else 97 _error_("Should not be requesting AD drivers when an AD library is not available!");97 _error_("Should not be requesting AD drivers when an AD library is not available!"); 98 98 #endif 99 99 }
Note:
See TracChangeset
for help on using the changeset viewer.
