Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 1208) +++ /issm/trunk/src/c/Makefile.am (revision 1209) @@ -574,4 +574,5 @@ ./parallel/control.cpp\ ./parallel/control_core.cpp\ + ./parallel/ControlPrepareStokes.cpp\ ./parallel/objectivefunctionC.cpp\ ./parallel/GradJCompute.cpp\ Index: /issm/trunk/src/c/parallel/ControlPrepareStokes.cpp =================================================================== --- /issm/trunk/src/c/parallel/ControlPrepareStokes.cpp (revision 1209) +++ /issm/trunk/src/c/parallel/ControlPrepareStokes.cpp (revision 1209) @@ -0,0 +1,136 @@ +/*!\file: ControlPrepareStokes.cpp + * \brief: ... + */ + +#undef __FUNCT__ +#define __FUNCT__ "ControlPrepareStokes" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void ControlPrepareStokes(FemModel* fems, ParameterInputs* inputs){ + + extern int my_rank; + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + + /*solutions: */ + Vec ug=NULL; + Vec ug_horiz=NULL; + Vec ug_vert=NULL; + Vec ug_stokes=NULL; + Vec pg=NULL; + Vec slopex=NULL; + Vec slopey=NULL; + + /*flags: */ + int debug=0; + int dim=-1; + int ishutter=0; + int ismacayealpattyn=0; + int isstokes=0; + int numberofdofspernode_sl; + int numberofdofspernode_dh; + int numberofdofspernode_ds; + int numberofnodes; + + double stokesreconditioning; + + /*dof recovery: */ + int dof01[2]={0,1}; + int dof2[1]={2}; + int dof012[3]={0,1,2}; + int dof3[1]={3}; + + /*recover fem models: */ + fem_dh=fems+0; + fem_dv=fems+1; + fem_ds=fems+2; + fem_dhu=fems+3; + fem_sl=fems+4; + + //first recover parameters common to all solutions + fem_dh->parameters->FindParam((void*)&debug,"debug"); + fem_dh->parameters->FindParam((void*)&dim,"dim"); + fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter"); + fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn"); + fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + fem_ds->parameters->FindParam((void*)&isstokes,"isstokes"); + fem_ds->parameters->FindParam((void*)&stokesreconditioning,"stokesreconditioning"); + + //specific parameters for specific models + fem_dh->parameters->FindParam((void*)&numberofdofspernode_dh,"numberofdofspernode"); + fem_sl->parameters->FindParam((void*)&numberofdofspernode_sl,"numberofdofspernode"); + fem_ds->parameters->FindParam((void*)&numberofdofspernode_ds,"numberofdofspernode"); + + /*1: compute slopes once for all*/ + + //compute slopes + if(debug)_printf_("%s\n","computing bed slope (x and y derivatives)..."); + diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum()); + diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum()); + FieldExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopex",0); + FieldExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopey",0); + + //Add in inputs + inputs->Add("bedslopex",slopex,numberofdofspernode_sl,numberofnodes); + inputs->Add("bedslopey",slopey,numberofdofspernode_sl,numberofnodes); + + /*2: run a complete diagnostic to update spcs*/ + + //horizontal velocity + if(debug)_printf_("%s\n"," computing horizontal velocities..."); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + if(debug)_printf_("%s\n"," extruding horizontal velocities..."); + VecDuplicatePatch(&ug_horiz,ug); FieldExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials,"velocity",1); + + //pressure + if(debug)_printf_("%s\n"," computing pressure according to MacAyeal..."); + ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->loads, fem_dh->materials, numberofnodes); + + //vertical velocity + if(debug)_printf_("%s\n"," computing vertical velocities..."); + inputs->Add("velocity",ug_horiz,numberofdofspernode_dh,numberofnodes); + diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticAnalysisEnum(),VertAnalysisEnum()); + + //Create 3d u_g + if(debug)_printf_("%s\n"," combining horizontal and vertical velocities..."); + VecFree(&ug); ug=NewVec(numberofnodes*3); + VecMerge(ug,ug_horiz,dofsetgen(2,&dof01[0],3,numberofnodes*3),numberofnodes*2); + VecMerge(ug,ug_vert,dofsetgen(1,&dof2[0],3,numberofnodes*3),numberofnodes*1); + + //Create 4d u_g + if(debug)_printf_("%s\n"," computing pressure according to Pattyn..."); + ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->loads, fem_dh->materials, numberofnodes); + VecScale(pg,1.0/stokesreconditioning); + ug_stokes=NewVec(fem_ds->nodesets->GetGSize()); + VecMerge(ug_stokes,ug,dofsetgen(3,dof012,4,numberofnodes*4),numberofnodes*3); + VecMerge(ug_stokes,pg,dofsetgen(1,dof3,4,numberofnodes*4),numberofnodes); + + //Add in inputs + inputs->Add("velocity",ug_stokes,numberofdofspernode_ds,numberofnodes); + + //update spcs + if(debug)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); + VecMerge(fem_ds->yg,ug,dofsetgen(3,dof012,4,numberofnodes*4),3*numberofnodes); + VecFree(&fem_ds->ys); VecFree(&fem_ds->ys0); + Reducevectorgtosx(&fem_ds->ys,&fem_ds->ys0, fem_ds->yg,fem_ds->nodesets); + + //Compute Stokes velocities to speed up later runs + if(debug)_printf_("%s\n"," computing stokes velocities and pressure ..."); + VecFree(&ug); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); + VecFree(&pg); + + //Add in inputs + inputs->Add("velocity",ug,numberofdofspernode_ds,numberofnodes); +} Index: /issm/trunk/src/c/parallel/GradJCompute.cpp =================================================================== --- /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 1208) +++ /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 1209) @@ -57,5 +57,5 @@ _printf_("%s\n"," recover solution for this stiffness and right hand side:"); - diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type,sub_analysis_type); + diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes); Index: /issm/trunk/src/c/parallel/ProcessResults.cpp =================================================================== --- /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 1208) +++ /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 1209) @@ -101,5 +101,4 @@ /*Recover femmodels first: */ if(analysis_type==DiagnosticAnalysisEnum()){ - fem_dh=fems+0; fem_dv=fems+1; @@ -116,5 +115,4 @@ if(analysis_type==TransientAnalysisEnum()){ - fem_dh=fems+0; fem_dv=fems+1; @@ -144,5 +142,9 @@ if(analysis_type==ControlAnalysisEnum()){ - fem_dh=fems+0; //load u_g + fem_dh=fems+0; + fem_dv=fems+1; + fem_ds=fems+2; + fem_dhu=fems+3; + fem_sl=fems+4; fem_c=fems+0; //load param_g Index: /issm/trunk/src/c/parallel/control_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/control_core.cpp (revision 1208) +++ /issm/trunk/src/c/parallel/control_core.cpp (revision 1209) @@ -15,4 +15,8 @@ /*fem models: */ FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; /*output: */ @@ -57,4 +61,8 @@ /*recover fem models: */ fem_dh=fems+0; + fem_dv=fems+1; + fem_ds=fems+2; + fem_dhu=fems+3; + fem_sl=fems+4; //first recover parameters common to all solutions @@ -91,4 +99,9 @@ fem_dh->parameters->DeleteObject((Object*)param); + /*take care of Stokes: compute slope and get spc once for all*/ + if (isstokes){ + ControlPrepareStokes(fems,inputs); + } + /*Start looping: */ for(n=0;nelements,fem_dh->nodes,fem_dh->loads, fem_dh->materials,inputs); + if (isstokes){ + UpdateFromInputsx(fem_ds->elements,fem_ds->nodes,fem_ds->loads, fem_ds->materials,inputs); + } + else{ + UpdateFromInputsx(fem_dh->elements,fem_dh->nodes,fem_dh->loads, fem_dh->materials,inputs); + } _printf_("%s\n"," computing gradJ..."); - grad_g=GradJCompute(inputs,fem_dh); + if (isstokes){ + grad_g=GradJCompute(inputs,fem_ds); + } + else{ + grad_g=GradJCompute(inputs,fem_dh); + } _printf_("%s\n"," done."); @@ -113,5 +136,11 @@ _printf_("%s\n"," optimizing along gradient direction..."); - optargs.femmodel=fem_dh; optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n; + if (isstokes){ + optargs.femmodel=fem_ds; + } + else{ + optargs.femmodel=fem_dh; + } + optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n; optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n]; BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs); Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 1208) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 1209) @@ -77,5 +77,5 @@ //Run diagnostic with updated parameters. - diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type,sub_analysis_type); + diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 1208) +++ /issm/trunk/src/c/parallel/parallel.h (revision 1209) @@ -40,4 +40,6 @@ void WriteLockFile(char* filename); +void ControlPrepareStokes(FemModel* fems, ParameterInputs* inputs); + void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type); //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename);