Index: /issm/trunk-jpl/src/c/cores/controlad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlad_core.cpp (revision 18941) +++ /issm/trunk-jpl/src/c/cores/controlad_core.cpp (revision 18942) @@ -25,8 +25,11 @@ /*Cost function prototype*/ void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs); +FemModel* presimulad(int* pintn, double** pX, FemModel* femmodel); +void postsimulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs); /*}}}*/ void controlad_core(FemModel* femmodel){ /*{{{*/ /*Intermediaries*/ + FemModel* femmodelad=NULL; int i; long omode; @@ -34,22 +37,16 @@ IssmDouble dxmind,gttold; int maxsteps,maxiter; - int intn,numberofvertices,num_controls,solution_type; - IssmDouble *scaling_factors = NULL; + int intn,solution_type; IssmPDouble *X = NULL; IssmDouble *Xd = NULL; - IssmDouble *Gd = NULL; IssmPDouble *G = NULL; - bool onsid=true; /*Recover some parameters*/ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum); femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum); femmodel->parameters->FindParam(&dxmind,InversionDxminEnum); dxmin=reCast(dxmind); femmodel->parameters->FindParam(>told,InversionGttolEnum); gttol=reCast(gttold); - femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); femmodel->parameters->SetParam(false,SaveResultsEnum); - numberofvertices=femmodel->vertices->NumberOfVertices(); /*Initialize M1QN3 parameters*/ @@ -71,26 +68,10 @@ long nsim = long(maxiter);/*Maximum number of function calls*/ - /*Get initial guess*/ - Vector *Xad = NULL; - GetVectorFromControlInputsx(&Xad,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value",onsid); - Xd = Xad->ToMPISerial(); - Xad->GetSize(&intn); - X=xNew(intn); - for(i=0;i(Xd[i]); - delete Xad; - _assert_(intn==numberofvertices*num_controls); - - /*Get problem dimension and initialize gradient and initial guess*/ + /*Run the first part of simulad, in order to get things started!:*/ + femmodelad=presimulad(&intn,&X,femmodel); + + /*Get problem dimension and initialize gradient: */ long n = long(intn); G = xNew(n); - Gd = xNew(n); - - /*Scale control for M1QN3*/ - for(int i=0;i(scaling_factors[c]); - } - } /*Allocate m1qn3 working arrays (see doc)*/ @@ -104,6 +85,6 @@ _printf0_("____________________________________________________________________\n"); - //first run before firing up the control optimization - simulad(&indic,&n,X,&f,G,izs,rzs,(void*)femmodel); + //run post simular phase, to fire up the control optimization + postsimulad(&indic,&n,X,&f,G,izs,rzs,(void*)femmodelad); double f1=f; @@ -125,18 +106,4 @@ } - /*Constrain solution vector*/ - IssmDouble *XL = NULL; - IssmDouble *XU = NULL; - GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound",onsid); - GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound",onsid); - for(int i=0;i(scaling_factors[c]); - if(X[index]>XU[index]) X[index]=reCast(XU[index]); - if(X[index](XL[index]); - } - } - /*Save results:*/ femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,1,0.0)); @@ -147,7 +114,6 @@ xDelete(X); xDelete(dz); - xDelete(scaling_factors); }/*}}}*/ -void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/ +FemModel* presimulad(int* pintn, double** pX, FemModel* femmodel){ /*{{{*/ /*Intermediaries:*/ @@ -157,37 +123,11 @@ char* toolkitsfilename=NULL; char* lockfilename=NULL; - IssmPDouble* G2=NULL; - bool onsid=true; - IssmDouble *XL = NULL; - IssmDouble *XU = NULL; int solution_type; - FemModel *femmodel = NULL; - FemModel *femmodelad = NULL; IssmDouble pfd; + IssmDouble* Xd=NULL; + int intn; + IssmPDouble* X=NULL; int i; - /*Recover Femmodel*/ - femmodel = (FemModel*)dzs; - - /*Recover number of cost functions responses*/ - int num_responses,num_controls,numberofvertices; - IssmDouble* scaling_factors = NULL; - femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); - numberofvertices=femmodel->vertices->NumberOfVertices(); - - /*Constrain input vector*/ - GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound",onsid); - GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound",onsid); - for(int i=0;i(scaling_factors[c]); - if(X[index]>reCast(XU[index])) X[index]=reCast(XU[index]); - if(X[index](XL[index])) X[index]=reCast(XL[index]); - } - } - /*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient @@ -199,7 +139,45 @@ femmodel->parameters->FindParam(&lockfilename,LockFileNameEnum); - femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, femmodel->comm, femmodel->solution_type,X); - femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it - + femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, femmodel->comm, femmodel->solution_type,NULL); + + + /*Get initial guess:*/ + femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum); + X=xNew(intn); for(i=0;i(Xd[i]); + + xDelete(rootpath); + xDelete(inputfilename); + xDelete(outputfilename); + xDelete(toolkitsfilename); + xDelete(lockfilename); + xDelete(Xd); + + *pintn=intn; + *pX=X; + + return femmodel; + +} /*}}}*/ +void postsimulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/ + + /*Intermediaries:*/ + char* rootpath=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* toolkitsfilename=NULL; + char* lockfilename=NULL; + IssmPDouble* G2=NULL; + int solution_type; + FemModel *femmodel = NULL; + IssmDouble pfd; + int i; + + /*Recover Femmodel*/ + femmodel = (FemModel*)dzs; + + /*Recover number of cost functions responses*/ + int num_responses; + femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + /*Recover some parameters*/ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); @@ -236,14 +214,5 @@ /*Constrain X and G*/ IssmDouble Gnorm = 0.; - for(int i=0;i=XU[index]) G[index]=0.; - if(X[index]<=XL[index]) G[index]=0.; - G[index] = G[index]*reCast(scaling_factors[c]); - X[index] = X[index]/reCast(scaling_factors[c]); - Gnorm += G[index]*G[index]; - } - } + for(long i=0;i<*n;i++) Gnorm += G[i]*G[i]; Gnorm = sqrt(Gnorm); @@ -255,8 +224,92 @@ /*Clean-up and return*/ xDelete(Jlist); - xDelete(XU); - xDelete(XL); xDelete(G2); - xDelete(scaling_factors); + + xDelete(rootpath); + xDelete(inputfilename); + xDelete(outputfilename); + xDelete(toolkitsfilename); + xDelete(lockfilename); + +} /*}}}*/ +void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/ + + /*Intermediaries:*/ + char* rootpath=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* toolkitsfilename=NULL; + char* lockfilename=NULL; + IssmPDouble* G2=NULL; + int solution_type; + FemModel *femmodel = NULL; + FemModel *femmodelad = NULL; + IssmDouble pfd; + int i; + + /*Recover Femmodel*/ + femmodel = (FemModel*)dzs; + + /*Recover number of cost functions responses*/ + int num_responses; + femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + + /*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate + *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient + in one run of the solution core. So first recover the filenames required for the FemModel constructor, then call a new ad tailored constructor:*/ + femmodel->parameters->FindParam(&rootpath,RootPathEnum); + femmodel->parameters->FindParam(&inputfilename,InputFileNameEnum); + femmodel->parameters->FindParam(&outputfilename,OutputFileNameEnum); + femmodel->parameters->FindParam(&toolkitsfilename,ToolkitsFileNameEnum); + femmodel->parameters->FindParam(&lockfilename,LockFileNameEnum); + + femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, femmodel->comm, femmodel->solution_type,X); + femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it + + /*Recover some parameters*/ + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + /*Compute solution:*/ + void (*solutioncore)(FemModel*)=NULL; + CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + solutioncore(femmodel); + + /*Compute objective function*/ + IssmDouble* Jlist = NULL; + femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast(pfd); + _printf0_("f(x) = "<* gradient=(GenericExternalResult*)femmodel->results->FindResult(AutodiffJacobianEnum); _assert_(gradient); + G_temp=gradient->GetValues(); + /*copy onto G2, to avoid a leak: */ + G2=xNew(*n); + xMemCpy(G2,G_temp,*n); + } + else G2=xNew(*n); + + /*MPI broadcast results:*/ + ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + + /*Send gradient to m1qn3 core: */ + for(long i=0;i<*n;i++) G[i] = G2[i]; + + /*Recover Gnorm: */ + IssmDouble Gnorm = 0.; + for(int i=0;i<*n;i++) Gnorm += G[i]*G[i]; + Gnorm = sqrt(Gnorm); + + /*Print info*/ + _printf0_(" "<(Jlist); + xDelete(G2); xDelete(rootpath);