Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 45) +++ /issm/trunk/src/c/Makefile.am (revision 46) @@ -71,8 +71,8 @@ ./shared/Numerics/numerics.h\ ./shared/Numerics/GaussPoints.h\ + ./shared/Numerics/GaussPoints.cpp\ ./shared/Numerics/BrentSearch.cpp\ ./shared/Numerics/OptFunc.cpp\ ./shared/Numerics/extrema.cpp\ - ./shared/Numerics/GaussPoints.cpp\ ./shared/Exceptions/exceptions.h\ ./shared/Exceptions/Exceptions.cpp\ @@ -419,5 +419,9 @@ ./parallel/CreateFemModel.cpp\ ./parallel/OutputDiagnostic.cpp\ - ./parallel/WriteLockFile.cpp + ./parallel/WriteLockFile.cpp\ + ./parallel/control.cpp\ + ./parallel/OutputControl.cpp\ + ./parallel/objectivefunctionC.cpp\ + ./parallel/GradJCompute.cpp libpISSM_a_CXXFLAGS = -fPIC -D_PARALLEL_ -D_C_ @@ -426,6 +430,6 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe -#control.exe thermalsteady.exe +bin_PROGRAMS = diagnostic.exe control.exe +#thermalsteady.exe endif Index: /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 45) +++ /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 46) @@ -64,8 +64,6 @@ double tria_q; int tria_shelf; - double tria_fit; double tria_meanvel;/*!scaling ratio for velocities*/ double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ - int tria_acceleration; /*matice constructor input: */ @@ -90,8 +88,6 @@ int penta_onbed; int penta_onsurface; - double penta_fit; double penta_meanvel;/*!scaling ratio for velocities*/ double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ - int penta_acceleration; int penta_collapse; double penta_melting[6]; @@ -151,4 +147,6 @@ int grid_id; + /*Get analysis_type: */ + analysis_type=AnalysisTypeAsEnum(model->analysis_type); /*First create the elements, nodes and material properties: */ @@ -156,5 +154,4 @@ nodes = new DataSet(NodesEnum()); materials = new DataSet(MaterialsEnum()); - /*Width of elements: */ @@ -226,8 +223,6 @@ ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); - ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat"); ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat"); ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat"); - for (i=0;inumberofelements;i++){ @@ -275,14 +270,9 @@ tria_shelf=(int)*(model->elementoniceshelf+i); - /*type of fitting? : */ - tria_fit=model->fit[0]; tria_meanvel=model->meanvel; tria_epsvel=model->epsvel; - /*diverse:*/ - tria_acceleration=0; - /*Create tria element using its constructor:*/ - tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_fit, tria_meanvel, tria_epsvel, tria_acceleration); + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel); /*Add tria element to elements dataset: */ @@ -333,5 +323,4 @@ xfree((void**)&model->q); xfree((void**)&model->elementoniceshelf); - xfree((void**)&model->fit); xfree((void**)&model->B); xfree((void**)&model->n);*/ @@ -348,5 +337,4 @@ ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); - ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat"); ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat"); ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat"); @@ -384,10 +372,8 @@ /*diverse: */ penta_shelf=(int)*(model->elementoniceshelf+i); - penta_fit=(int)model->fit[0]; penta_onbed=(int)*(model->elementonbed+i); penta_onsurface=(int)*(model->elementonsurface+i); penta_meanvel=model->meanvel; penta_epsvel=model->epsvel; - penta_acceleration=0; if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base. @@ -401,6 +387,6 @@ /*Create Penta using its constructor:*/ penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, - penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_fit,penta_meanvel,penta_epsvel, - penta_acceleration,penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff, + penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel, + penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff, penta_thermal_steadystate); @@ -446,5 +432,5 @@ /*Free data: */ - /*xfree((void**)&model->elements); + xfree((void**)&model->elements); xfree((void**)&model->thickness); xfree((void**)&model->surface); @@ -454,10 +440,9 @@ xfree((void**)&model->q); xfree((void**)&model->elementoniceshelf); - xfree((void**)&model->fit); xfree((void**)&model->elementonbed); xfree((void**)&model->elementonsurface); xfree((void**)&model->elements_type); xfree((void**)&model->n); - xfree((void**)&model->B);*/ + xfree((void**)&model->B); } //if (strcmp(meshtype,"2d")==0) @@ -555,7 +540,5 @@ ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat"); - /*Get analysis_type: */ - analysis_type=AnalysisTypeAsEnum(model->analysis_type); - + /*Get number of dofs per node: */ DistributeNumDofs(&node_numdofs,analysis_type); Index: /issm/trunk/src/c/objects/OptArgs.h =================================================================== --- /issm/trunk/src/c/objects/OptArgs.h (revision 45) +++ /issm/trunk/src/c/objects/OptArgs.h (revision 46) @@ -21,6 +21,15 @@ }; #else + +#include "./FemModel.h" +#include "./ParameterInputs.h" + struct OptArgs{ - + FemModel* femmodel; + double* p_g; + double* u_g_obs; + double* grad_g; + ParameterInputs* inputs; + int n; }; #endif Index: /issm/trunk/src/c/objects/OptPars.h =================================================================== --- /issm/trunk/src/c/objects/OptPars.h (revision 45) +++ /issm/trunk/src/c/objects/OptPars.h (revision 46) @@ -16,5 +16,2 @@ #endif - - - Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 45) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 46) @@ -19,6 +19,6 @@ } Penta::Penta( int penta_id, int penta_mid, int penta_mparid, int penta_node_ids[6], double penta_h[6], double penta_s[6], double penta_b[6], double penta_k[6], int penta_friction_type, - double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_fit, double penta_meanvel,double penta_epsvel, - int penta_acceleration, int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], + double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel, + int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], int penta_artdiff, int penta_thermal_steadystate){ @@ -46,8 +46,6 @@ onbed = penta_onbed; onsurface = penta_onsurface; - fit = penta_fit; meanvel = penta_meanvel; epsvel = penta_epsvel; - acceleration = penta_acceleration; collapse = penta_collapse; artdiff = penta_artdiff; @@ -79,8 +77,6 @@ printf(" onbed: %i\n",onbed); printf(" onsurface: %i\n",onsurface); - printf(" fit: %g\n",fit); printf(" meanvel: %g\n",meanvel); printf(" epsvel: %g\n",epsvel); - printf(" acceleration: %i\n",acceleration); printf(" collapse: %i\n",collapse); @@ -122,8 +118,6 @@ memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed); memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); - memcpy(marshalled_dataset,&fit,sizeof(fit));marshalled_dataset+=sizeof(fit); memcpy(marshalled_dataset,&meanvel,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel); memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel); - memcpy(marshalled_dataset,&acceleration,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration); memcpy(marshalled_dataset,&collapse,sizeof(collapse));marshalled_dataset+=sizeof(collapse); memcpy(marshalled_dataset,&melting,sizeof(melting));marshalled_dataset+=sizeof(melting); @@ -139,5 +133,5 @@ int Penta::MarshallSize(){ - return sizeof(id)+sizeof(mid)+sizeof(mparid)+sizeof(node_ids)+sizeof(h)+sizeof(s)+sizeof(b)+sizeof(k)+sizeof(friction_type)+sizeof(p)+sizeof(q)+sizeof(shelf)+sizeof(onbed)+sizeof(onsurface)+sizeof(fit)+sizeof(meanvel)+sizeof(epsvel)+sizeof(acceleration)+sizeof(collapse)+sizeof(melting)+sizeof(accumulation)+sizeof(geothermalflux)+sizeof(artdiff)+sizeof(thermal_steadystate) +sizeof(int); //sizeof(int) for enum type + return sizeof(id)+sizeof(mid)+sizeof(mparid)+sizeof(node_ids)+sizeof(h)+sizeof(s)+sizeof(b)+sizeof(k)+sizeof(friction_type)+sizeof(p)+sizeof(q)+sizeof(shelf)+sizeof(onbed)+sizeof(onsurface)+sizeof(meanvel)+sizeof(epsvel)+sizeof(collapse)+sizeof(melting)+sizeof(accumulation)+sizeof(geothermalflux)+sizeof(artdiff)+sizeof(thermal_steadystate) +sizeof(int); //sizeof(int) for enum type } @@ -170,8 +164,6 @@ memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed); memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); - memcpy(&fit,marshalled_dataset,sizeof(fit));marshalled_dataset+=sizeof(fit); memcpy(&meanvel,marshalled_dataset,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel); memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel); - memcpy(&acceleration,marshalled_dataset,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration); memcpy(&collapse,marshalled_dataset,sizeof(collapse));marshalled_dataset+=sizeof(collapse); memcpy(&melting,marshalled_dataset,sizeof(melting));marshalled_dataset+=sizeof(melting); Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 45) +++ /issm/trunk/src/c/objects/Penta.h (revision 46) @@ -42,8 +42,6 @@ int onbed; int onsurface; - double fit; double meanvel;/*!scaling ratio for velocities*/ double epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ - int acceleration; int collapse; double melting[6]; @@ -58,6 +56,6 @@ Penta(); Penta( int id, int mid, int mparid, int g[6], double h[6], double s[6], double b[6], double k[6], int friction_type, - double p, double q, int shelf, int onbed, int onsurface, double fit, double meanvel,double epsvel, - int acceleration, int collapse, double melting[6], double accumulation[6], double geothermalflux[6], + double p, double q, int shelf, int onbed, int onsurface, double meanvel,double epsvel, + int collapse, double melting[6], double accumulation[6], double geothermalflux[6], int artdiff, int thermal_steadystate); ~Penta(); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 45) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 46) @@ -33,5 +33,5 @@ Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3], - int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_fit,double tria_meanvel,double tria_epsvel,int tria_acceleration){ + int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel){ int i; @@ -57,8 +57,6 @@ q=tria_q; shelf=tria_shelf; - fit=tria_fit; meanvel=tria_meanvel; epsvel=tria_epsvel; - acceleration=tria_acceleration; onbed=1; @@ -87,8 +85,6 @@ printf(" q: %g\n",q); printf(" shelf: %i\n",shelf); - printf(" fit: %g\n",fit); printf(" meanvel: %g\n",meanvel); printf(" epsvel: %g\n",epsvel); - printf(" acceleration: %i\n",acceleration); printf(" onbed: %i\n",onbed); printf(" nodes: \n"); @@ -136,8 +132,6 @@ memcpy(marshalled_dataset,&q,sizeof(q));marshalled_dataset+=sizeof(q); memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf); - memcpy(marshalled_dataset,&fit,sizeof(fit));marshalled_dataset+=sizeof(fit); memcpy(marshalled_dataset,&meanvel,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel); memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel); - memcpy(marshalled_dataset,&acceleration,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration); *pmarshalled_dataset=marshalled_dataset; @@ -165,8 +159,6 @@ +sizeof(q) +sizeof(shelf) - +sizeof(fit) +sizeof(meanvel) +sizeof(epsvel) - +sizeof(acceleration) +sizeof(int); //sizeof(int) for enum type } @@ -206,8 +198,6 @@ memcpy(&q,marshalled_dataset,sizeof(q));marshalled_dataset+=sizeof(q); memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf); - memcpy(&fit,marshalled_dataset,sizeof(fit));marshalled_dataset+=sizeof(fit); memcpy(&meanvel,marshalled_dataset,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel); memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel); - memcpy(&acceleration,marshalled_dataset,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration); /*nodes and materials are not pointing to correct objects anymore:*/ @@ -338,10 +328,4 @@ double MOUNTAINKEXPONENT=10; - - /*First check that acceleration is not on: */ - if (acceleration){ - /*do nothing: */ - return; - } /*recover extra inputs from users, at current convergence iteration: */ @@ -650,9 +634,4 @@ double thickness; - /*First check that acceleration is not on: */ - if (acceleration){ - return; - } - /* Get node coordinates and dof list: */ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); @@ -1313,4 +1292,5 @@ fit=*fit_param; } + else throw ErrorException(__FUNCT__," missing fit input parameter"); /*Initialize velocities: */ @@ -1844,4 +1824,5 @@ fit=*fit_param; } + else ErrorException(__FUNCT__," missing fit input parameter"); /*Initialize velocities: */ Index: /issm/trunk/src/c/objects/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Tria.h (revision 45) +++ /issm/trunk/src/c/objects/Tria.h (revision 46) @@ -40,8 +40,6 @@ double q; int shelf; - double fit; double meanvel;/*!scaling ratio for velocities*/ double epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ - int acceleration; int onbed; @@ -50,5 +48,5 @@ Tria(); Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3], - int friction_type,double p,double q,int shelf,double fit,double meanvel,double epsvel,int acceleration); + int friction_type,double p,double q,int shelf,double meanvel,double epsvel); ~Tria(); Index: /issm/trunk/src/c/parallel/CreateFemModel.cpp =================================================================== --- /issm/trunk/src/c/parallel/CreateFemModel.cpp (revision 45) +++ /issm/trunk/src/c/parallel/CreateFemModel.cpp (revision 46) @@ -68,4 +68,7 @@ ConfigureObjectsx(elements, loads, nodes, materials); + _printf_(" process parameters:\n"); + ProcessParamsx( parameters, partition); + _printf_(" free ressources:\n"); DeleteModel(&model); @@ -86,4 +89,5 @@ femmodel->ys=ys; femmodel->ys0=ys0; + } Index: sm/trunk/src/c/parallel/GradJCompute.c =================================================================== --- /issm/trunk/src/c/parallel/GradJCompute.c (revision 45) +++ (revision ) @@ -1,89 +1,0 @@ -/* - * GradJCompute.c: - */ - -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "GradJCompute" -#undef CLEANUP -#define CLEANUP GradJComputeLocalCleanup(); - -void GradJComputeLocalCleanup(void); - -int GradJCompute(ParameterInputs* inputs,FemModel* femmodel){ - - /*Error management: */ - int noerr=1; - int i; - int dummy; - - Vec* u_g=NULL; - Vec* du_g=NULL; - Vec* du_f=NULL; - Vec* lambda_f=NULL; - Vec* lambda_g=NULL; - double* lambda_g_double=NULL; - double* u_g_double=NULL; - Vec* gradj_g=NULL; - Mat* K_ff0=NULL; - Mat* K_fs0=NULL; - - //Recover uset: */ - uset=femmodel->uset; //external variable - - //Recover solution for this stiffness and right hand side: - cielodiagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,inputs,femmodel); - - //Buid Du, difference between observed velocity and model velocity. - VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); - - Dux(&du_g,femmodel->bgpdt, femmodel->bgpdtb, femmodel->est, femmodel->ept, femmodel->mpt,femmodel->lst, u_g_double,femmodel->workspaceparams->u_g_obs,inputs,femmodel->batchparams->analysis_type_enum); - - - //Reduce adjoint load from g-set to f-set - Reducerightside(&du_f, du_g, femmodel->G_mn, K_fs0, femmodel->y_s0, femmodel->flag_y_s0, uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, uset->msize,uset->ssize,uset->fsize); - VecFree(&du_g);MatFree(&K_fs0); - - - //Solve for adjoint vector: - lambda_f=xmalloc(sizeof(Vec)); - Solverx(lambda_f,&dummy,K_ff0,dummy,dummy,du_f,dummy,NULL,0,femmodel->batchparams->solverstring); - VecFree(&du_f); - MatFree(&K_ff0); - - - //Merge back to g set - Mergesolvec( &lambda_g, lambda_f, femmodel->G_mn, femmodel->y_s0, uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, uset->pv_s, uset->ssize, uset->gsize,uset->msize,uset->nsize); - VecFree(&lambda_f); - - //Compute gradJ for each parameter - VecToMPISerial(&lambda_g_double,lambda_g);VecFree(&lambda_g); - - for(i=0;iworkspaceparams->num_control_parameters;i++){ - char* control_type=femmodel->workspaceparams->control_types[i]; - - Gradjx( &gradj_g, femmodel->bgpdt, femmodel->bgpdtb, femmodel->est, femmodel->ept, femmodel->mpt,femmodel->lst, u_g_double,lambda_g_double,inputs,femmodel->batchparams->analysis_type_enum,control_type); - - WorkspaceParamsSetParameterGradient(femmodel->workspaceparams,gradj_g,control_type); - } - - /*free ressources: */ - xfree((void**)&lambda_g_double); - xfree((void**)&u_g_double); - - - EXIT(noerr); -} - -void GradJComputeLocalCleanup(void){ - return; -} -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: /issm/trunk/src/c/parallel/GradJCompute.cpp =================================================================== --- /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 46) +++ /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 46) @@ -0,0 +1,80 @@ +/*!\file: GradJCompute.cpp + * \brief compute inverse method gradient direction. + */ + +#include "../issm.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "GradJCompute" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel,double* u_g_obs){ + + + /*intermediary: */ + int analysis_type; + char* solver_string=NULL; + char* control_type=NULL; + + Vec u_g=NULL; + double* u_g_double=NULL; + + Vec du_g=NULL; + Vec du_f=NULL; + + Vec lambda_f=NULL; + Vec lambda_g=NULL; + double* lambda_g_double=NULL; + + Mat K_ff0=NULL; + Mat K_fs0=NULL; + + /*output: */ + Vec grad_g=NULL; + + + /*some parameters:*/ + femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type"); + femmodel->parameters->FindParam((void*)&solver_string,"solver_string"); + femmodel->parameters->FindParam((void*)&control_type,"control_type"); + + _printf_("%s\n"," Recover solution for this stiffness and right hand side: \n"); + diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,inputs,femmodel); + VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); + + _printf_("%s\n"," Buid Du, difference between observed velocity and model velocity:\n"); + Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,u_g_double,u_g_obs, inputs,analysis_type); + + _printf_("%s\n"," Reduce adjoint load from g-set to f-set:\n"); + Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys0, femmodel->nodesets); + VecFree(&du_g);MatFree(&K_fs0); + + _printf_("%s\n"," Solve for adjoint vector: \n"); + Solverx(&lambda_f, K_ff0, du_f, NULL, solver_string); + VecFree(&du_f); MatFree(&K_ff0); + + _printf_("%s\n"," Merge back to g set:\n"); + Mergesolutionfromftogx(&lambda_g, lambda_f,femmodel->Gmn,femmodel->ys0,femmodel->nodesets); + VecFree(&lambda_f); + + _printf_("%s\n"," Compute gradJ:\n"); + VecToMPISerial(&lambda_g_double,lambda_g);VecFree(&lambda_g); + + Gradjx( &grad_g, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, + u_g_double,lambda_g_double, inputs,analysis_type,control_type); + + /*Free ressources:*/ + xfree((void**)&solver_string); + xfree((void**)&control_type); + xfree((void**)&u_g_double); + xfree((void**)&lambda_g_double); + + return grad_g; + +} Index: sm/trunk/src/c/parallel/OutputControl.c =================================================================== --- /issm/trunk/src/c/parallel/OutputControl.c (revision 45) +++ (revision ) @@ -1,99 +1,0 @@ - -/* - OutputControl.c: output model results for control methods. -*/ -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" - -#undef __FUNCT__ -#define __FUNCT__ "OutputControl" - -#undef CLEANUP -#define CLEANUP OutputControlLocalCleanup(); - -void OutputControlLocalCleanup(void); - -int OutputControl(WorkspaceParams* workspaceparams,BatchParams* batchparams,Vec* u_g,Vec* partition, char* filename,char* analysis_type){ - - /* error handling: */ - int i; - int noerr=1; - - /* output: */ - FILE* fid=NULL; - - /* standard output: */ - double* serial_partition=NULL; - int serial_partition_rows; - int dummy=1; - int analysis_size; - - double* serial_u_g=NULL; - int serial_u_g_rows; - - double* parameter=NULL; - int parameter_rows; - - /*serialize outputs: */ - VecShift(partition,1.0); //matlab indexing - VecGetSize(*partition,&serial_partition_rows); - noerr=VecToMPISerial(&serial_partition,partition);TESTEXIT(noerr); - - VecGetSize(*u_g,&serial_u_g_rows); - noerr=VecToMPISerial(&serial_u_g,u_g);TESTEXIT(noerr); - - /* Open output file to write raw binary data: */ - if(my_rank==0){ - fid=fopen(filename,"wb"); - if(fid==NULL){ - _printf_("%s%s%s%s\n",__FUNCT__," error message: could not open file ",filename," for binary writing."); - noerr=0; goto cleanup_and_return; - } - - /*Write solution type: */ - WriteDataToDisk(analysis_type,NULL,NULL,"String",fid); - - /*Write uset.gsize: */ - WriteDataToDisk(&uset->gsize,NULL,NULL,"Integer",fid); - - /*Write partition: */ - WriteDataToDisk(serial_partition,&serial_partition_rows,&dummy,"Mat",fid); - - /*Write misfit J vector: */ - WriteDataToDisk(&workspaceparams->num_control_parameters,NULL,NULL,"Integer",fid); - WriteDataToDisk(&workspaceparams->nsteps,NULL,NULL,"Integer",fid); - WriteDataToDisk(workspaceparams->J,&workspaceparams->nsteps,&dummy,"Mat",fid); - - /*Go through all control parameters, and write optimized parameters: */ - for(i=0;inum_control_parameters;i++){ - char* control_type=workspaceparams->control_types[i]; - parameter=WorkspaceParamsGetParameter(workspaceparams,control_type); - parameter_rows=uset->gsize; - WriteDataToDisk(parameter,¶meter_rows,&dummy,"Mat",fid); - } - - /*Write solution: */ - WriteDataToDisk(serial_u_g,&serial_u_g_rows,&dummy,"Mat",fid); - - /*Close file: */ - if(fclose(fid)!=0){ - _printf_("%s%s%s\n",__FUNCT__," error message: could not close file ",filename); - noerr=0; goto cleanup_and_return; - } - } - - cleanup_and_return: - TESTEXIT(noerr); - - EXIT(noerr); -} -void OutputControlLocalCleanup(void){ - return; -} - - -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: /issm/trunk/src/c/parallel/OutputControl.cpp =================================================================== --- /issm/trunk/src/c/parallel/OutputControl.cpp (revision 46) +++ /issm/trunk/src/c/parallel/OutputControl.cpp (revision 46) @@ -0,0 +1,63 @@ + +/* + OutputControl.c: output model results for control solution. +*/ +#undef __FUNCT__ +#define __FUNCT__ "OutputControl" + +#include "../toolkits/toolkits.h" +#include "../shared/shared.h" +#include "../io/io.h" +#include "../objects/objects.h" + +void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets){ + + int i; + extern int my_rank; + + /* output: */ + FILE* fid=NULL; + + /* standard output: */ + double* serial_partition=NULL; + + double* serial_u_g=NULL; + int one=1; + int gsize; + + /*serialize outputs: */ + VecShift(partition,1.0); //matlab indexing + VecToMPISerial(&serial_partition,partition); + + VecToMPISerial(&serial_u_g,u_g); + + /* Open output file to write raw binary data: */ + if(my_rank==0){ + fid=fopen(filename,"wb"); + if(fid==NULL)throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",filename," for binary writing.")); + + /*Write solution type: */ + WriteDataToDisk((void*)"control",NULL,NULL,"String",fid); + + /*Write gsize: */ + gsize=nodesets->GetGSize(); + WriteDataToDisk(&gsize,NULL,NULL,"Integer",fid); + + /*Write partition: */ + WriteDataToDisk(serial_partition,&gsize,&one,"Mat",fid); + + /*Write solution to disk: */ + WriteDataToDisk(serial_u_g,&gsize,&one,"Mat",fid); + + /*Write parameter to disk: */ + WriteDataToDisk(p_g,&gsize,&one,"Mat",fid); + + /*Write J to disk: */ + WriteDataToDisk(&nsteps,NULL,NULL,"Integer",fid); + WriteDataToDisk(J,&nsteps,&one,"Mat",fid); + + /*Close file: */ + if(fclose(fid)!=0)throw ErrorException(__FUNCT__,exprintf("%s%s","could not close file ",filename)); + } + +} Index: /issm/trunk/src/c/parallel/OutputDiagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/OutputDiagnostic.cpp (revision 45) +++ /issm/trunk/src/c/parallel/OutputDiagnostic.cpp (revision 46) @@ -13,7 +13,5 @@ void OutputDiagnostic(Vec u_g,Vec partition,char* filename,NodeSets* nodesets,char* analysis_type){ - /* error handling: */ int i; - extern int my_rank; @@ -23,18 +21,14 @@ /* standard output: */ double* serial_partition=NULL; - int serial_partition_rows; int analysis_size; double* serial_u_g=NULL; - int serial_u_g_rows; - int dummy=1; + int one=1; int gsize; /*serialize outputs: */ VecShift(partition,1.0); //matlab indexing - VecGetSize(partition,&serial_partition_rows); VecToMPISerial(&serial_partition,partition); - VecGetSize(u_g,&serial_u_g_rows); VecToMPISerial(&serial_u_g,u_g); @@ -52,8 +46,8 @@ /*Write partition: */ - WriteDataToDisk(serial_partition,&serial_partition_rows,&dummy,"Mat",fid); + WriteDataToDisk(serial_partition,&one,&one,"Mat",fid); /*Write solution to disk: */ - WriteDataToDisk(serial_u_g,&serial_u_g_rows,&dummy,"Mat",fid); + WriteDataToDisk(serial_u_g,&gsize,&one,"Mat",fid); /*Close file: */ Index: /issm/trunk/src/c/parallel/control.cpp =================================================================== --- /issm/trunk/src/c/parallel/control.cpp (revision 45) +++ /issm/trunk/src/c/parallel/control.cpp (revision 46) @@ -1,17 +1,23 @@ -/* - * cielocontrol.c: - */ +/*!\file: control.cpp + * \brief: control solution + */ -#include "../include/cielo.h" -#include "../modules.h" +#include "../issm.h" #include "./parallel.h" #undef __FUNCT__ -#define __FUNCT__ "cielocontrol" +#define __FUNCT__ "control" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif int main(int argc,char* *argv){ - - int i,n; - + + int n; + int i; + /*I/O: */ FILE* fid=NULL; @@ -21,148 +27,149 @@ char* analysis_type="control"; - /*Finite element model: */ - #if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + /*Intermediary: */ FemModel femmodel; + Vec u_g=NULL; + ParameterInputs* inputs=NULL; + int waitonlock=0; + double search_scalar; + char* control_type=NULL; + int gsize; + double* fit=NULL; + double* optscal=NULL; + double* u_g_obs=NULL; + int* maxiter=NULL; + double tolx; + double* p_g=NULL; + Vec grad_g=NULL; + Vec new_grad_g=NULL; + Vec grad_g_old=NULL; + double* grad_g_double=NULL; + double mincontrolconstraint; + double maxcontrolconstraint; + int nsteps; + double* J=NULL; + OptArgs optargs; + OptPars optpars; + + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!"); #endif - /*control : */ - #if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Vec* u_g=NULL; - #endif - double* search_vector=NULL; - int status; - ParameterInputs* inputs=NULL; - int reloop; + PetscInitialize(&argc,&argv,(char *)0,""); + + /*Size and rank: */ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); + + _printf_("recover , input file name and output file name:\n"); + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Open handle to data on disk: */ + fid=fopen(inputfilename,"rb"); + if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading")); - #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) - _printf_("%s%s\n",__FUNCT__," error message: parallel executable was compiled without support of parallel libraries!"); - return 1; - #else + _printf_("read and create finite element model:\n"); + CreateFemModel(&femmodel,fid,analysis_type); + + /*Recover parameters used throughout the solution:*/ + femmodel.parameters->FindParam((void*)&nsteps,"nsteps"); + femmodel.parameters->FindParam((void*)&control_type,"control_type"); + femmodel.parameters->FindParam((void*)&fit,"fit"); + femmodel.parameters->FindParam((void*)&optscal,"optscal"); + femmodel.parameters->FindParam((void*)&maxiter,"maxiter"); + femmodel.parameters->FindParam((void*)&tolx,"tolx"); + femmodel.parameters->FindParam((void*)&waitonlock,"waitonlock"); + femmodel.parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint"); + femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint"); + gsize=femmodel.nodes->NumberOfDofs(); - /*Initialize MPI environment: */ - PetscInitialize(&argc,&argv,(char *)0,""); + femmodel.parameters->FindParam((void*)&p_g,"p_g"); + femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs"); - /*Size and rank: */ - MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); - MPI_Comm_size(MPI_COMM_WORLD,&num_procs); - /*Some checks on size of cluster*/ - if (num_procs<=1){ - _printf_("\nSize of MPI COMM WORLD is 1, needs to be at least 2. Include more nodes\n"); - PetscFinalize(); - return 0; + /*Start looping: */ + for(n=0;nnum_control_parameters>1){ - _printf_("%s%s\n",__FUNCT__," error message: multiple control parameters not implemented yet!"); - PetscFinalize(); - return 1; - } - - /*Initialize search vector: */ - search_vector=xmalloc(femmodel.workspaceparams->num_control_parameters*sizeof(double)); - - /*Start looping: */ - for(n=0;nnsteps;n++){ - - _printf_("\n%s%i%s%i\n"," control method step ",(n+1),"/",femmodel.workspaceparams->nsteps); - - //initialize inputs, ie parameters on which we invert. - DeleteParameterInputs(&inputs); inputs=NewParameterInputs(); - - for(i=0;inum_control_parameters;i++){ - char* control_type=femmodel.workspaceparams->control_types[i]; - ParameterInputsAddFromMat(inputs,WorkspaceParamsGetParameter(femmodel.workspaceparams,control_type),femmodel.workspaceparams->gsize,control_type); - } - ParameterInputsAddFromDouble(inputs,femmodel.workspaceparams->fit[n],"fit"); - - //update velocity: - _printf_("%s\n"," updating velocity..."); - VecFree(&u_g); cielodiagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel); - _printf_("%s\n"," done."); - - //call gradJ module - _printf_("%s\n"," computing gradJ..."); - GradJCompute(inputs,&femmodel); - _printf_("%s\n"," done."); - - //normalize directions - _printf_("%s\n"," normalizing directions..."); - GradJOrth(femmodel.workspaceparams); - _printf_("%s\n"," done."); + /*Close MPI libraries: */ + PetscFinalize(); - //search along the direction for the parameter vector that minimizes the misfit J - _printf_("%s\n"," optimizing..."); - status=GradJSearch(search_vector,&femmodel,n); - _printf_("%s\n"," done."); - - _printf_("%s%g\n"," misfit: ",femmodel.workspaceparams->J[n]); - - //Check the GradJSearch actually resulting in improving the misfit. Otherwise, reloop. - reloop=GradJCheck(femmodel.workspaceparams,n,status); - if (reloop){ - printf("reloop\n",reloop); - n=n-1; - continue; - } - - //update parameters with new optimized values - _printf_("%s\n"," updating parameters..."); - ParameterUpdate(search_vector,n,femmodel.workspaceparams,femmodel.batchparams); - _printf_("%s\n"," done."); - - //temporary saving - if ((n%1)==0){ - _printf_("%s\n"," saving temporary results..."); - OutputControl(femmodel.workspaceparams,femmodel.batchparams,u_g,femmodel.tpartition,outputfilename,"control"); - _printf_("%s\n"," done."); - } - - } - - /*Write results to disk: */ - _printf_("%s\n"," saving final results..."); - OutputControl(femmodel.workspaceparams,femmodel.batchparams,u_g,femmodel.tpartition,outputfilename,"control"); - _printf_("%s\n"," done."); - - /*Write lock file if requested: */ - if (femmodel.batchparams->waitonlock){ - WriteLockFile(lockname); - } - - cleanup_and_return: - _printf_("exiting.\n"); - - /*Synchronize everyone before exiting: */ - MPI_Barrier(MPI_COMM_WORLD); - - /*Close MPI libraries: */ - PetscFinalize(); - - return 0; //unix success return; - #endif + /*end module: */ + MODULEEND(); + + return 0; //unix success return; } Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 45) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 46) @@ -59,5 +59,5 @@ _printf_("call computational core:\n"); - diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,femmodel); + diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel); _printf_("write results to disk:\n"); Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 45) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 46) @@ -11,5 +11,5 @@ #include "../issm.h" -void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, ParameterInputs* inputs,FemModel fem){ +void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, ParameterInputs* inputs,FemModel* fem){ @@ -45,15 +45,15 @@ /*Recover parameters: */ kflag=1; pflag=1; - fem.parameters->FindParam((void*)&connectivity,"connectivity"); - fem.parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode"); - fem.parameters->FindParam((void*)&solver_string,"solverstring"); - fem.parameters->FindParam((void*)&eps_rel,"eps_rel"); - fem.parameters->FindParam((void*)&debug,"debug"); + fem->parameters->FindParam((void*)&connectivity,"connectivity"); + fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode"); + fem->parameters->FindParam((void*)&solver_string,"solverstring"); + fem->parameters->FindParam((void*)&eps_rel,"eps_rel"); + fem->parameters->FindParam((void*)&debug,"debug"); - fem.parameters->FindParam((void*)&analysis_type_string,"analysis_type"); + fem->parameters->FindParam((void*)&analysis_type_string,"analysis_type"); analysis_type=AnalysisTypeAsEnum(analysis_type_string); /*Copy loads for backup: */ - loads=fem.loads->Copy(); + loads=fem->loads->Copy(); count=1; @@ -71,17 +71,17 @@ /*Update parameters: */ - UpdateFromInputsx(fem.elements,fem.nodes,loads, fem.materials,inputs); + UpdateFromInputsx(fem->elements,fem->nodes,loads, fem->materials,inputs); if (debug) _printf_(" Generating matrices\n"); //*Generate system matrices - SystemMatricesx(&Kgg, &pg,fem.elements,fem.nodes,loads,fem.materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); if (debug) _printf_(" Generating penalty matrices\n"); //*Generate penalty system matrices - PenaltySystemMatricesx(Kgg, pg,fem.elements,fem.nodes,loads,fem.materials,kflag,pflag,inputs,analysis_type); + PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type); if (debug) _printf_(" reducing matrix from g to f set\n"); /*!Reduce matrix from g to f size:*/ - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem.Gmn,fem.nodesets); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); /*Free ressources: */ @@ -90,5 +90,5 @@ if (debug) _printf_(" reducing load from g to f set\n"); /*!Reduce load from g to f size: */ - Reduceloadfromgtofx(&pf, pg, fem.Gmn, Kfs, fem.ys, fem.nodesets); + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); //no need for pg and Kfs anymore @@ -106,5 +106,5 @@ if (debug) _printf_(" merging solution from f to g set\n"); //Merge back to g set - Mergesolutionfromftogx(&ug, uf,fem.Gmn,fem.ys,fem.nodesets); + Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets); //Deal with penalty loads @@ -112,5 +112,5 @@ ParameterInputsAddFromVec(inputs,ug,"velocity"); - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem.elements,fem.nodes,loads,fem.materials,inputs,analysis_type); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type); //Figure out if convergence is reached. @@ -143,7 +143,7 @@ kflag=1; pflag=0; //stiffness generation only - SystemMatricesx(&Kgg, &pg,fem.elements,fem.nodes,fem.loads,fem.materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem.Gmn,fem.nodesets); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg);VecFree(&pg); Index: sm/trunk/src/c/parallel/objectivefunctionC.c =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.c (revision 45) +++ (revision ) @@ -1,113 +1,0 @@ -/* - * objectivefunctionC.c: - */ -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "objectivefunctionC" - -double objectivefunctionC(double* search_vector,double fit,double optscal,FemModel* femmodel,ParameterInputs* inputs){ - - - /*Error management: */ - int noerr=1; - int i; - int j; - double* parameter=NULL; - double* newparameter=NULL; - Vec* vec_gradient=NULL; - double* gradient=NULL; - Vec* u_g=NULL; - double* u_g_double=NULL; - - /*output: */ - double J; - - /*inputs recovered from the FemModel: */ - WorkspaceParams* workspaceparams=NULL; - BatchParams* batchparams=NULL; - DataSet* bgpdt=NULL; - DataSet* bgpdtb=NULL; - DataSet* est=NULL; - DataSet* lst=NULL; - DataSet* ept=NULL; - DataSet* mpt=NULL; - DataSet* geom3=NULL; - Mat* G_mn=NULL; - Vec* y_s=NULL; - int flag_y_s=-1; - int analysis_type=-1; - - /*Recover model inputs: */ - workspaceparams=femmodel->workspaceparams; - batchparams=femmodel->batchparams; - bgpdt=femmodel->bgpdt; - bgpdtb=femmodel->bgpdtb; - est=femmodel->est; - lst=femmodel->lst; - ept=femmodel->ept; - mpt=femmodel->mpt; - geom3=femmodel->geom3; - G_mn=femmodel->G_mn; - y_s=femmodel->y_s; - flag_y_s=femmodel->flag_y_s; - analysis_type=femmodel->batchparams->analysis_type_enum; - - //Go through parameters, and update along gradients, multiplying by search_vector. - for(i=0;inum_control_parameters;i++){ - - - char* control_type=workspaceparams->control_types[i]; - - /*Get parameter, and gradient for the parameter: */ - parameter=WorkspaceParamsGetParameter(workspaceparams,control_type); - vec_gradient=WorkspaceParamsGetParameterGradient(workspaceparams,control_type); - - - /*serialize gradient: */ - VecToMPISerial(&gradient,vec_gradient); - - - /*Ok, for this parameter, we have a direction. Update the parameter along - * the direction, using the search_vector value. Now, because parameter comes - * directly from a workspaceparams pointer, we do not want to modify it. Make a copy - * before updating it.*/ - newparameter=xmalloc(workspaceparams->gsize*sizeof(double)); - memcpy(newparameter,parameter,workspaceparams->gsize*sizeof(double)); - - - for(j=0;j<(int)(workspaceparams->gsize/6);j++){ - newparameter[6*j+0]=parameter[6*j+0]-search_vector[i]*optscal*gradient[6*j+0]; - } - - - /*Add newparameter to inputs: */ - ParameterInputsAddFromMat(inputs,newparameter,workspaceparams->gsize,control_type); - - /*Now, call on each parameter's self constraining routine: */ - WorkspaceParamsConstrain(workspaceparams,control_type); - } - - //Run diagnostic with updated parameters. - cielodiagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,femmodel); - VecToMPISerial(&u_g_double,u_g); - - //Compute misfit for this velocity field. - ParameterInputsAddFromDouble(inputs,fit,"fit"); - Misfitx( &J, femmodel->bgpdt, femmodel->bgpdtb, femmodel->est, femmodel->ept, femmodel->mpt,femmodel->geom3, u_g_double, workspaceparams->u_g_obs,inputs, femmodel->batchparams->analysis_type_enum); - - /*Free ressources: */ - xfree((void**)&gradient); - xfree((void**)&newparameter); - xfree((void**)&u_g_double); - - return J; -} -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 46) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 46) @@ -0,0 +1,91 @@ +/*!\file: objectivefunctionC + * \brief objective function that returns a misfit, for a certain parameter. + */ + +#include "../issm.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "diagnostic" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +double objectivefunctionC(double search_scalar,OptArgs* optargs){ + + int i; + + /*output: */ + double J; + + /*parameters: */ + FemModel* femmodel=NULL; + double* p_g=NULL; + double* u_g_obs=NULL; + double* grad_g=NULL; + ParameterInputs* inputs=NULL; + int n; + + /*intermediary:*/ + int gsize; + double* optscal=NULL; + double* fit=NULL; + double mincontrolconstraint; + double maxcontrolconstraint; + char* control_type=NULL; + double* p_g_copy=NULL; + int analysis_type; + Vec u_g=NULL; + double* u_g_double=NULL; + + + /*Recover parameters: */ + femmodel=optargs->femmodel; + p_g=optargs->p_g; + u_g_obs=optargs->u_g_obs; + grad_g=optargs->grad_g; + inputs=optargs->inputs; + n=optargs->n; + + gsize=femmodel->nodesets->GetGSize(); + femmodel->parameters->FindParam((void*)&optscal,"optscal"); + femmodel->parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint"); + femmodel->parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint"); + femmodel->parameters->FindParam((void*)&control_type,"control_type"); + femmodel->parameters->FindParam((void*)&fit,"fit"); + femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type"); + + /*First copy p_g so we don't modify it: */ + p_g_copy=(double*)xmalloc(gsize*sizeof(double)); + memcpy(p_g_copy,p_g,gsize*sizeof(double)); + + /*First, update p_g using search_scalar: */ + for(i=0;ielements,femmodel->nodes, femmodel->loads, femmodel->materials, + u_g_double,u_g_obs, inputs,analysis_type); + + + /*Free ressources:*/ + xfree((void**)&optscal); + xfree((void**)&control_type); + xfree((void**)&p_g_copy); + xfree((void**)&u_g_double); + + return J; +} Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 45) +++ /issm/trunk/src/c/parallel/parallel.h (revision 46) @@ -8,7 +8,7 @@ #include "../objects/objects.h" -int GradJCompute(ParameterInputs* inputs,FemModel* femmodel); +Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel,double* u_g_obs); -void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, ParameterInputs* inputs,FemModel fem); +void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, ParameterInputs* inputs,FemModel* fem); int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel); @@ -19,5 +19,5 @@ int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*,ParameterInputs*),FemModel* femmodel); -double objectivefunctionC(double* search_vector,double fit,double optscal,FemModel* femmodel,ParameterInputs* inputs); +double objectivefunctionC(double search_scalar,OptArgs* optargs); int GradJSearch(double* search_vector,FemModel* femmodel,int step); @@ -27,5 +27,5 @@ void OutputDiagnostic(Vec u_g,Vec tpartition,char* filename,NodeSets* nodesets,char* analysis_type); int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type); -//int OutputControl(WorkspaceParams* workspaceparams,BatchParams* batchparams,Vec* u_g,Vec* tpartition, char* filename,char* analysis_type); +void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets); void WriteLockFile(char* filename);