Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 1670) +++ /issm/trunk/src/c/Makefile.am (revision 1671) @@ -587,5 +587,5 @@ ./parallel/control.cpp\ ./parallel/control_core.cpp\ - ./parallel/ControlPrepareStokes.cpp\ + ./parallel/ControlInitialization.cpp\ ./parallel/ControlTemporaryResults.cpp\ ./parallel/objectivefunctionC.cpp\ Index: /issm/trunk/src/c/parallel/ControlInitialization.cpp =================================================================== --- /issm/trunk/src/c/parallel/ControlInitialization.cpp (revision 1671) +++ /issm/trunk/src/c/parallel/ControlInitialization.cpp (revision 1671) @@ -0,0 +1,145 @@ +/*!\file: ControlInitialization.cpp + * \brief: ... + */ + +#undef __FUNCT__ +#define __FUNCT__ "ControlInitialization" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void ControlInitialization(FemModel** pmodel, 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"); + + /*if no Stokes, assign output and return*/ + if (!isstokes){ + *pmodel=fem_dh; + return; + } + + /*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); + + /*Assign output*/ + *pmodel=fem_ds; +} Index: sm/trunk/src/c/parallel/ControlPrepareStokes.cpp =================================================================== --- /issm/trunk/src/c/parallel/ControlPrepareStokes.cpp (revision 1670) +++ (revision ) @@ -1,136 +1,0 @@ -/*!\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/control_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/control_core.cpp (revision 1670) +++ /issm/trunk/src/c/parallel/control_core.cpp (revision 1671) @@ -13,10 +13,6 @@ 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; + /*fem model: */ + FemModel* fem_model=NULL; /*output: */ @@ -47,45 +43,26 @@ /*flags: */ - int analysis_type; - int sub_analysis_type; + int analysis_type; + int sub_analysis_type; 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; - /*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_dh->parameters->FindParam((void*)&ishutter,"ishutter"); - fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn"); - fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); - fem_dh->parameters->FindParam((void*)&isstokes,"isstokes"); + /*Process models*/ + ControlInitialization(&fem_model,fems,inputs); /*Recover parameters used throughout the solution:*/ - fem_dh->parameters->FindParam((void*)&nsteps,"nsteps"); - fem_dh->parameters->FindParam((void*)&control_type,"control_type"); - fem_dh->parameters->FindParam((void*)&fit,"fit"); - fem_dh->parameters->FindParam((void*)&optscal,"optscal"); - fem_dh->parameters->FindParam((void*)&maxiter,"maxiter"); - fem_dh->parameters->FindParam((void*)&tolx,"tolx"); - fem_dh->parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint"); - fem_dh->parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint"); - fem_dh->parameters->FindParam((void*)¶m_g,"param_g"); - fem_dh->parameters->FindParam((void*)&analysis_type,"analysis_type"); - fem_dh->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); - fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); - gsize=fem_dh->nodes->NumberOfDofs(); + fem_model->parameters->FindParam((void*)&nsteps,"nsteps"); + fem_model->parameters->FindParam((void*)&control_type,"control_type"); + fem_model->parameters->FindParam((void*)&fit,"fit"); + fem_model->parameters->FindParam((void*)&optscal,"optscal"); + fem_model->parameters->FindParam((void*)&maxiter,"maxiter"); + fem_model->parameters->FindParam((void*)&tolx,"tolx"); + fem_model->parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint"); + fem_model->parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint"); + fem_model->parameters->FindParam((void*)¶m_g,"param_g"); + fem_model->parameters->FindParam((void*)&analysis_type,"analysis_type"); + fem_model->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); + fem_model->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + gsize=fem_model->nodes->NumberOfDofs(); /*Initialize misfit: */ @@ -93,12 +70,9 @@ /*erase useless parameters: */ - param=(Param*)fem_dh->parameters->FindParamObject("param_g"); - fem_dh->parameters->DeleteObject((Object*)param); + param=(Param*)fem_model->parameters->FindParamObject("param_g"); + fem_model->parameters->DeleteObject((Object*)param); - param=(Param*)fem_dh->parameters->FindParamObject("u_g"); - fem_dh->parameters->DeleteObject((Object*)param); - - /*take care of Stokes: compute slope and get spc once for all*/ - if (isstokes) ControlPrepareStokes(fems,inputs); + param=(Param*)fem_model->parameters->FindParamObject("u_g"); + fem_model->parameters->DeleteObject((Object*)param); /*Start looping: */ @@ -110,10 +84,8 @@ /*Update parameters: */ - 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); + UpdateFromInputsx(fem_model->elements,fem_model->nodes,fem_model->loads, fem_model->materials,inputs); _printf_("%s\n"," computing gradJ..."); - if (isstokes) grad_g=GradJCompute(inputs,fem_ds); - else grad_g=GradJCompute(inputs,fem_dh); + grad_g=GradJCompute(inputs,fem_model); _printf_("%s\n"," done."); @@ -126,6 +98,5 @@ _printf_("%s\n"," optimizing along gradient direction..."); - if (isstokes) optargs.femmodel=fem_ds; - else optargs.femmodel=fem_dh; + optargs.femmodel=fem_model; 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]; Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 1670) +++ /issm/trunk/src/c/parallel/parallel.h (revision 1671) @@ -41,5 +41,5 @@ void WriteLockFile(char* filename); -void ControlPrepareStokes(FemModel* fems, ParameterInputs* inputs); +void ControlInitialization(FemModel** pmodel, FemModel* fems, ParameterInputs* inputs); void ControlTemporaryResults(FemModel* fems,double* param_g,double* J,int n,ParameterInputs* inputs);