Changeset 22612
- Timestamp:
- 03/22/18 15:38:35 (7 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 2 added
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/Makefile.am
r22580 r22612 77 77 ./classes/Misfit.cpp\ 78 78 ./classes/Cfsurfacesquare.cpp\ 79 ./classes/Cfdragcoeffabsgrad.cpp\ 79 80 ./classes/Cfsurfacelogvel.cpp\ 80 81 ./classes/Regionaloutput.cpp\ -
issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp
r22507 r22612 33 33 this->datatime=0.; 34 34 this->timepassedflag = false; 35 this->last_time=0.; 35 36 36 37 } … … 105 106 /*recover time parameters: */ 106 107 femmodel->parameters->FindParam(&time,TimeEnum); 107 108 /*ELLEN UNFINISHED*/ 108 if(time < last_time) timepassedflag = false; 109 last_time = time; 110 109 111 int i; 110 112 IssmDouble J=0.; -
issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.h
r22507 r22612 3 3 */ 4 4 5 #ifndef _CFSURFACE SQUARE_H_6 #define _CFSURFACE SQUARE_H_5 #ifndef _CFSURFACELOGVEL_H_ 6 #define _CFSURFACELOGVEL_H_ 7 7 8 8 /*Headers:*/ … … 21 21 IssmDouble datatime; 22 22 bool timepassedflag; 23 IssmDouble last_time; 23 24 24 25 int lock; // if lock is on, we just return the value stored in "misfit". this is used so we don't compute misfit past the final_time … … 44 45 IssmDouble Cfsurfacelogvel_Calculation(Element* element, int definitionenum); 45 46 }; 46 #endif /* _CFSURFACE SQUARE_H_ */47 #endif /* _CFSURFACELOGVEL_H_ */ -
issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp
r22517 r22612 36 36 this->datatime=0.; 37 37 this->timepassedflag = false; 38 this->last_time = 0.; 38 39 39 40 } … … 108 109 /*}}}*/ 109 110 IssmDouble Cfsurfacesquare::Response(FemModel* femmodel){/*{{{*/ 110 111 111 /*diverse: */ 112 112 IssmDouble time; … … 114 114 /*recover time parameters: */ 115 115 femmodel->parameters->FindParam(&time,TimeEnum); 116 if(time < last_time) timepassedflag = false; 117 last_time = time; 116 118 117 119 int i; 118 120 IssmDouble J=0.; 119 121 IssmDouble J_sum=0.; 122 120 123 if(datatime<=time && !timepassedflag){ 121 124 for(i=0;i<femmodel->elements->Size();i++){ … … 130 133 timepassedflag = true; 131 134 return J; 132 135 } 133 136 else return J; 134 137 } -
issm/trunk-jpl/src/c/classes/Cfsurfacesquare.h
r22507 r22612 3 3 */ 4 4 5 #ifndef _ TIMEMISFIT_H_6 #define _ TIMEMISFIT_H_5 #ifndef _CFSURFACESQUARE_H_ 6 #define _CFSURFACESQUARE_H_ 7 7 8 8 /*Headers:*/ … … 26 26 IssmDouble datatime; 27 27 bool timepassedflag; 28 IssmDouble last_time; 28 29 29 30 int lock; // if lock is on, we just return the value stored in "misfit". this is used so we don't compute misfit past the final_time … … 49 50 IssmDouble Cfsurfacesquare_Calculation(Element* element, int model_enum, int observation_enum, int weights_enum); 50 51 }; 51 #endif /* _ TIMEMISFIT_H_ */52 #endif /* _CFSURFACESQUARE_H_ */ -
issm/trunk-jpl/src/c/classes/classes.h
r22507 r22612 21 21 #include "./Numberedcostfunction.h" 22 22 #include "./Cfsurfacesquare.h" 23 #include "./Cfdragcoeffabsgrad.h" 23 24 #include "./Cfsurfacelogvel.h" 24 25 #include "./Masscon.h" -
issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp
r22515 r22612 65 65 /*Recover Arguments*/ 66 66 m1qn3_struct *input_struct = (m1qn3_struct*)dzs; 67 FemModel *femmodel = input_struct->femmodel; 67 68 FemModel* femmodel = input_struct->femmodel; 69 int num_responses,num_controls,numberofvertices,solution_type; 70 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); 71 72 if (solution_type == TransientSolutionEnum){ 73 femmodel = input_struct->femmodel->copy(); 74 } 75 68 76 IssmPDouble *Jlist = input_struct->Jlist; 69 77 int JlistM = input_struct->M; … … 73 81 74 82 /*Recover some parameters*/ 75 int num_responses,num_controls,numberofvertices,solution_type;76 83 IssmDouble* scaling_factors = NULL; 77 84 femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); 78 85 femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); 79 86 femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); 80 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);81 87 numberofvertices=femmodel->vertices->NumberOfVertices(); 82 88 … … 98 104 simul_starttrace(femmodel); 99 105 100 /*Set X as our new control input a bd as INDEPENDENT!!*/106 /*Set X as our new control input and as INDEPENDENT*/ 101 107 #ifdef _HAVE_AD_ 102 108 IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices,"t"); … … 112 118 SetControlInputsFromVectorx(femmodel,aX); 113 119 xDelete<IssmDouble>(aX); 114 120 115 121 /*Compute solution (forward)*/ 116 122 void (*solutioncore)(FemModel*)=NULL; 117 123 CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); 118 124 solutioncore(femmodel); 125 126 127 /*Reset the time to zero for next optimization*/ 128 if(solution_type==TransientSolutionEnum){ 129 IssmDouble restart_time; 130 131 femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum); 132 femmodel->parameters->SetParam(restart_time,TimeEnum); 133 134 } 119 135 120 136 /*Get Dependents*/ … … 123 139 IssmPDouble *dependents; 124 140 DataSet* dependent_objects=NULL; 125 141 IssmDouble J=0.; 142 126 143 femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); 127 144 femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); … … 131 148 for(int i=0;i<dependent_objects->Size();i++){ 132 149 DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i); 133 dep->Responsex(&output_value,femmodel); 150 if(solution_type==TransientSolutionEnum) output_value = dep->GetValue(); 151 if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel); 134 152 if (my_rank==0) { 135 153 output_value>>=dependents[i]; 154 J+=output_value; 136 155 } 137 156 } … … 218 237 /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so 219 238 * as to generate num_dependents gradients: */ 220 totalgradient=xNewZeroInit<IssmPDouble>(num_independents );239 totalgradient=xNewZeroInit<IssmPDouble>(num_independents_old); 221 240 222 241 for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){ … … 253 272 /*Broadcast gradient to other ranks*/ 254 273 ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 255 256 274 /*Check size of Jlist to avoid crashes*/ 257 275 _assert_((*Jlisti)<JlistM); … … 259 277 260 278 /*Compute objective function*/ 261 IssmDouble* Jtemp = NULL;262 IssmDouble J;263 femmodel->CostFunctionx(&J,&Jtemp,NULL);264 279 *pf = reCast<double>(J); 265 280 _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | "); 266 281 267 282 /*Record cost function values and delete Jtemp*/ 268 for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = reCast<IssmPDouble>( Jtemp[i]);283 for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = reCast<IssmPDouble>(dependents[i]); 269 284 Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J); 270 xDelete<IssmDouble>(Jtemp);271 285 272 286 if(*indic==0){ … … 285 299 286 300 /*Compute gradient*/ 287 for(long i=0;i<num_independents_old;i++) 301 for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i]; 288 302 289 303 /*Constrain Gradient*/ … … 359 373 Vector<double> *Xpetsc = NULL; 360 374 361 /*THIS IS WHERE IT FAILS*/362 375 GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); 363 376 X = Xpetsc->ToMPISerial(); … … 397 410 /*Initialize Gradient and cost function of M1QN3*/ 398 411 indic = 4; /*gradient required*/ 399 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__);412 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); 400 413 simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct); 414 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); 401 415 /*Estimation of the expected decrease in f during the first iteration*/ 402 416 double df1=f; 403 417 418 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); 404 419 /*Call M1QN3 solver*/ 405 420 m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_, … … 407 422 >tol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz, 408 423 &reverse,&indic,izs,rzs,(void*)&mystruct); 424 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); 409 425 switch(int(omode)){ 410 426 case 0: _printf0_(" Stop requested (indic = 0)\n"); break; … … 420 436 /*Constrain solution vector*/ 421 437 double *XL = NULL; 438 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); 422 439 double *XU = NULL; 423 440 GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); 441 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); 424 442 GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); 425 443 444 printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); 426 445 for(int i=0;i<numberofvertices;i++){ 427 446 for(int c=0;c<num_controls;c++){ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
r22607 r22612 86 86 87 87 parameters->AddObject(new IntVecParam(InversionControlParametersEnum,ind_enums,num_independent_objects)); 88 //iomodel->FetchData(&num_costfunc,"md.numberedcostfunction.num_cost_functions");88 iomodel->FindConstant(&cm_responses,&num_costfunc,"md.autodiff.dependent_object_names"); 89 89 _assert_(num_costfunc>0); 90 iomodel->FindConstant(&cm_responses,&num_costfunc,"md.autodiff.dependent_object_names");91 90 if(num_costfunc<1) _error_ ("no cost functions found"); 92 91 int* costfunc_enums=xNew<int>(num_costfunc); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp
r22515 r22612 249 249 /*}}}*/ 250 250 } 251 else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){ 252 /*Deal with cfdragcoeffabsgrad: {{{*/ 253 254 /*cfdragcoeffabsgrad variables: */ 255 int num_cfdragcoeffabsgrads; 256 char** cfdragcoeffabsgrad_name_s = NULL; 257 char** cfdragcoeffabsgrad_definitionstring_s = NULL; 258 IssmDouble** cfdragcoeffabsgrad_weights_s = NULL; 259 int* cfdragcoeffabsgrad_weights_M_s = NULL; 260 int* cfdragcoeffabsgrad_weights_N_s = NULL; 261 char** cfdragcoeffabsgrad_weights_string_s = NULL; 262 int* cfdragcoeffabsgrad_datatime_s = NULL; 263 264 /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfdragcoeffabsgrad.m): */ 265 iomodel->FetchMultipleData(&cfdragcoeffabsgrad_name_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.name"); 266 iomodel->FetchMultipleData(&cfdragcoeffabsgrad_definitionstring_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.definitionstring"); 267 iomodel->FetchMultipleData(&cfdragcoeffabsgrad_weights_s,&cfdragcoeffabsgrad_weights_M_s,&cfdragcoeffabsgrad_weights_N_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.weights"); 268 iomodel->FetchMultipleData(&cfdragcoeffabsgrad_weights_string_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.weights_string"); 269 iomodel->FetchMultipleData(&cfdragcoeffabsgrad_datatime_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.datatime"); 270 271 for(j=0;j<num_cfdragcoeffabsgrads;j++){ 272 273 int weight_vector_type=0; 274 if ((cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofvertices) || (cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofvertices+1)){ 275 weight_vector_type=1; 276 } 277 else if ((cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofelements) || (cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofelements+1)){ 278 weight_vector_type=2; 279 } 280 else 281 _error_("cfdragcoeffabsgrad weight size not supported yet"); 282 283 /*First create a cfdragcoeffabsgrad object for that specific string (cfdragcoeffabsgrad_model_string_s[j]):*/ 284 output_definitions->AddObject(new Cfdragcoeffabsgrad(cfdragcoeffabsgrad_name_s[j],StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),cfdragcoeffabsgrad_datatime_s[j],false)); 285 286 /*Now, for this particular cfdragcoeffabsgrad object, make sure we plug into the elements: the observation, and the weights.*/ 287 for(int k=0;k<elements->Size();k++){ 288 289 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k)); 290 291 element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j], iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),7,WeightsSurfaceObservationEnum); 292 293 } 294 295 } 296 297 /*Free ressources:*/ 298 for(j=0;j<num_cfdragcoeffabsgrads;j++){ 299 char* string=NULL; 300 IssmDouble* matrix = NULL; 301 302 string = cfdragcoeffabsgrad_definitionstring_s[j]; xDelete<char>(string); 303 string = cfdragcoeffabsgrad_weights_string_s[j]; xDelete<char>(string); 304 string = cfdragcoeffabsgrad_name_s[j]; xDelete<char>(string); 305 matrix = cfdragcoeffabsgrad_weights_s[j]; xDelete<IssmDouble>(matrix); 306 } 307 xDelete<char*>(cfdragcoeffabsgrad_name_s); 308 xDelete<char*>(cfdragcoeffabsgrad_definitionstring_s); 309 xDelete<IssmDouble*>(cfdragcoeffabsgrad_weights_s); 310 xDelete<int>(cfdragcoeffabsgrad_weights_M_s); 311 xDelete<int>(cfdragcoeffabsgrad_weights_N_s); 312 xDelete<char*>(cfdragcoeffabsgrad_weights_string_s); 313 xDelete<int>(cfdragcoeffabsgrad_datatime_s); 314 /*}}}*/ 315 } 251 316 else if (output_definition_enums[i]==CfsurfacelogvelEnum){ 252 317 /*Deal with cfsurfacelogvel: {{{*/ … … 552 617 for(int k=0;k<elements->Size();k++){ 553 618 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k)); 554 element->DatasetInputCreate(cost_functions_weights[ i],cost_functions_weights_M[i],cost_functions_weights_N[i],cost_function_enums,num_cost_functions,iomodel,InversionCostFunctionsCoefficientsEnum);619 element->DatasetInputCreate(cost_functions_weights[j],cost_functions_weights_M[j],cost_functions_weights_N[j],cost_function_enums,num_cost_functions,iomodel,InversionCostFunctionsCoefficientsEnum); 555 620 } 556 621 output_definitions->AddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums)); -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r22611 r22612 636 636 CfsurfacelogvelEnum, 637 637 CfsurfacesquareEnum, 638 CfdragcoeffabsgradEnum, 638 639 ClosedEnum, 639 640 ColinearEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r22611 r22612 640 640 case CfsurfacelogvelEnum : return "Cfsurfacelogvel"; 641 641 case CfsurfacesquareEnum : return "Cfsurfacesquare"; 642 case CfdragcoeffabsgradEnum : return "Cfdragcoeffabsgrad"; 642 643 case ClosedEnum : return "Closed"; 643 644 case ColinearEnum : return "Colinear"; -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r22611 r22612 655 655 else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; 656 656 else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; 657 else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; 657 658 else if (strcmp(name,"Closed")==0) return ClosedEnum; 658 659 else if (strcmp(name,"Colinear")==0) return ColinearEnum; … … 751 752 else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; 752 753 else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum; 753 else if (strcmp(name,"HydrologydcEplColapseThickness")==0) return HydrologydcEplColapseThicknessEnum;754 754 else stage=7; 755 755 } 756 756 if(stage==7){ 757 if (strcmp(name,"HydrologydcEplCompressibility")==0) return HydrologydcEplCompressibilityEnum; 757 if (strcmp(name,"HydrologydcEplColapseThickness")==0) return HydrologydcEplColapseThicknessEnum; 758 else if (strcmp(name,"HydrologydcEplCompressibility")==0) return HydrologydcEplCompressibilityEnum; 758 759 else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum; 759 760 else if (strcmp(name,"HydrologydcEplInitialThickness")==0) return HydrologydcEplInitialThicknessEnum; … … 874 875 else if (strcmp(name,"Mumps")==0) return MumpsEnum; 875 876 else if (strcmp(name,"NewDamage")==0) return NewDamageEnum; 876 else if (strcmp(name,"Nodal")==0) return NodalEnum;877 877 else stage=8; 878 878 } 879 879 if(stage==8){ 880 if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; 880 if (strcmp(name,"Nodal")==0) return NodalEnum; 881 else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; 881 882 else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; 882 883 else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; … … 997 998 else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; 998 999 else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; 999 else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;1000 1000 else stage=9; 1001 1001 } 1002 1002 if(stage==9){ 1003 if (strcmp(name,"TransientParam")==0) return TransientParamEnum; 1003 if (strcmp(name,"TransientInput")==0) return TransientInputEnum; 1004 else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; 1004 1005 else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; 1005 1006 else if (strcmp(name,"Tria")==0) return TriaEnum; … … 1120 1121 else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; 1121 1122 else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; 1122 else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;1123 1123 else stage=10; 1124 1124 } 1125 1125 if(stage==10){ 1126 if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; 1126 if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; 1127 else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; 1127 1128 else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum; 1128 1129 else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; -
issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp
r22515 r22612 13 13 int input_enum = -1; 14 14 15 if(strcmp(string_in,"Thickness")==0 ){15 if(strcmp(string_in,"Thickness")==0 || strcmp(string_in,"md.geometry.thickness")==0){ 16 16 const char* field = "md.geometry.thickness"; 17 17 input_enum = ThicknessEnum;
Note:
See TracChangeset
for help on using the changeset viewer.