Changeset 2268


Ignore:
Timestamp:
09/21/09 15:52:16 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added regularization term in control methods. Misfit needs to be completed

Location:
issm/trunk/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp

    r2267 r2268  
    8686                param= new Param(count,"eps_cm",DOUBLE);
    8787                param->SetDouble(iomodel->eps_cm);
     88                parameters->AddObject(param);
     89
     90                /*cm_noisedampening: */
     91                count++;
     92                param= new Param(count,"cm_noisedampening",DOUBLE);
     93                param->SetDouble(iomodel->cm_noisedampening);
    8894                parameters->AddObject(param);
    8995

  • issm/trunk/src/c/ModelProcessorx/IoModel.cpp

    r2267 r2268  
    116116        iomodel->tolx=0;
    117117        iomodel->maxiter=NULL;
     118        iomodel->cm_noisedampening=0;
    118119        iomodel->mincontrolconstraint=0;
    119120        iomodel->maxcontrolconstraint=0;
    … …  
    331332        IoModelFetchData((void**)&iomodel->eps_cm,NULL,NULL,iomodel_handle,"eps_cm","Scalar",NULL);
    332333        IoModelFetchData((void**)&iomodel->tolx,NULL,NULL,iomodel_handle,"tolx","Scalar",NULL);
     334        IoModelFetchData((void**)&iomodel->cm_noisedampening,NULL,NULL,iomodel_handle,"cm_noisedampening","Scalar",NULL);
    333335        IoModelFetchData((void**)&iomodel->mincontrolconstraint,NULL,NULL,iomodel_handle,"mincontrolconstraint","Scalar",NULL);
    334336        IoModelFetchData((void**)&iomodel->maxcontrolconstraint,NULL,NULL,iomodel_handle,"maxcontrolconstraint","Scalar",NULL);

  • issm/trunk/src/c/ModelProcessorx/IoModel.h

    r2267 r2268  
    114114        double  tolx;
    115115        double* maxiter;
     116        double  cm_noisedampening;
    116117        double  mincontrolconstraint;
    117118        double  maxcontrolconstraint;

  • issm/trunk/src/c/objects/Tria.cpp

    r2112 r2268  
    22272227        int          doflist1[numgrids];
    22282228        int          numberofdofspernode;
     2229        double       dh1dh2dh3_basic[NDOF2][numgrids];
    22292230
    22302231        /* grid data: */
    … …  
    22502251
    22512252        /* parameters: */
    2252         double  dvx[NDOF2];
    2253         double  dvy[NDOF2];
    2254         double  dadjx[NDOF2];
    2255         double  dadjy[NDOF2];
     2253        double  dk[NDOF2];
    22562254        double  vx,vy;
    22572255        double  lambda,mu;
    22582256        double  bed,thickness,Neff;
     2257        double  cm_noisedampening;
    22592258       
    22602259        /*drag: */
    … …  
    22802279        /*recover pointers: */
    22812280        inputs=(ParameterInputs*)vinputs;
     2281
     2282        /*Get out if shelf*/
     2283        if(shelf) return;
    22822284
    22832285        /* Get node coordinates and dof list: */
    … …  
    22982300        if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
    22992301                throw ErrorException(__FUNCT__,"missing adjoint input parameter");
     2302        }
     2303        if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){
     2304                throw ErrorException(__FUNCT__,"missing cm_noisedampening");
    23002305        }
    23012306
    … …  
    23922397                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    23932398
     2399                /*Get nodal functions derivatives*/
     2400                GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
     2401
     2402                /*Get k derivative: dk/dx */
     2403                GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
     2404
    23942405                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    23952406                for (i=0;i<numgrids;i++){
    2396                         grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
     2407                        grade_g_gaussian[i]=
     2408                          -2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i]                         //standard term dJ/dki
     2409                          +cm_noisedampening*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);  // regularization term d/dki(1/2*(dk/dx)^2)
    23972410                }
    23982411               
    … …  
    24242437        int          doflist1[numgrids];
    24252438        int          numberofdofspernode;
     2439        double       dh1dh2dh3_basic[NDOF2][numgrids];
    24262440
    24272441        /* grid data: */
    … …  
    24542468        double  surface_normal[3];
    24552469        double  bed_normal[3];
     2470        double  dk[NDOF2];
     2471        double  cm_noisedampening;
    24562472
    24572473        /*drag: */
    … …  
    24942510        if(!inputs->Recover("adjoint",&adjxyz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
    24952511                throw ErrorException(__FUNCT__,"missing adjoint input parameter");
     2512        }
     2513        if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){
     2514                throw ErrorException(__FUNCT__,"missing cm_noisedampening");
    24962515        }
    24972516
    … …  
    25992618                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    26002619
     2620                /*Get nodal functions derivatives*/
     2621                GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
     2622
     2623                /*Get k derivative: dk/dx */
     2624                GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
     2625
    26012626                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    26022627                for (i=0;i<numgrids;i++){
     2628                        //standard gradient dJ/dki
    26032629                        grade_g_gaussian[i]=(
    26042630                                                -lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
    … …  
    26062632                                                -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
    26072633                                                )*Jdet*gauss_weight*l1l2l3[i];
     2634
     2635                        //Add regularization term
     2636                        grade_g_gaussian[i]+=cm_noisedampening*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
    26082637                }
    26092638
    … …  
    26662695        int          doflist1[numgrids];
    26672696        int          numberofdofspernode;
     2697        double       dh1dh2dh3_basic[NDOF2][numgrids];
    26682698
    26692699        /* grid data: */
    … …  
    27102740        double  thickness;
    27112741        int     dofs[1]={0};
     2742        double  dk[NDOF2];
     2743        double  cm_noisedampening;
    27122744       
    27132745        ParameterInputs* inputs=NULL;
    … …  
    27312763        if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
    27322764                throw ErrorException(__FUNCT__,"missing adjoint input parameter");
     2765        }
     2766        if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){
     2767                throw ErrorException(__FUNCT__,"missing cm_noisedampening");
    27332768        }
    27342769
    … …  
    27812816                #endif
    27822817
     2818                /*Get nodal functions derivatives*/
     2819                GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
     2820
     2821                /*Get k derivative: dk/dx */
     2822                GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
     2823
    27832824                /*Build gradje_g_gaussian vector (actually -dJ/dB): */
    27842825                for (i=0;i<numgrids;i++){
     2826                        //standard gradient dJ/dki
    27852827                        grade_g_gaussian[i]=-viscosity_complement*thickness*( (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1])*Jdet*gauss_weight*l1l2l3[i];
     2828
     2829                        //Add regularization term
     2830                        grade_g_gaussian[i]+=cm_noisedampening*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
    27862831                }
    27872832
    … …  
    28292874        double relative_list[numgrids];
    28302875        double logarithmic_list[numgrids];
     2876        double B[numgrids];
    28312877
    28322878        /* gaussian points: */
    … …  
    28422888        double  velocity_mag,obs_velocity_mag;
    28432889        double  absolute,relative,logarithmic;
     2890        double  dk[NDOF2];
     2891        double  dB[NDOF2];
     2892        double  cm_noisedampening;
    28442893
    28452894        /* Jacobian: */
    … …  
    28702919        if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
    28712920                throw ErrorException(__FUNCT__,"missing velocity input parameter");
     2921        }
     2922        if(!inputs->Recover("cm_noisedampening",&cm_noisedampening)){
     2923                throw ErrorException(__FUNCT__,"missing cm_noisedampening");
    28722924        }
    28732925
    … …  
    29332985                #endif
    29342986               
     2987
     2988
     2989                /*Initialize Jelem with dampening term*/
     2990                /* NOT working for NOW
     2991                if (strcmp(control_type,"drag")==0 & !shelf){
     2992                        GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
     2993                        Jelem+=cm_noisedampening*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
     2994                        if (id==1){
     2995                                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]);
     2996                        }
     2997                        if ((pow(dk[0],2)+pow(dk[1],2))>pow(10,-20)){
     2998                                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]);
     2999                        }
     3000                }
     3001                else if (strcmp(control_type,"B")==0){
     3002                        B=matice->GetB();
     3003                        GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3);
     3004                        Jelem+=cm_noisedampening*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
     3005                }
     3006                */
     3007
    29353008                /*Differents misfits are allowed: */
    29363009                if(fit==0){

  • issm/trunk/src/c/parallel/diagnostic.cpp

    r2234 r2268  
    3333        DataSet* results=NULL;
    3434        DataSet* processed_results=NULL;
    35         Result* result=NULL;
     35        Result*  result=NULL;
    3636       
    3737        ParameterInputs* inputs=NULL;
    3838        int waitonlock=0;
    39         int   numberofnodes;
     39        int numberofnodes;
    4040       
    4141        double* u_g_initial=NULL;
    … …  
    4444        int      count;
    4545        DataSet* parameters=NULL;
     46        double cm_noisedampening;
    4647
    4748        MODULEBOOT();
    … …  
    99100        inputs->Add("velocity",u_g_initial,3,numberofnodes);
    100101        if(control_analysis){
     102                model->FindParam(&cm_noisedampening,"cm_noisedampening");
     103                inputs->Add("cm_noisedampening",cm_noisedampening);
    101104                model->FindParam(&u_g_obs,"u_g_obs",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    102105                inputs->Add("velocity_obs",u_g_obs,2,numberofnodes);

  • issm/trunk/src/c/parallel/steadystate.cpp

    r2234 r2268  
    4545        double* u_g_obs=NULL;
    4646        double  dt;
     47        double  cm_noisedampening;
    4748        Param*  param=NULL;
    4849
    … …  
    109110
    110111        if(control_analysis){
     112                model->FindParam(&cm_noisedampening,"cm_noisedampening");
     113                inputs->Add("cm_noisedampening",cm_noisedampening);
    111114                model->FindParam(&u_g_obs,"u_g_obs",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    112115                inputs->Add("velocity_obs",u_g_obs,2,numberofnodes);

  • issm/trunk/src/m/classes/@model/model.m

    r2267 r2268  
    182182        md.tolx=0;
    183183        md.optscal=[];
     184        md.cm_noisedampening=0;
    184185        md.mincontrolconstraint=0;
    185186        md.maxcontrolconstraint=0;

  • issm/trunk/src/m/classes/public/display/displaycontrol.m

    r2267 r2268  
    1919        fielddisplay(md,'maxiter','maximum iterations during each optimization step');
    2020        fielddisplay(md,'tolx','minimum tolerance which will stop one optimization search');
     21        fielddisplay(md,'cm_noisedampening','noise dampening coefficient');
    2122        fielddisplay(md,'mincontrolconstraint','minimum contraint for the controlled parameters');
    2223        fielddisplay(md,'maxcontrolconstraint','maximum contraint for the controlled parameters');

  • issm/trunk/src/m/classes/public/marshall.m

    r2267 r2268  
    112112WriteData(fid,md.tolx,'Scalar','tolx');
    113113WriteData(fid,md.maxiter,'Mat','maxiter');
     114WriteData(fid,md.cm_noisedampening,'Scalar','cm_noisedampening');
    114115WriteData(fid,md.mincontrolconstraint,'Scalar','mincontrolconstraint');
    115116WriteData(fid,md.maxcontrolconstraint,'Scalar','maxcontrolconstraint');

  • issm/trunk/src/m/solutions/cielo/diagnostic.m

    r1963 r2268  
    3434        inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
    3535        if md.control_analysis,
     36                inputs=add(inputs,'cm_noisedampening',models.dh.parameters.cm_noisedampening,'double');
    3637                inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes);
    3738        end

  • issm/trunk/src/m/solutions/cielo/steadystate.m

    r2027 r2268  
    4444        inputs=add(inputs,'dt',models.t.parameters.dt*models.t.parameters.yts,'double');
    4545        if md.control_analysis,
     46                inputs=add(inputs,'cm_noisedampening',models.dh.parameters.cm_noisedampening,'double');
    4647                inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes);
    4748        end
Note: See TracChangeset for help on using the changeset viewer.