Index: /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 2267) +++ /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 2268) @@ -86,4 +86,10 @@ param= new Param(count,"eps_cm",DOUBLE); param->SetDouble(iomodel->eps_cm); + parameters->AddObject(param); + + /*cm_noisedampening: */ + count++; + param= new Param(count,"cm_noisedampening",DOUBLE); + param->SetDouble(iomodel->cm_noisedampening); parameters->AddObject(param); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 2267) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 2268) @@ -116,4 +116,5 @@ iomodel->tolx=0; iomodel->maxiter=NULL; + iomodel->cm_noisedampening=0; iomodel->mincontrolconstraint=0; iomodel->maxcontrolconstraint=0; @@ -331,4 +332,5 @@ IoModelFetchData((void**)&iomodel->eps_cm,NULL,NULL,iomodel_handle,"eps_cm","Scalar",NULL); IoModelFetchData((void**)&iomodel->tolx,NULL,NULL,iomodel_handle,"tolx","Scalar",NULL); + IoModelFetchData((void**)&iomodel->cm_noisedampening,NULL,NULL,iomodel_handle,"cm_noisedampening","Scalar",NULL); IoModelFetchData((void**)&iomodel->mincontrolconstraint,NULL,NULL,iomodel_handle,"mincontrolconstraint","Scalar",NULL); IoModelFetchData((void**)&iomodel->maxcontrolconstraint,NULL,NULL,iomodel_handle,"maxcontrolconstraint","Scalar",NULL); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 2267) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 2268) @@ -114,4 +114,5 @@ double tolx; double* maxiter; + double cm_noisedampening; double mincontrolconstraint; double maxcontrolconstraint; Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 2267) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 2268) @@ -2227,4 +2227,5 @@ int doflist1[numgrids]; int numberofdofspernode; + double dh1dh2dh3_basic[NDOF2][numgrids]; /* grid data: */ @@ -2250,11 +2251,9 @@ /* parameters: */ - double dvx[NDOF2]; - double dvy[NDOF2]; - double dadjx[NDOF2]; - double dadjy[NDOF2]; + double dk[NDOF2]; double vx,vy; double lambda,mu; double bed,thickness,Neff; + double cm_noisedampening; /*drag: */ @@ -2280,4 +2279,7 @@ /*recover pointers: */ inputs=(ParameterInputs*)vinputs; + + /*Get out if shelf*/ + if(shelf) return; /* Get node coordinates and dof list: */ @@ -2298,4 +2300,7 @@ if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ throw ErrorException(__FUNCT__,"missing adjoint input parameter"); + } + if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){ + throw ErrorException(__FUNCT__,"missing cm_noisedampening"); } @@ -2392,7 +2397,15 @@ GetNodalFunctions(l1l2l3, gauss_l1l2l3); + /*Get nodal functions derivatives*/ + GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get k derivative: dk/dx */ + GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ for (i=0;iRecover("adjoint",&adjxyz_list[0][0],3,dofs3,numgrids,(void**)nodes)){ throw ErrorException(__FUNCT__,"missing adjoint input parameter"); + } + if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){ + throw ErrorException(__FUNCT__,"missing cm_noisedampening"); } @@ -2599,6 +2618,13 @@ GetNodalFunctions(l1l2l3, gauss_l1l2l3); + /*Get nodal functions derivatives*/ + GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get k derivative: dk/dx */ + GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ for (i=0;iRecover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ throw ErrorException(__FUNCT__,"missing adjoint input parameter"); + } + if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){ + throw ErrorException(__FUNCT__,"missing cm_noisedampening"); } @@ -2781,7 +2816,17 @@ #endif + /*Get nodal functions derivatives*/ + GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get k derivative: dk/dx */ + GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); + /*Build gradje_g_gaussian vector (actually -dJ/dB): */ for (i=0;iRecover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ throw ErrorException(__FUNCT__,"missing velocity input parameter"); + } + if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){ + throw ErrorException(__FUNCT__,"missing cm_noisedampening"); } @@ -2933,4 +2985,25 @@ #endif + + + /*Initialize Jelem with dampening term*/ + /* NOT working for NOW + if (strcmp(control_type,"drag")==0 & !shelf){ + GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); + Jelem+=cm_noisedampening*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; + if (id==1){ + printf("id=%i value=%g k=[%g %g %g]\n",id,(pow(dk[0],2)+pow(dk[1],2)),k[0],k[1],k[2]); + } + if ((pow(dk[0],2)+pow(dk[1],2))>pow(10,-20)){ + printf("id=%i value=%g k=[%g %g %g]\n",id,(pow(dk[0],2)+pow(dk[1],2)),k[0],k[1],k[2]); + } + } + else if (strcmp(control_type,"B")==0){ + B=matice->GetB(); + GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); + Jelem+=cm_noisedampening*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; + } + */ + /*Differents misfits are allowed: */ if(fit==0){ Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 2267) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 2268) @@ -33,9 +33,9 @@ DataSet* results=NULL; DataSet* processed_results=NULL; - Result* result=NULL; + Result* result=NULL; ParameterInputs* inputs=NULL; int waitonlock=0; - int numberofnodes; + int numberofnodes; double* u_g_initial=NULL; @@ -44,4 +44,5 @@ int count; DataSet* parameters=NULL; + double cm_noisedampening; MODULEBOOT(); @@ -99,4 +100,6 @@ inputs->Add("velocity",u_g_initial,3,numberofnodes); if(control_analysis){ + model->FindParam(&cm_noisedampening,"cm_noisedampening"); + inputs->Add("cm_noisedampening",cm_noisedampening); model->FindParam(&u_g_obs,"u_g_obs",DiagnosticAnalysisEnum(),HorizAnalysisEnum()); inputs->Add("velocity_obs",u_g_obs,2,numberofnodes); Index: /issm/trunk/src/c/parallel/steadystate.cpp =================================================================== --- /issm/trunk/src/c/parallel/steadystate.cpp (revision 2267) +++ /issm/trunk/src/c/parallel/steadystate.cpp (revision 2268) @@ -45,4 +45,5 @@ double* u_g_obs=NULL; double dt; + double cm_noisedampening; Param* param=NULL; @@ -109,4 +110,6 @@ if(control_analysis){ + model->FindParam(&cm_noisedampening,"cm_noisedampening"); + inputs->Add("cm_noisedampening",cm_noisedampening); model->FindParam(&u_g_obs,"u_g_obs",DiagnosticAnalysisEnum(),HorizAnalysisEnum()); inputs->Add("velocity_obs",u_g_obs,2,numberofnodes); Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 2267) +++ /issm/trunk/src/m/classes/@model/model.m (revision 2268) @@ -182,4 +182,5 @@ md.tolx=0; md.optscal=[]; + md.cm_noisedampening=0; md.mincontrolconstraint=0; md.maxcontrolconstraint=0; Index: /issm/trunk/src/m/classes/public/display/displaycontrol.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displaycontrol.m (revision 2267) +++ /issm/trunk/src/m/classes/public/display/displaycontrol.m (revision 2268) @@ -19,4 +19,5 @@ fielddisplay(md,'maxiter','maximum iterations during each optimization step'); fielddisplay(md,'tolx','minimum tolerance which will stop one optimization search'); + fielddisplay(md,'cm_noisedampening','noise dampening coefficient'); fielddisplay(md,'mincontrolconstraint','minimum contraint for the controlled parameters'); fielddisplay(md,'maxcontrolconstraint','maximum contraint for the controlled parameters'); Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 2267) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 2268) @@ -112,4 +112,5 @@ WriteData(fid,md.tolx,'Scalar','tolx'); WriteData(fid,md.maxiter,'Mat','maxiter'); +WriteData(fid,md.cm_noisedampening,'Scalar','cm_noisedampening'); WriteData(fid,md.mincontrolconstraint,'Scalar','mincontrolconstraint'); WriteData(fid,md.maxcontrolconstraint,'Scalar','maxcontrolconstraint'); Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic.m (revision 2267) +++ /issm/trunk/src/m/solutions/cielo/diagnostic.m (revision 2268) @@ -34,4 +34,5 @@ inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes); if md.control_analysis, + inputs=add(inputs,'cm_noisedampening',models.dh.parameters.cm_noisedampening,'double'); inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes); end Index: /issm/trunk/src/m/solutions/cielo/steadystate.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/steadystate.m (revision 2267) +++ /issm/trunk/src/m/solutions/cielo/steadystate.m (revision 2268) @@ -44,4 +44,5 @@ inputs=add(inputs,'dt',models.t.parameters.dt*models.t.parameters.yts,'double'); if md.control_analysis, + inputs=add(inputs,'cm_noisedampening',models.dh.parameters.cm_noisedampening,'double'); inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes); end