Index: /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp =================================================================== --- /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp (revision 3906) +++ /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp (revision 3906) @@ -0,0 +1,79 @@ +/*!\file NodeConnectivityx + * \brief: compute node connectivity table, using elements connectivity table. + * + * For each node, we want to know how many elements are connected to this element, and which they are. + * Given that the 2d meshes we create in ISSM are triangular for now, and they are delaunay conforming, + * each triangle has a minimum angle of 30 degrees, which implies a connectivity <=6. We therefore return + * a nods x 7 connectivity table, with the 7'th column giving us the number of elements connected to each + * row node, and the first 6 columns giving us the elements numbers. + * Amend that: sounds like some triangles get up to 9 connectivity. Take 10 to be on the safe side. + * In order to be compatible with matlab output, the connectivity table is given in matlab indexing (starts at 1). + */ + +#include "./NodeConnectivityx.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void NodeConnectivityx( double** pconnectivity, int* pwidth, double* elements, int nel, int nods){ + + int i,j,n; + const int maxels=25; + const int width=maxels+1; + + /*intermediary: */ + int node; + int index; + int num_elements; + int already_plugged=0; + double element; + + /*output: */ + double* connectivity=NULL; + + + + /*Allocate connectivity: */ + connectivity=(double*)xcalloc(nods*width,sizeof(double)); + + /*Go through all elements, and for each elements, plug into the connectivity, all the nodes. + * If nodes are already plugged into the connectivity, skip them.: */ + for(n=0;nmaxels)ISSMERROR("%s%g%s"," max connectivity width reached (",*(connectivity+width*i+maxels),")! increase width of connectivity table"); + } + + /*Assign output pointers: */ + *pconnectivity=connectivity; + *pwidth=width; +} Index: /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h =================================================================== --- /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h (revision 3906) +++ /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h (revision 3906) @@ -0,0 +1,12 @@ +/*!\file: NodeConnectivityx.h + * \brief header file for node connectivity computation + */ + +#ifndef _NODECONNECTIVITYX_H +#define _NODECONNECTIVITYX_H + +/* local prototypes: */ +void NodeConnectivityx( double** pconnectivity, int* pwidth,double* elements, int nel, int nods); + +#endif /* _NODECONNECTIVITYX_H */ + Index: /issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.cpp =================================================================== --- /issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.cpp (revision 3906) +++ /issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.cpp (revision 3906) @@ -0,0 +1,66 @@ +/*!\file NormalizeConstraintsx + * \brief normalize and reduce constraints matrix Rmg to Gmn: + */ + +#include "./NormalizeConstraintsx.h" + +#include "../shared/shared.h" + +void NormalizeConstraintsx( Mat* pGmn, Mat Rmg, NodeSets* nodesets){ + + /*output: */ + Mat Gmn=NULL; + + /*intermediary: */ + int i; + Mat Rmm=NULL; + Mat Rmn=NULL; + Mat InvRmm=NULL; + int msize,nsize; + double* pv_m=NULL; + double* row_m=NULL; + double* pv_n=NULL; + int null=0; + + if(nodesets){ + /*Recover data: */ + msize=nodesets->GetMSize(); + nsize=nodesets->GetNSize(); + pv_m=nodesets->GetPV_M(); + pv_n=nodesets->GetPV_N(); + + + if (msize){ + + /*Build row partitioning vector for Rmm and Rmn: */ + row_m=(double*)xmalloc(msize*sizeof(double)); + for(i=0;iOutputRifts(riftproperties); + + /*Assemble vector: */ + VecAssemblyBegin(riftproperties); + VecAssemblyEnd(riftproperties); + + /*Assign output pointers: */ + *priftproperties=riftproperties; +} Index: /issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.h =================================================================== --- /issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.h (revision 3906) +++ /issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.h (revision 3906) @@ -0,0 +1,14 @@ +/*!\file: OutputRiftsx.h + * \brief header file for rift results output. + */ + +#ifndef _OUTPUTRIFTSX_H +#define _OUTPUTRIFTSX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void OutputRiftsx( Vec* priftproperties, DataSet* loads, int numrifts); + +#endif /* _OUTPUTRIFTSX_H */ + Index: /issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp =================================================================== --- /issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 3906) +++ /issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 3906) @@ -0,0 +1,50 @@ +/*!\file PenaltyConstraintsx + * \brief: set up penalty constraints on loads + */ + +#include "./PenaltyConstraintsx.h" +#include "./RiftConstraints.h" +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes, DataSet* vertices, + DataSet* loads,DataSet* materials, Parameters* parameters,int analysis_type,int sub_analysis_type){ + + int i; + + extern int num_procs; + extern int my_rank; + + /*output: */ + int converged=0; + int num_unstable_constraints=0; + int min_mechanical_constraints=0; + + /*recover parameters: */ + parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + + /*First, get nodes and loads configured: */ + elements->Configure(elements, loads, nodes,vertices, materials,parameters); + nodes->Configure(elements, loads, nodes,vertices, materials,parameters); + loads->Configure(elements, loads, nodes,vertices, materials,parameters); + + /*Do we have penalties linked to rifts? In this case, run our special rifts penalty + * management routine, otherwise, skip : */ + if (RiftIsPresent(loads)){ + RiftConstraints(&converged,&num_unstable_constraints,loads,min_mechanical_constraints,analysis_type,sub_analysis_type); + } + else if(loads->MeltingIsPresent()){ + loads->MeltingConstraints(&converged,&num_unstable_constraints,analysis_type,sub_analysis_type); + } + else{ + /*Do nothing, no constraints management!:*/ + num_unstable_constraints=0; + converged=1; + } + + /*Assign output pointers: */ + *pconverged=converged; + *pnum_unstable_constraints=num_unstable_constraints; +} Index: /issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.h =================================================================== --- /issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.h (revision 3906) +++ /issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.h (revision 3906) @@ -0,0 +1,16 @@ +/*!\file: PenaltyConstraintsx.h + * \brief header file for penalty constraints module + */ + +#ifndef _PENALTYCONSTRAINTSX_H +#define _PENALTYCONSTRAINTSX_H + +#include "../DataSet/DataSet.h" +#include "../objects/objects.h" + +/* local prototypes: */ +void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes, DataSet* vertices, + DataSet* loads,DataSet* materials, Parameters* parameters,int analysis_type,int sub_analysis_type); + +#endif /* _PENALTYCONSTRAINTSX_H */ + Index: /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp =================================================================== --- /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp (revision 3906) +++ /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp (revision 3906) @@ -0,0 +1,345 @@ +/*!\file RiftConstraints.cpp + * \brief: manage penalties for rifts + */ + +#include "./RiftConstraints.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../shared/shared.h" + +#define _ZIGZAGCOUNTER_ + +int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,int min_mechanical_constraints,int analysis_type,int sub_analysis_type){ + + int num_unstable_constraints=0; + int converged=0; + int potential; + extern int my_rank; + + Constrain(&num_unstable_constraints,loads,analysis_type); + if(num_unstable_constraints==0)converged=1; + + + if(IsFrozen(loads)){ + converged=1; + num_unstable_constraints=0; + } + else if(num_unstable_constraints<=min_mechanical_constraints){ + _printf_(" freezing constraints\n"); + FreezeConstraints(loads,analysis_type); + } + + /*Assign output pointers: */ + *pconverged=converged; + *pnum_unstable_constraints=num_unstable_constraints; +} + +int IsMaterialStable(DataSet* loads,int analysis_type){ + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and if non-linearity of the material has converged, let all penalties know: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + if (riftfront->IsMaterialStable(analysis_type)){ + found=1; + /*do not break! all penalties should get informed the non-linearity converged!*/ + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} + +int RiftIsPresent(DataSet* loads){ + + + int i; + + int found=0; + int mpi_found=0; + + /*go though loads, and figure out if one of the loads is a Riftfront: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum==loads->GetEnum(i)){ + found=1; + break; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} + +int IsPreStable(DataSet* loads){ + + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and figure out if one of the penpair loads is still not stable: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + if (riftfront->PreStable()==0){ + found=1; + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + if (found){ + /*We found an unstable constraint. : */ + return 0; + } + else{ + return 1; + } +} + +int SetPreStable(DataSet* loads){ + + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and set loads to pre stable.:*/ + for (i=0;iSize();i++){ + + if(RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + riftfront->SetPreStable(); + } + } +} + +int PreConstrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + int unstable; + int sum_num_unstable_constraints; + int num_unstable_constraints=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->PreConstrain(&unstable,analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; + +} + +int Constrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + int unstable; + int sum_num_unstable_constraints; + int num_unstable_constraints=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->Constrain(&unstable,analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; + +} + +void FreezeConstraints(DataSet* loads,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->FreezeConstraints(analysis_type); + + } + } + +} + +int IsFrozen(DataSet* loads){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + if (riftfront->IsFrozen()){ + found=1; + break; + } + } + } + + /*Is there just one found? that would mean we have frozen! : */ + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} + +int MaxPenetrationInInputs(DataSet* loads,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*rift penetration: */ + double max_penetration=0; + double mpi_max_penetration; + double penetration; + + /*Ok, we are going to find the grid pairs which are not penetrating, even though they + * are penalised. We will release only the one with has least <0 penetration. : */ + + max_penetration=0; + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->MaxPenetration(&penetration,analysis_type); + + if (penetration>max_penetration)max_penetration=penetration; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&max_penetration,&mpi_max_penetration,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_max_penetration,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + max_penetration=mpi_max_penetration; + #endif + + /*feed max_penetration to inputs: */ + loads->UpdateInputsFromVector(&max_penetration,MaxPenetrationEnum,ConstantEnum); +} + +int PotentialUnstableConstraints(DataSet* loads,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*Ok, we are going to find the grid pairs which are not penetrating, even though they + * are penalised. We will release only the one with has least <0 penetration. : */ + int unstable=0; + int sum_num_unstable_constraints=0; + int num_unstable_constraints=0; + + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->PotentialUnstableConstraint(&unstable,analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + return num_unstable_constraints; +} Index: /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.h =================================================================== --- /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.h (revision 3906) +++ /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.h (revision 3906) @@ -0,0 +1,33 @@ +/*!\file RiftConstraints.h + * \brief: manage penalties for rifts + */ + + +#ifndef _RIFTCONSTRAINTS_H_ +#define _RIFTCONSTRAINTS_H_ + +#include "../objects/objects.h" +#include "../DataSet/DataSet.h" + +int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,int min_mechanical_constraints,int analysis_type,int sub_analysis_type); + +int RiftIsPresent(DataSet* loads); + +int IsPreStable(DataSet* loads); + +int SetPreStable(DataSet* loads); + +int PreConstrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type_enum); + +int Constrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type_enum); + +void FreezeConstraints(DataSet* loads,int analysis_type); + +int MaxPenetrationInInputs(DataSet* loads,int analysis_type_enum); + +int PotentialUnstableConstraints(DataSet* loads,int analysis_type_enum); + +int IsMaterialStable(DataSet* loads,int analysis_type_enum); + +int IsFrozen(DataSet* loads); +#endif Index: /issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 3906) +++ /issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 3906) @@ -0,0 +1,48 @@ +/*!\file PenaltySystemMatricesx + * \brief: add penalties to system matrices + */ + +#include "./PenaltySystemMatricesx.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, Parameters* parameters, + int kflag,int pflag,int analysis_type,int sub_analysis_type){ + + int i; + + extern int num_procs; + extern int my_rank; + double kmax; + + /*First, get elements and loads configured: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + nodes->Configure(elements,loads, nodes,vertices, materials,parameters); + loads->Configure(elements, loads, nodes,vertices, materials,parameters); + parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */ + MatNorm(Kgg,NORM_INFINITY,&kmax); + + /*Add penalties to stiffnesses, from loads: */ + if(kflag)loads->PenaltyCreateKMatrix(Kgg,kmax,analysis_type,sub_analysis_type); + if(pflag)loads->PenaltyCreatePVector(pg,kmax,analysis_type,sub_analysis_type); + + /*Assemble matrices: */ + if(kflag){ + MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY); + MatCompress(Kgg); + } + if(pflag){ + VecAssemblyBegin(pg); + VecAssemblyEnd(pg); + } + + /*Assign output pointers:*/ + if(pkmax)*pkmax=kmax; + +} Index: /issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.h =================================================================== --- /issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 3906) +++ /issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 3906) @@ -0,0 +1,16 @@ +/*!\file: PenaltySystemMatricesx.h + * \brief + */ + +#ifndef _PENALTYSYSTEMMATRICESX_H +#define _PENALTYSYSTEMMATRICESX_H + +#include "../DataSet/DataSet.h" +#include "../objects/objects.h" + +/* local prototypes: */ +void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, Parameters* parameters, + int kflag,int pflag,int analysis_type,int sub_analysis_type); + +#endif /* _PENALTYSYSTEMMATRICESX_H */ + Index: /issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.cpp =================================================================== --- /issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.cpp (revision 3906) +++ /issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.cpp (revision 3906) @@ -0,0 +1,41 @@ +/*!\file ProcessParamsx + * \brief: process parameters using partitioning vector. + * Go through all parameters in the 'parameters' dataset. For each parameter that holds a doublevec (ie, a double* vector synchronized across + * the MPI ring of a cluster), partition the vector so that the new node partitioning decided in ModelProcessor is applied. Otherwise, + * parameters coming directly from Matlab would be serially partitioned, which means could not be recognized by each individual node or element. + * The partition needs to be the parallel partitionting. + */ + +#include "./ProcessParamsx.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void ProcessParamsx( Parameters* parameters, Vec part){ + + + int i; + double *partition = NULL; + int numberofvertices; + Param *param = NULL; + + /*Need number of vertices to repartition DoubleVecParam objects: */ + parameters->FindParam(&numberofvertices,NumberOfVerticesEnum); + + /*serialize partition vector: */ + if(part)VecToMPISerial(&partition,part); + + for(i=0;iSize();i++){ + + param=(Param*)parameters->GetObjectByOffset(i); + param->Process(partition,numberofvertices); + + } + + /*Free ressources:*/ + xfree((void**)&partition); + +} Index: /issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.h =================================================================== --- /issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.h (revision 3906) +++ /issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.h (revision 3906) @@ -0,0 +1,14 @@ +/*!\file: ProcessParamsx.h + * \brief header file for processing parameters + */ + +#ifndef _PROCESSPARAMSX_H +#define _PROCESSPARAMSX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void ProcessParamsx( Parameters* parameters, Vec partition); + +#endif /* _PROCESSPARAMSX_H */ + Index: /issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp =================================================================== --- /issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp (revision 3906) +++ /issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp (revision 3906) @@ -0,0 +1,38 @@ +/*!\file PropagateFlagsFromConnectivityx + */ + +#include "./PropagateFlagsFromConnectivityx.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" + +void RecursivePropagation(double* pool, double* connectivity,int index, double* flags); + +void PropagateFlagsFromConnectivityx( double* pool, double* connectivity,int index, double* flags){ + + /*Call recursive propagation routine: */ + RecursivePropagation(pool, connectivity,index, flags); +} + + +void RecursivePropagation(double* pool, double* connectivity, int index, double* flags){ + + int i; + int newel; + + /*if this element (index) belongs to the pool already, skip: */ + if(pool[index-1])return; + + /*if this element does not belong to the flags set, skip: */ + if(flags[index-1]==0)return; + + /*put this element (index), which belongs to the flags, into the pool: */ + pool[index-1]=1; + + /*now, propagate recursively using connectivity of this element: */ + for(i=0;i<3;i++){ + newel=(int)*(connectivity+(index-1)*3+i); + RecursivePropagation(pool, connectivity, newel, flags); + } +} Index: /issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h =================================================================== --- /issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h (revision 3906) +++ /issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h (revision 3906) @@ -0,0 +1,13 @@ +/*!\file: PropagateFlagsFromConnectivityx.h + */ + +#ifndef _PROPAGATEFLAGSFROMCONNECTIVITYX_H +#define _PROPAGATEFLAGSFROMCONNECTIVITYX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void PropagateFlagsFromConnectivityx( double* pool, double* connectivity,int index, double* flags); + +#endif /* _PROPAGATEFLAGSFROMCONNECTIVITYX_H */ + Index: /issm/trunk/src/c/modules/Qmux/DakotaResponses.cpp =================================================================== --- /issm/trunk/src/c/modules/Qmux/DakotaResponses.cpp (revision 3906) +++ /issm/trunk/src/c/modules/Qmux/DakotaResponses.cpp (revision 3906) @@ -0,0 +1,338 @@ +/*!\file: DakotaResponses.cpp + * \brief compute dakota responses, using a list of response descriptors. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../DataSet/DataSet.h" +#include "../shared/shared.h" +#include "../include/include.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../MassFluxx/MassFluxx.h" +#include "../Misfitx/Misfitx.h" +#include "../modules.h" + +void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model,DataSet* results,DataSet* processed_results,int analysis_type,int sub_analysis_type){ + + int i,j; + int found=0; + char* response_descriptor=NULL; + int numberofnodes; + extern int my_rank; + + /*some data needed across the responses: */ + model->FindParam(&numberofnodes,NumberOfNodesEnum); + + + for(i=0;iFindResult((void*)&vel,"vel"); + if(!found)ISSMERROR(" could not find vel to compute min_vel"); + + min_vel=vel[0]; + for(j=1;jFindResult((void*)&vel,"vel"); + if(!found)ISSMERROR(" could not find vel to compute max_vel"); + + max_vel=vel[0]; + for(j=1;jmax_vel)max_vel=vel[j]; + } + if(my_rank==0)responses[i]=max_vel; + + /*Free ressources:*/ + xfree((void**)&vel); + + } + else if(strcmp(response_descriptor,"min_vx")==0){ + double* vx=NULL; + double min_vx=0; + + found=processed_results->FindResult((void*)&vx,"vx"); + if(!found)ISSMERROR(" could not find vx to compute min_vx"); + + min_vx=vx[0]; + for(j=1;jFindResult((void*)&vx,"vx"); + if(!found)ISSMERROR(" could not find vx to compute max_vx"); + + max_vx=vx[0]; + for(j=1;jmax_vx)max_vx=vx[j]; + } + if(my_rank==0)responses[i]=max_vx; + + /*Free ressources:*/ + xfree((void**)&vx); + + } + else if(strcmp(response_descriptor,"max_abs_vx")==0){ + double* vx=NULL; + double max_abs_vx=0; + + found=processed_results->FindResult((void*)&vx,"vx"); + if(!found)ISSMERROR(" could not find vx to compute max_abs_vx"); + + max_abs_vx=fabs(vx[0]); + for(j=1;jmax_abs_vx)max_abs_vx=fabs(vx[j]); + } + if(my_rank==0)responses[i]=max_abs_vx; + + /*Free ressources:*/ + xfree((void**)&vx); + + } + else if(strcmp(response_descriptor,"min_vy")==0){ + double* vy=NULL; + double min_vy=0; + + found=processed_results->FindResult((void*)&vy,"vy"); + if(!found)ISSMERROR(" could not find vy to compute min_vy"); + + min_vy=vy[0]; + for(j=1;jFindResult((void*)&vy,"vy"); + if(!found)ISSMERROR(" could not find vy to compute max_vy"); + + max_vy=vy[0]; + for(j=1;jmax_vy)max_vy=vy[j]; + } + if(my_rank==0)responses[i]=max_vy; + + /*Free ressources:*/ + xfree((void**)&vy); + + } + else if(strcmp(response_descriptor,"max_abs_vy")==0){ + double* vy=NULL; + double max_abs_vy=0; + + found=processed_results->FindResult((void*)&vy,"vy"); + if(!found)ISSMERROR(" could not find vy to compute max_abs_vy"); + + max_abs_vy=fabs(vy[0]); + for(j=1;jmax_abs_vy)max_abs_vy=fabs(vy[j]); + } + if(my_rank==0)responses[i]=max_abs_vy; + + /*Free ressources:*/ + xfree((void**)&vy); + + } + else if(strcmp(response_descriptor,"min_vz")==0){ + double* vz=NULL; + double min_vz=0; + + found=processed_results->FindResult((void*)&vz,"vz"); + if(!found)ISSMERROR(" could not find vz to compute min_vz"); + + min_vz=vz[0]; + for(j=1;jFindResult((void*)&vz,"vz"); + if(!found)ISSMERROR(" could not find vz to compute max_vz"); + + max_vz=vz[0]; + for(j=1;jmax_vz)max_vz=vz[j]; + } + if(my_rank==0)responses[i]=max_vz; + + /*Free ressources:*/ + xfree((void**)&vz); + + } + else if(strcmp(response_descriptor,"max_abs_vz")==0){ + double* vz=NULL; + double max_abs_vz=0; + + found=processed_results->FindResult((void*)&vz,"vz"); + if(!found)ISSMERROR(" could not find vz to compute max_abs_vz"); + + max_abs_vz=fabs(vz[0]); + for(j=1;jmax_abs_vz)max_abs_vz=fabs(vz[j]); + } + if(my_rank==0)responses[i]=max_abs_vz; + + /*Free ressources:*/ + xfree((void**)&vz); + } + else if(strcmp(response_descriptor,"misfit")==0){ + + int isstokes,ismacayealpattyn,ishutter; + FemModel* femmodel=NULL; + double J=0; + int numberofdofspernode,numberofnodes; + Vec u_g=NULL; + double* u_g_double=NULL; + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + double* fit=NULL; + + /*retrieve active fem model: */ + model->FindParam(&isstokes,IsStokesEnum); + model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); + model->FindParam(&ishutter,IsHutterEnum); + + if(isstokes){ + femmodel=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + } + if(ismacayealpattyn){ + femmodel=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + } + if(ishutter){ + femmodel=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + } + + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + femmodel->parameters->FindParam(&fit,NULL,NULL,FitEnum); + + /*Recover velocity: */ + found=results->FindResult(&u_g,"u_g"); + VecToMPISerial(&u_g_double,u_g); + if(!found)ISSMERROR(" could not find velocity to compute misfit"); + + SplitSolutionVectorx(u_g,numberofnodes,numberofdofspernode,&vx,&vy,&vz); + + /*Add to inputs: */ + femmodel->elements->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + femmodel->elements->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + femmodel->elements->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + femmodel->elements->UpdateInputsFromVector(&fit[0],FitEnum,ConstantEnum); + + /*Compute misfit: */ + Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,analysis_type,sub_analysis_type); + + + if(my_rank==0)responses[i]=J; + + /*Some cleanup: */ + VecFree(&u_g); + xfree((void**)&u_g_double); + xfree((void**)&fit); + + } + else if(strcmp(response_descriptor,"mass_flux")==0){ + + int isstokes,ismacayealpattyn,ishutter; + Vec ug=NULL; + double* ug_serial=NULL; + double mass_flux=0; + double* segments=NULL; + int num_segments; + DataSet* elements=NULL; + DataSet* nodes=NULL; + DataSet* materials=NULL; + DataSet* parameters=NULL; + FemModel* femmodel=NULL; + Param* param=NULL; + + /*retrieve velocities: */ + found=results->FindResult(&ug,"u_g"); + if(!found)ISSMERROR(" could not find velocity to compute mass_flux"); + VecToMPISerial(&ug_serial,ug); + + /*retrieve active fem model: */ + model->FindParam(&isstokes,IsStokesEnum); + model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); + model->FindParam(&ishutter,IsHutterEnum); + + if(isstokes){ + femmodel=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + } + if(ismacayealpattyn){ + femmodel=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + } + if(ishutter){ + femmodel=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + } + + /*retrieve qmu_mass_flux_segments: */ + femmodel->parameters->FindParam(&segments,&num_segments,QmuMassFluxSegmentsEnum); + + /*call mass flux module: */ + MassFluxx(&mass_flux,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,segments,num_segments,ug_serial); + + if(my_rank==0)responses[i]=mass_flux; + + /*Free ressources:*/ + VecFree(&ug); + xfree((void**)&ug_serial); + xfree((void**)&segments); + } + else{ + if(my_rank==0)printf("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!"); + ISSMERROR("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!"); + } + } + +} Index: /issm/trunk/src/c/modules/Qmux/Qmux.cpp =================================================================== --- /issm/trunk/src/c/modules/Qmux/Qmux.cpp (revision 3906) +++ /issm/trunk/src/c/modules/Qmux/Qmux.cpp (revision 3906) @@ -0,0 +1,145 @@ +/*!\file: Qmux.cpp + * \brief: wrapper to the Dakota capabilities. qmu fires up Dakota, and registers a Dakota Pluggin + * which will be in charge of running the solution sequences repeteadly, to garner statistics. + * + * This routine deals with running ISSM and Dakota in library mode. In library mode, Dakota does not + * run as an execuatble. Its capabilities are linked into the ISSM software. ISSM calls dakota routines + * directly from the dakota library. qmu.cpp is the code that is in charge of calling those routines. + * + * Dakota has its own way of running in parallel (for embarassingly parallel jobs). We do not want that, + * as ISSM knows exactly how to run "really parallel" jobs that use all CPUS. To bypass Dakota's parallelism, + * we overloaded the constructor for the parallel library (see the Dakota patch in the externalpackages/dakota + * directory). This overloaded constructor fires up Dakota serially on CPU 0 only! We take care of broadcasting + * to the other CPUS, hence ISSM is running in parallel, and Dakota serially on CPU0. + * + * Now, how does CPU 0 drive all other CPUS to carry out sensitivity analysese? By synchronizing its call to + * our ISSM cores (diagnostic_core, thermal_core, transient_core, etc ...) on CPU 0 with all other CPUS. + * This explains the structure of qmu.cpp, where cpu 0 runs Dakota, the Dakota pluggin fires up SpawnCore.cpp, + * while the other CPUS are waiting for a broadcast from CPU0, once they get it, they also fire up + * SpawnCore. In the end, SpawnCore is fired up on all CPUS, with CPU0 having Dakota inputs, that it will + * broacast to other CPUS. + * + * Now, how does dakota call the SpawnCore routine? The SpawnCore is embedded into the DakotaPlugin object + * which is derived from the Direct Interface Dakota objct. This is the only way to run Dakota in library + * mode (see their developper guide for more info). Dakota registers the DakotaPlugin object into its own + * database, and calls on the embedded SpawnCore from CPU0. + * + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./Qmux.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../objects/DakotaPlugin.h" + +#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. +#include "ParallelLibrary.H" +#include "ProblemDescDB.H" +#include "DakotaStrategy.H" +#include "DakotaModel.H" +#include "DakotaInterface.H" + +#endif + +#ifdef _SERIAL_ +void Qmux(mxArray* model,int analysis_type,int sub_analysis_type,char* dakota_input_file,char* dakota_output_file,char* dakota_error_file){ +#else +void Qmux(Model* model,int analysis_type,int sub_analysis_type){ +#endif + + + #ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. + + #ifdef _PARALLEL_ + char* dakota_input_file=NULL; + char* dakota_output_file=NULL; + char* dakota_error_file=NULL; + extern int my_rank; + #endif + int status=0; + Dakota::ModelLIter ml_iter; + + #ifdef _PARALLEL_ + /*Recover dakota_input_file, dakota_output_file and dakota_error_file, in the parameters dataset in parallel */ + model->FindParam(&dakota_input_file,QmuInNameEnum); + model->FindParam(&dakota_output_file,QmuOutNameEnum); + model->FindParam(&dakota_error_file,QmuErrNameEnum); + #endif + + #ifdef _PARALLEL_ + if(my_rank==0){ + #endif + + // Instantiate/initialize the parallel library and problem description + // database objects. + #ifdef _SERIAL_ + Dakota::ParallelLibrary parallel_lib; //use Dakota's standard library mode constructor + #else + Dakota::ParallelLibrary parallel_lib("serial"); //use our own ISSM Dakota library mode constructor, which only fires up Dakota on CPU 0. + #endif + Dakota::ProblemDescDB problem_db(parallel_lib); + + // Manage input file parsing, output redirection, and restart processing + // without a CommandLineHandler. This version relies on parsing of an + // input file. + problem_db.manage_inputs(dakota_input_file); + // specify_outputs_restart() is only necessary if specifying non-defaults + parallel_lib.specify_outputs_restart(dakota_output_file,dakota_error_file,NULL,NULL); + + // Instantiate the Strategy object (which instantiates all Model and + // Iterator objects) using the parsed information in problem_db. + Dakota::Strategy selected_strategy(problem_db); + + // convenience function for iterating over models and performing any + // interface plug-ins + Dakota::ModelList& models = problem_db.model_list(); + + for (ml_iter = models.begin(); ml_iter != models.end(); ml_iter++) { + + Dakota::Interface& interface = ml_iter->interface(); + + //set DB nodes to the existing Model specification + problem_db.set_db_model_nodes(ml_iter->model_id()); + + // Serial case: plug in derived Interface object without an analysisComm + interface.assign_rep(new SIM::DakotaPlugin(problem_db,(void*)model,analysis_type,sub_analysis_type), false); + } + + // Execute the strategy + problem_db.lock(); // prevent run-time DB queries + selected_strategy.run_strategy(); + + #ifdef _PARALLEL_ + //Warn other cpus that we are done running the dakota iterator, by setting the counter to -1: + SpawnCore(NULL,0, NULL,NULL,0,model,analysis_type,sub_analysis_type,-1); + #endif + + #ifdef _PARALLEL_ + } + else{ + + for(;;){ + if(!SpawnCore(NULL,0, NULL,NULL,0,model,analysis_type,sub_analysis_type,0))break; //counter came in at -1 on cpu0, bail out. + } + } + #endif + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&dakota_input_file); + xfree((void**)&dakota_error_file); + xfree((void**)&dakota_output_file); + #endif + + + #endif + +} Index: /issm/trunk/src/c/modules/Qmux/Qmux.h =================================================================== --- /issm/trunk/src/c/modules/Qmux/Qmux.h (revision 3906) +++ /issm/trunk/src/c/modules/Qmux/Qmux.h (revision 3906) @@ -0,0 +1,24 @@ +/*!\file: Qmux.h + * \brief header file for Qmu engine using Dakota + */ + +#ifndef _QMUX_H +#define _QMUX_H + +#include "../DataSet/DataSet.h" +#include "../objects/objects.h" + +/* local prototypes: */ +int SpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* model,int analysis_type,int sub_analysis_type,int counter); +#ifdef _SERIAL_ +void Qmux(mxArray* model,int analysis_type,int sub_analysis_type,char* dakota_input_file,char* dakota_output_file,char* dakota_error_file); +void SpawnCoreSerial(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, mxArray* model,int analysis_type,int sub_analysis_type,int counter); +#else +void Qmux(Model* model,int analysis_type,int sub_analysis_type); +void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, Model* model,int analysis_type,int sub_analysis_type,int counter); +void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model, DataSet* results,DataSet* processed_results,int analysis_type,int sub_analysis_type); +#endif + + +#endif /* _QMUX_H */ + Index: /issm/trunk/src/c/modules/Qmux/SpawnCore.cpp =================================================================== --- /issm/trunk/src/c/modules/Qmux/SpawnCore.cpp (revision 3906) +++ /issm/trunk/src/c/modules/Qmux/SpawnCore.cpp (revision 3906) @@ -0,0 +1,32 @@ +/*!\file: SpawnCore.cpp + * \brief: branch into SpawnCoreMatlab and SpawnCoreParallel. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +#include "../objects/objects.h" +#include "../io/io.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "./Qmux.h" +#include "../include/include.h" + +int SpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* model,int analysis_type,int sub_analysis_type,int counter){ + + /*Branch into a serial SpawnCore and a parallel SpawnCore: */ + #ifdef _SERIAL_ + SpawnCoreSerial(responses, numresponses, variables, variables_descriptors,numvariables, (mxArray*)model, analysis_type,sub_analysis_type,counter); + #else + /*Call SpawnCoreParallel unless counter=-1 on cpu0, in which case, bail out and return 0: */ + MPI_Bcast(&counter,1,MPI_INT,0,MPI_COMM_WORLD); if(counter==-1)return 0; + + SpawnCoreParallel(responses, numresponses, variables, variables_descriptors,numvariables, (Model*)model, analysis_type,sub_analysis_type,counter); + #endif + + return 1; +} Index: /issm/trunk/src/c/modules/Qmux/SpawnCoreParallel.cpp =================================================================== --- /issm/trunk/src/c/modules/Qmux/SpawnCoreParallel.cpp (revision 3906) +++ /issm/trunk/src/c/modules/Qmux/SpawnCoreParallel.cpp (revision 3906) @@ -0,0 +1,158 @@ +/*!\file: SpawnCoreParallel.cpp + * \brief: run core ISSM solution using Dakota inputs coming from CPU 0. + * \sa qmu.cpp DakotaPlugin.cpp + * + * This routine needs to be understood simultaneously with qmu.cpp and DakotaPlugin. + * SpawnCoreParallel is called by all CPUS, with CPU 0 holding Dakota variable values, along + * with variable descriptors. + * + * SpawnCoreParallel takes care of broadcasting the variables and their descriptors across the MPI + * ring. Once this is done, we use the variables to modify the inputs for the solution core. + * For ex, if "rho_ice" is provided, for ex 920, we include "rho_ice" in the inputs, then + * call the core with the modified inputs. This is the way we get Dakota to explore the parameter + * spce of the core. + * + * Once the core is called, we process the results of the core, and using the processed results, + * we compute response functions. The responses are computed on all CPUS, but they are targeted + * for CPU 0, which will get these values back to the Dakota engine. + * + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +#include "../objects/objects.h" +#include "../io/io.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "./Qmux.h" +#include "../include/include.h" +#include "../solutions/solutions.h" + +void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, Model* model,int analysis_type,int sub_analysis_type,int counter){ + + int i; + + /*output from core solutions: */ + DataSet* results=NULL; + DataSet* processed_results=NULL; + + char** responses_descriptors=NULL; + int num_responses_descriptors; + char* string=NULL; + int string_length; + double* qmu_part=NULL; + int qmu_npart; + int verbose=0; + int dummy; + + extern int my_rank; + + + /*synchronize all cpus, as CPU 0 is probably late (it is starting the entire dakota strategy!) : */ + MPI_Barrier(MPI_COMM_WORLD); + + /*some parameters needed: */ + model->FindParam(&verbose,VerboseEnum); + + /*First off, recover the response descriptors for the response functions: */ + model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum)->parameters->FindParam(&responses_descriptors,&num_responses_descriptors,ResponseDescriptorsEnum); + + /*Recover partitioning for dakota: */ + model->FindParam(&qmu_npart,QmuNPartEnum); + model->FindParam(&qmu_part,&dummy,QmuPartEnum); + + /*broadcast variables: only cpu 0 has correct values*/ + MPI_Bcast(&numvariables,1,MPI_INT,0,MPI_COMM_WORLD); + if(my_rank!=0)variables=(double*)xmalloc(numvariables*sizeof(double)); + MPI_Bcast(variables,numvariables,MPI_DOUBLE,0,MPI_COMM_WORLD); + + /*broadcast variables_descriptors: */ + if(my_rank!=0){ + variables_descriptors=(char**)xmalloc(numvariables*sizeof(char*)); + } + for(i=0;iUpdateFromDakota(variables,variables_descriptors,numvariables,model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum)->parameters,qmu_part,qmu_npart); //diagnostic horiz model is the one holding the parameters for Dakota. + + /*Run the analysis core solution sequence: */ + if(analysis_type==DiagnosticAnalysisEnum){ + + if(verbose)_printf_("Starting diagnostic core\n"); + + results=diagnostic_core(model); + + } + else if(analysis_type==ThermalAnalysisEnum){ + + if(verbose)_printf_("Starting thermal core\n"); + results=thermal_core(model); + + } + else if(analysis_type==PrognosticAnalysisEnum){ + + if(verbose)_printf_("Starting prognostic core\n"); + results=prognostic_core(model); + + } + else if(analysis_type==TransientAnalysisEnum){ + + if(verbose)_printf_("Starting transient core\n"); + results=transient_core(model); + + } + else ISSMERROR("%s%i%s%i%s"," analysis_type ",analysis_type," and sub_analysis_type ",sub_analysis_type," not supported yet!"); + + + + /*Now process the outputs, before computing the dakota responses: */ + if(verbose)_printf_("process results:\n"); + + ProcessResults(&processed_results,results,model,analysis_type); + + /*compute responses on cpu 0: dummy for now! */ + if(verbose)_printf_("compute dakota responses:\n"); + DakotaResponses(responses,responses_descriptors,numresponses,model,results,processed_results,analysis_type,sub_analysis_type); + + /*Free ressources:*/ + delete results; + delete processed_results; + + //variables only on cpu != 0 + if(my_rank!=0){ + xfree((void**)&variables); + for(i=0;iGetMSize()){ + + VecPartition(&pn, pg, nodesets->GetPV_N(), nodesets->GetNSize()); + VecPartition(&pm, pg, nodesets->GetPV_M(), nodesets->GetMSize()); + + /*Create G_mn'*/ + #if _PETSC_VERSION_ == 2 + MatTranspose(G_mn,&tG_mn); + #else + MatTranspose(G_mn,MAT_INITIAL_MATRIX,&tG_mn); + #endif + + /*Create G_mn' * pm*/ + MatGetLocalSize(tG_mn,&tG_mnm,&tG_mnn); + tG_mnpm=NewVecFromLocalSize(tG_mnm); + MatMultPatch(tG_mn,pm,tG_mnpm); + + /*add pn to G_mn' * pm*/ + a=1; + VecAXPY(pn,a,tG_mnpm); + } + else{ + VecDuplicate(pg,&pn); + VecCopy(pg,pn); + } + + /*% Reduce pn to pf:*/ + + if (nodesets->GetSSize()){ + + VecPartition(&pf, pn, nodesets->GetPV_F(),nodesets->GetFSize()); + + /*pf = pf - Kfs * y_s;*/ + MatGetLocalSize(Kfs,&Kfsm,&Kfsn); + Kfsy_s=NewVecFromLocalSize(Kfsm); + MatMultPatch(Kfs,y_s,Kfsy_s); + + a=-1; + VecAXPY(pf,a,Kfsy_s); + } + else{ + /*pf=pn*/ + VecDuplicate(pn,&pf); + VecCopy(pn,pf); + } + } + + + /*Assign correct pointer*/ + *ppf=pf; + + /*Free ressources and return*/ + VecFree(&pn); + VecFree(&pm); + MatFree(&tG_mn); + VecFree(&tG_mnpm); + VecFree(&Kfsy_s); + +} Index: /issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h =================================================================== --- /issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h (revision 3906) +++ /issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h (revision 3906) @@ -0,0 +1,14 @@ +/*!\file: Reduceloadfromgtofx.h + * \brief reduce load from g set to f set + */ + +#ifndef _REDUCELOADFROMGTOFX_H +#define _REDUCELOADFROMGTOFX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +void Reduceloadfromgtofx( Vec* ppf, Vec pg, Mat Gmn, Mat Kfs, Vec ys, NodeSets* nodesets); + +#endif /* _REDUCELOADFROMGTOFX_H */ +