Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 585) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 586) @@ -1229,5 +1229,53 @@ } - +void DataSet::ThicknessExtrude(Vec tg,double* tg_serial){ + + vector::iterator object; + Penta* penta=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PentaEnum()){ + + penta=(Penta*)(*object); + penta->ThicknessExtrude(tg,tg_serial); + + } + } + +} +void DataSet::VelocityExtrudeAllElements(Vec ug,double* ug_serial){ + + vector::iterator object; + Penta* penta=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PentaEnum()){ + + penta=(Penta*)(*object); + penta->VelocityExtrudeAllElements(ug,ug_serial); + + } + } + +} +void DataSet::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){ + + + vector::iterator object; + Penta* penta=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PentaEnum()){ + + penta=(Penta*)(*object); + penta->VelocityDepthAverageAtBase(ug,ug_serial); + + } + } + +} void DataSet::SlopeExtrude(Vec sg,double* sg_serial){ Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 585) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 586) @@ -75,4 +75,7 @@ void Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type); void VelocityExtrude(Vec ug,double* ug_serial); + void ThicknessExtrude(Vec ug,double* ug_serial); + void VelocityExtrudeAllElements(Vec ug,double* ug_serial); + void VelocityDepthAverageAtBase(Vec ug,double* ug_serial); void SlopeExtrude(Vec sg,double* sg_serial); int DeleteObject(Object* object); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 585) +++ /issm/trunk/src/c/Makefile.am (revision 586) @@ -1,3 +1,3 @@ -INCLUDES = @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@ @METISINCL@ @PLAPACKINCL@ @BLASLAPACKINCL@ @MUMPSINCL@ @TRIANGLEINCL@ +INCLUDES = @DAKOTAINCL@ @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@ @METISINCL@ @PLAPACKINCL@ @BLASLAPACKINCL@ @MUMPSINCL@ @TRIANGLEINCL@ #Compile serial library, and then try and compile parallel library @@ -126,4 +126,5 @@ ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\ ./toolkits/petsc/patches/NewMat.cpp\ + ./toolkits/petsc/patches/SerialToVec.cpp\ ./toolkits/petsc/patches/VecFree.cpp\ ./toolkits/petsc/patches/VecDuplicatePatch.cpp\ @@ -196,4 +197,8 @@ ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\ ./ModelProcessorx/Melting/CreateParametersMelting.cpp\ + ./ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp\ + ./ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp\ + ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\ + ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\ ./Dofx/Dofx.h\ ./Dofx/Dofx.cpp\ @@ -252,6 +257,10 @@ ./ProcessParamsx/ProcessParamsx.cpp\ ./ProcessParamsx/ProcessParamsx.h\ + ./ThicknessExtrudex/ThicknessExtrudex.cpp\ + ./ThicknessExtrudex/ThicknessExtrudex.h\ ./VelocityExtrudex/VelocityExtrudex.cpp\ ./VelocityExtrudex/VelocityExtrudex.h\ + ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\ + ./VelocityDepthAveragex/VelocityDepthAveragex.h\ ./SlopeExtrudex/SlopeExtrudex.cpp\ ./SlopeExtrudex/SlopeExtrudex.h @@ -288,4 +297,6 @@ ./objects/Node.h\ ./objects/Node.cpp\ + ./objects/DakotaPlugin.h\ + ./objects/DakotaPlugin.cpp\ ./objects/Tria.h\ ./objects/Tria.cpp\ @@ -380,4 +391,5 @@ ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\ ./toolkits/petsc/patches/NewMat.cpp\ + ./toolkits/petsc/patches/SerialToVec.cpp\ ./toolkits/petsc/patches/VecFree.cpp\ ./toolkits/petsc/patches/VecDuplicatePatch.cpp\ @@ -450,4 +462,8 @@ ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\ ./ModelProcessorx/Melting/CreateParametersMelting.cpp\ + ./ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp\ + ./ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp\ + ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\ + ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\ ./Dofx/Dofx.h\ ./Dofx/Dofx.cpp\ @@ -502,4 +518,6 @@ ./Solverx/Solverx.cpp\ ./Solverx/Solverx.h\ + ./ThicknessExtrudex/ThicknessExtrudex.cpp\ + ./ThicknessExtrudex/ThicknessExtrudex.h\ ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\ ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\ @@ -508,4 +526,6 @@ ./VelocityExtrudex/VelocityExtrudex.cpp\ ./VelocityExtrudex/VelocityExtrudex.h\ + ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\ + ./VelocityDepthAveragex/VelocityDepthAveragex.h\ ./SlopeExtrudex/SlopeExtrudex.cpp\ ./SlopeExtrudex/SlopeExtrudex.h\ @@ -520,6 +540,10 @@ ./parallel/OutputControl.cpp\ ./parallel/OutputThermal.cpp\ + ./parallel/OutputPrognostic.cpp\ ./parallel/objectivefunctionC.cpp\ - ./parallel/GradJCompute.cpp + ./parallel/GradJCompute.cpp\ + ./parallel/qmu.cpp\ + ./parallel/SpawnCore.cpp\ + ./parallel/DakotaResponses.cpp libpISSM_a_CXXFLAGS = -fPIC -D_PARALLEL_ -D_C_ @@ -528,8 +552,9 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe control.exe thermal.exe +bin_PROGRAMS = diagnostic.exe +#control.exe thermal.exe prognostic.exe endif -LDADD = ./libpISSM.a $(TRIANGLELIB) $(METISLIB) $(PETSCLIB) $(SLEPCLIB) $(MUMPSLIB) $(PLAPACKLIB) $(MPILIB) $(X_LIBS) -lX11 $(BLASLAPACKLIB) $(SCALAPACKLIB) $(BLACSLIB) $(FLIBS) +LDADD = ./libpISSM.a $(TRIANGLELIB) $(METISLIB) $(PETSCLIB) $(DAKOTALIB) $(SLEPCLIB) $(MUMPSLIB) $(PLAPACKLIB) $(MPILIB) $(X_LIBS) -lX11 $(BLASLAPACKLIB) $(SCALAPACKLIB) $(BLACSLIB) $(FLIBS) diagnostic_exe_SOURCES = parallel/diagnostic.cpp @@ -541,2 +566,6 @@ thermal_exe_SOURCES = parallel/thermal.cpp thermal_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + +prognostic_exe_SOURCES = parallel/prognostic.cpp +prognostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 586) @@ -87,8 +87,10 @@ } else if (strcmp(model->analysis_type,"prognostic")==0){ - - //CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, model,model_handle); - //CreateConstraintsPrognostic(pconstraints,model,model_handle); - //CreateLoadsPrognostic(ploads,model,model_handle); + + CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsPrognostic(pconstraints,model,model_handle); + CreateLoadsPrognostic(ploads,model,model_handle); + CreateParametersPrognostic(pparameters,model,model_handle); + } else{ Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 586) @@ -24,4 +24,7 @@ int numberofdofspernode; int dim; + double* epart=NULL; + double* part=NULL; + /*Initialize dataset: */ @@ -42,4 +45,15 @@ parameters->AddObject(param); + //qmu analysis? + count++; + param= new Param(count,"qmu_analysis",INTEGER); + param->SetInteger(model->qmu_analysis); + parameters->AddObject(param); + + count++; + param= new Param(count,"qmu_npart",INTEGER); + param->SetInteger(model->qmu_npart); + parameters->AddObject(param); + //dimension 2d or 3d: if (strcmp(model->meshtype,"2d")==0)dim=2; @@ -201,4 +215,59 @@ parameters->AddObject(param); + /*Deal with responses and partition for qmu modeling: */ + if(model->qmu_analysis){ + char** descriptors=NULL; + char* descriptor=NULL; + char* tag=NULL; + + descriptors=(char**)xmalloc(model->numberofresponses*sizeof(char*)); + tag=(char*)xmalloc((strlen("descriptori")+1)*sizeof(char)); + + /*Fetch descriptors: */ + for(i=0;inumberofresponses;i++){ + sprintf(tag,"%s%i","descriptor",i); + ModelFetchData((void**)&descriptor,NULL,NULL,model_handle,tag,"String",NULL); + descriptors[i]=descriptor; + } + + /*Ok, we have all the descriptors. Build a parameter with it: */ + count++; + param= new Param(count,"descriptors",STRINGARRAY); + param->SetStringArray(descriptors,model->numberofresponses); + parameters->AddObject(param); + + /*Free data: */ + xfree((void**)&tag); + for(i=0;inumberofresponses;i++){ + char* descriptor=descriptors[i]; + xfree((void**)&descriptor); + } + xfree((void**)&descriptors); + + #ifdef _PARALLEL_ + /*partition grids in model->qmu_npart parts: */ + + if(strcmp(model->meshtype,"2d")==0){ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + } + else{ + ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); + } + + MeshPartitionx(&epart, &part,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,model->qmu_npart); + + count++; + param= new Param(count,"qmu_part",DOUBLEVEC); + param->SetDoubleVec(part,model->numberofnodes); + parameters->AddObject(param); + + /*Free elements and elements2d: */ + xfree((void**)&model->elements); + xfree((void**)&model->elements2d); + xfree((void**)&epart); + xfree((void**)&part); + + #endif + } Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 586) @@ -84,9 +84,9 @@ - cleanup_and_return: /*Free data: */ xfree((void**)&gridondirichlet_diag); xfree((void**)&dirichletvalues_diag); + cleanup_and_return: /*Assign output pointer: */ *pconstraints=constraints; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 586) @@ -69,6 +69,7 @@ double tria_meanvel;/*!scaling ratio for velocities*/ double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ - double tria_viscosity_overshoot; - + double tria_viscosity_overshoot; + int tria_artdiff; + /*matice constructor input: */ int matice_mid; @@ -294,5 +295,5 @@ /*Create tria element using its constructor:*/ - tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff); /*Add tria element to elements dataset: */ @@ -653,10 +654,4 @@ xfree((void**)&model->uppernodes); - cleanup_and_return: - - /*Assign output pointer: */ - *pelements=elements; - *pnodes=nodes; - *pmaterials=materials; /*Keep partitioning information into model*/ @@ -672,3 +667,10 @@ #endif + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + } Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 586) @@ -265,5 +265,4 @@ loads->Presort(); - cleanup_and_return: /*Free ressources:*/ @@ -276,5 +275,6 @@ xfree((void**)&riftsfill); xfree((void**)&riftsfriction); - + + cleanup_and_return: /*Assign output pointer: */ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 586) @@ -357,4 +357,9 @@ xfree((void**)&model->uppernodes); + + /*Keep partitioning information into model*/ + xfree((void**)&epart); + model->npart=npart; + cleanup_and_return: @@ -364,6 +369,3 @@ *pmaterials=materials; - /*Keep partitioning information into model*/ - xfree((void**)&epart); - model->npart=npart; } Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp (revision 586) @@ -88,7 +88,8 @@ constraints->Presort(); - cleanup_and_return: /*Free data: */ xfree((void**)&gridonstokes); + + cleanup_and_return: /*Assign output pointer: */ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 586) @@ -517,10 +517,4 @@ xfree((void**)&model->borderstokes); - cleanup_and_return: - - /*Assign output pointer: */ - *pelements=elements; - *pnodes=nodes; - *pmaterials=materials; /*Keep partitioning information into model*/ @@ -536,3 +530,10 @@ #endif + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + } Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 586) @@ -146,10 +146,4 @@ sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); - /*Width of elements: */ - if(strcmp(model->meshtype,"2d")==0){ - throw ErrorException(__FUNCT__," 2d mesh not supported for vertical velocity"); - } - - #ifdef _PARALLEL_ /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ @@ -420,4 +414,17 @@ xfree((void**)&model->uppernodes); + + /*Keep partitioning information into model*/ + model->epart=epart; + model->my_grids=my_grids; + model->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + xfree((void**)&npart); + VecFree(&gridborder); + #endif + cleanup_and_return: @@ -427,15 +434,3 @@ *pmaterials=materials; - /*Keep partitioning information into model*/ - model->epart=epart; - model->my_grids=my_grids; - model->my_bordergrids=my_bordergrids; - - /*Free ressources:*/ - #ifdef _PARALLEL_ - xfree((void**)&all_numgrids); - xfree((void**)&npart); - VecFree(&gridborder); - #endif - } Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 586) @@ -476,4 +476,17 @@ xfree((void**)&model->uppernodes); + + /*Keep partitioning information into model*/ + model->epart=epart; + model->my_grids=my_grids; + model->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + xfree((void**)&npart); + VecFree(&gridborder); + #endif + cleanup_and_return: @@ -483,15 +496,3 @@ *pmaterials=materials; - /*Keep partitioning information into model*/ - model->epart=epart; - model->my_grids=my_grids; - model->my_bordergrids=my_bordergrids; - - /*Free ressources:*/ - #ifdef _PARALLEL_ - xfree((void**)&all_numgrids); - xfree((void**)&npart); - VecFree(&gridborder); - #endif - } Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 586) @@ -32,10 +32,12 @@ /*!initialize all pointers to 0: */ - - model->repository=NULL; model->meshtype=NULL; model->analysis_type=NULL; model->sub_analysis_type=NULL; + model->qmu_analysis=0; + model->solverstring=NULL; + model->numberofresponses=0; + model->qmu_npart=0; model->numberofelements=0; model->numberofnodes=0; @@ -75,4 +77,5 @@ model->drag=NULL; model->p=NULL; + model->q=NULL; @@ -100,6 +103,6 @@ model->rho_water=0; model->rho_ice=0; - model->g=0; - model->n=0; + model->g=0; + model->n=NULL; model->B=NULL; @@ -165,4 +168,10 @@ model->ismacayealpattyn=0; model->isstokes=0; + + + model->epart=NULL; + model->npart=NULL; + model->my_grids=NULL; + model->my_bordergrids=NULL; return model; @@ -295,6 +304,8 @@ /*In ModelInit, we get all the data that is not difficult to get, and that is small: */ + ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL); ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL); + ModelFetchData((void**)&model->qmu_analysis,NULL,NULL,model_handle,"qmu_analysis","Integer",NULL); ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL); @@ -302,5 +313,4 @@ ModelFetchData((void**)&model->numberofnodes,NULL,NULL,model_handle,"numberofgrids","Integer",NULL); ModelFetchData((void**)&model->numberofelements,NULL,NULL,model_handle,"numberofelements","Integer",NULL); - /*!In case we are running 3d, we are going to need the collapsed and non-collapsed 2d meshes, from which the 3d mesh was extruded: */ if (strcmp(model->meshtype,"3d")==0){ @@ -311,4 +321,5 @@ ModelFetchData((void**)&model->numlayers,NULL,NULL,model_handle,"numlayers","Integer",NULL); } + /*elements type: */ @@ -367,4 +378,10 @@ ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL); + /*qmu: */ + if(model->qmu_analysis){ + ModelFetchData((void**)&model->numberofresponses,NULL,NULL,model_handle,"numberofresponses","Integer",NULL); + ModelFetchData((void**)&model->qmu_npart,NULL,NULL,model_handle,"npart","Integer",NULL); + } + /*Assign output pointers: */ *pmodel=model; Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 586) @@ -17,4 +17,5 @@ char* analysis_type; char* sub_analysis_type; + int qmu_analysis; char* solverstring; @@ -52,4 +53,8 @@ double* vx_obs; double* vy_obs; + + /*qmu: */ + int numberofresponses; + int qmu_npart; /*geometry: */ @@ -221,6 +226,12 @@ void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); void CreateConstraintsMelting(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); - void CreateLoadsMelting(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle); + void CreateLoadsMelting(DataSet** ploads, Model* model, ConstDataHandle model_handle); + void CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle); + + /*prognostic:*/ + void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); + void CreateConstraintsPrognostic(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); + void CreateLoadsPrognostic(DataSet** ploads, Model* model, ConstDataHandle model_handle); + void CreateParametersPrognostic(DataSet** pparameters,Model* model,ConstDataHandle model_handle); Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 586) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 586) @@ -0,0 +1,81 @@ +/* + * CreateConstraintsPrognostic.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateConstraintsPrognostic" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + + +void CreateConstraintsPrognostic(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ + + + int i; + int count; + + DataSet* constraints = NULL; + Spc* spc = NULL; + + /*spc intermediary data: */ + int spc_sid; + int spc_node; + int spc_dof; + double spc_value; + + double* dirichletvalues_prog=NULL; + double* gridondirichlet_prog=NULL; + + /*Create constraints: */ + constraints = new DataSet(ConstraintsEnum()); + + /*Fetch data: */ + ModelFetchData((void**)&gridondirichlet_prog,NULL,NULL,model_handle,"gridondirichlet_prog","Matrix","Mat"); + ModelFetchData((void**)&dirichletvalues_prog,NULL,NULL,model_handle,"dirichletvalues_prog","Matrix","Mat"); + + count=0; + + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((model->my_grids[i]==1)){ + #endif + + if ((int)gridondirichlet_prog[i]){ + + /*This grid needs to be spc'd: */ + + spc_sid=count; + spc_node=i+1; + spc_dof=1; //we enforce first translation degree of freedom, for temperature + spc_value=dirichletvalues_prog[i]; + + spc = new Spc(spc_sid,spc_node,spc_dof,spc_value); + constraints->AddObject(spc); + count++; + } + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + constraints->Presort(); + + /*Free data: */ + xfree((void**)&gridondirichlet_prog); + xfree((void**)&dirichletvalues_prog); + + cleanup_and_return: + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 586) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 586) @@ -0,0 +1,546 @@ +/* + * CreateElementsNodesAndMaterialsPrognostic.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterialsPrognostic" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../Model.h" + + +void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ + + + /*output: int* epart, int* my_grids, double* my_bordergrids*/ + + + int i,j,k,n; + extern int my_rank; + extern int num_procs; + + /*DataSets: */ + DataSet* elements = NULL; + DataSet* nodes = NULL; + DataSet* materials = NULL; + + /*Objects: */ + Node* node = NULL; + Tria* tria = NULL; + Penta* penta = NULL; + Matice* matice = NULL; + Matpar* matpar = NULL; + + int analysis_type; + int sub_analysis_type; + + /*output: */ + int* epart=NULL; //element partitioning. + int* npart=NULL; //node partitioning. + int* my_grids=NULL; + double* my_bordergrids=NULL; + + + /*intermediary: */ + int elements_width; //size of elements + double B_avg; + + /*tria constructor input: */ + int tria_id; + int tria_mid; + int tria_mparid; + int tria_g[3]; + double tria_h[3]; + double tria_s[3]; + double tria_b[3]; + double tria_k[3]; + double tria_melting[3]; + double tria_accumulation[3]; + double tria_geothermalflux[3]; + int tria_friction_type; + double tria_p; + double tria_q; + int tria_shelf; + double tria_meanvel;/*!scaling ratio for velocities*/ + double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ + double tria_viscosity_overshoot; + int tria_artdiff; + + /*matice constructor input: */ + int matice_mid; + double matice_B; + double matice_n; + + /*penta constructor input: */ + + int penta_id; + int penta_mid; + int penta_mparid; + int penta_g[6]; + double penta_h[6]; + double penta_s[6]; + double penta_b[6]; + double penta_k[6]; + int penta_friction_type; + double penta_p; + double penta_q; + int penta_shelf; + int penta_onbed; + int penta_onsurface; + double penta_meanvel;/*!scaling ratio for velocities*/ + double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ + int penta_collapse; + double penta_melting[6]; + double penta_accumulation[6]; + double penta_geothermalflux[6]; + int penta_artdiff; + int penta_thermal_steadystate; + double penta_viscosity_overshoot; + double penta_stokesreconditioning; + + /*matpar constructor input: */ + int matpar_mid; + double matpar_rho_ice; + double matpar_rho_water; + double matpar_heatcapacity; + double matpar_thermalconductivity; + double matpar_latentheat; + double matpar_beta; + double matpar_meltingpoint; + double matpar_mixed_layer_capacity; + double matpar_thermal_exchange_velocity; + double matpar_g; + + /* node constructor input: */ + int node_id; + int node_partitionborder=0; + double node_x[3]; + int node_onbed; + int node_onsurface; + int node_upper_node_id; + int node_numdofs; + + + #ifdef _PARALLEL_ + /*Metis partitioning: */ + int range; + Vec gridborder=NULL; + int my_numgrids; + int* all_numgrids=NULL; + int gridcount; + int count; + #endif + int first_grid_index; + + /*Rifts:*/ + int* riftsnumpenaltypairs; + double** riftspenaltypairs; + int* riftsfill; + double* riftsfriction; + double* riftpenaltypairs=NULL; + int el1,el2; + + /*First create the elements, nodes and material properties: */ + elements = new DataSet(ElementsEnum()); + nodes = new DataSet(NodesEnum()); + materials = new DataSet(MaterialsEnum()); + + + /*Get analysis_type: */ + analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); + + + /*Width of elements: */ + if(strcmp(model->meshtype,"2d")==0){ + elements_width=3; //tria elements + } + else{ + elements_width=6; //penta elements + } + + #ifdef _PARALLEL_ + /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ + if(strcmp(model->meshtype,"2d")==0){ + /*load elements: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + } + else{ + /*load elements2d: */ + ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); + } + + + MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&model->elements); + xfree((void**)&model->elements2d); + + + /*Deal with rifts, they have to be included into one partition only, not several: */ + FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); + + for(i=0;inumrifts;i++){ + riftpenaltypairs=model->riftspenaltypairs[i]; + for(j=0;jriftsnumpenaltypairs[i];j++){ + el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing + el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing + epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; + } + } + /*Free rifts: */ + xfree((void**)&riftsnumpenaltypairs); + for(i=0;inumrifts;i++){ + double* temp=riftspenaltypairs[i]; + xfree((void**)&temp); + } + xfree((void**)&riftspenaltypairs); + xfree((void**)&riftsfill); + xfree((void**)&riftsfriction); + + /*Used later on: */ + my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int)); + #endif + + + + /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ + + /*2d mesh: */ + if (strcmp(model->meshtype,"2d")==0){ + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); + + for (i=0;inumberofelements;i++){ + + #ifdef _PARALLEL_ + /*!All elements have been partitioned above, only create elements for this CPU: */ + if(my_rank==epart[i]){ + #endif + + + /*ids: */ + tria_id=i+1; //matlab indexing. + tria_mid=-1; //no need for materials + tria_mparid=-1; //no need for materials + + /*vertices offsets: */ + tria_g[0]=(int)*(model->elements+elements_width*i+0); + tria_g[1]=(int)*(model->elements+elements_width*i+1); + tria_g[2]=(int)*(model->elements+elements_width*i+2); + + /*thickness,surface and bed:*/ + tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. + tria_h[1]=*(model->thickness+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_h[2]=*(model->thickness+ ((int)*(model->elements+elements_width*i+2)-1)) ; + + tria_s[0]=*(model->surface+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_s[1]=*(model->surface+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_s[2]=*(model->surface+ ((int)*(model->elements+elements_width*i+2)-1)); + + tria_b[0]=*(model->bed+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_b[1]=*(model->bed+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_b[2]=*(model->bed+ ((int)*(model->elements+elements_width*i+2)-1)); + + /*element on iceshelf?:*/ + tria_shelf=(int)*(model->elementoniceshelf+i); + tria_artdiff=model->artificial_diffusivity; + + /*Create tria element using its constructor:*/ + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff); + + /*Add tria element to elements dataset: */ + elements->AddObject(tria); + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(model->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&model->thickness); + xfree((void**)&model->surface); + xfree((void**)&model->bed); + xfree((void**)&model->elementoniceshelf); + + } + else{ // if (strcmp(meshtype,"2d")==0) + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); + ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat"); + ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat"); + + for (i=0;inumberofelements;i++){ + #ifdef _PARALLEL_ + /*We are using our element partition to decide which elements will be created on this node: */ + if(my_rank==epart[i]){ + #endif + + + /*name and id: */ + penta_id=i+1; //matlab indexing. + penta_mid=-1; + penta_mparid=-1; //no need for materials + + /*vertices,thickness,surface,bed and drag: */ + for(j=0;j<6;j++){ + penta_g[j]=(int)*(model->elements+elements_width*i+j); + penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1)); + } + + /*diverse: */ + penta_shelf=(int)*(model->elementoniceshelf+i); + penta_onbed=(int)*(model->elementonbed+i); + penta_onsurface=(int)*(model->elementonsurface+i); + penta_collapse=1; + penta_artdiff=model->artificial_diffusivity; + + + /*Create Penta using its constructor:*/ + penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, + penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel, + penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff, + penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning); + + /*Add penta element to elements dataset: */ + elements->AddObject(penta); + + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(model->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+3)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+4)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+5)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&model->thickness); + xfree((void**)&model->surface); + xfree((void**)&model->bed); + xfree((void**)&model->elementoniceshelf); + xfree((void**)&model->elementonbed); + xfree((void**)&model->elementonsurface); + + } //if (strcmp(meshtype,"2d")==0) + + #ifdef _PARALLEL_ + /*From the element partitioning, we can determine which grids are on the inside of this cpu's + *element partition, and which are on its border with other nodes:*/ + gridborder=NewVec(model->numberofnodes); + + for (i=0;inumberofnodes;i++){ + if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + } + VecAssemblyBegin(gridborder); + VecAssemblyEnd(gridborder); + + #ifdef _DEBUG_ + VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); + #endif + + VecToMPISerial(&my_bordergrids,gridborder); + + #ifdef _DEBUG_ + if(my_rank==0){ + for (i=0;inumberofnodes;i++){ + printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); + } + } + #endif + #endif + + /*Partition penalties in 3d: */ + if(strcmp(model->meshtype,"3d")==0){ + + /*Get penalties: */ + ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat"); + + if(model->numpenalties){ + + model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int)); + #ifdef _SERIAL_ + for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=1; + #else + for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=-1; + + for(i=0;inumpenalties;i++){ + first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1); + if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids + /*All grids that are being penalised belong to this node's internal grid partition.:*/ + model->penaltypartitioning[i]=1; + } + if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border + model->penaltypartitioning[i]=0; + } + } + #endif + } + + /*Free penalties: */ + xfree((void**)&model->penalties); + } + + /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. + We can therefore determine which grids are internal to this node's partition + and which ones are shared with other nodes because they are on the border of this node's partition. Knowing + that, go and create the grids*/ + + /*Create nodes from x,y,z, as well as the spc values on those grids: */ + + /*First fetch data: */ + if (strcmp(model->meshtype,"3d")==0){ + ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat"); + ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat"); + } + ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat"); + ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat"); + ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat"); + ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat"); + ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat"); + + + /*Get number of dofs per node: */ + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); + + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((my_grids[i]==1)){ + #endif + + node_id=i+1; //matlab indexing + + + + #ifdef _PARALLEL_ + if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border + node_partitionborder=1; + } + else{ + node_partitionborder=0; + } + #else + node_partitionborder=0; + #endif + + + node_x[0]=model->x[i]; + node_x[1]=model->y[i]; + node_x[2]=model->z[i]; + + + node_onbed=(int)model->gridonbed[i]; + node_onsurface=(int)model->gridonsurface[i]; + if (strcmp(model->meshtype,"3d")==0){ + if (isnan(model->uppernodes[i])){ + node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + node_upper_node_id=(int)model->uppernodes[i]; + } + } + else{ + /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ + node_upper_node_id=node_id; + } + + /*Create node using its constructor: */ + node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id); + + /*set single point constraints.: */ + if (strcmp(model->meshtype,"3d")==0){ + /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ + if (!model->gridonbed[i]){ + for(k=1;k<=node_numdofs;k++){ + node->FreezeDof(k); + } + } + } + /*Add node to nodes dataset: */ + nodes->AddObject(node); + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + elements->Presort(); + nodes->Presort(); + materials->Presort(); + + /*Clean fetched data: */ + xfree((void**)&model->deadgrids); + xfree((void**)&model->x); + xfree((void**)&model->y); + xfree((void**)&model->z); + xfree((void**)&model->gridonbed); + xfree((void**)&model->gridonsurface); + xfree((void**)&model->uppernodes); + + + /*Keep partitioning information into model*/ + model->epart=epart; + model->my_grids=my_grids; + model->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + xfree((void**)&npart); + VecFree(&gridborder); + #endif + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + +} Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp (revision 586) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp (revision 586) @@ -0,0 +1,35 @@ +/*! \file CreateLoadsPrognostic.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateLoadsPrognostic" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../Model.h" + + +void CreateLoadsPrognostic(DataSet** ploads, Model* model,ConstDataHandle model_handle){ + + DataSet* loads = NULL; + + /*Create loads: */ + loads = new DataSet(LoadsEnum()); + + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + loads->Presort(); + + cleanup_and_return: + + /*Assign output pointer: */ + *ploads=loads; + +} + + Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp (revision 586) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp (revision 586) @@ -0,0 +1,111 @@ +/*!\file: CreateParametersPrognostic.cpp + * \brief driver for creating parameters dataset, for prognostic analysis. + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateParametersPrognostic" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + +void CreateParametersPrognostic(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ + + Param* param = NULL; + DataSet* parameters=NULL; + int count; + int i; + + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + double* pressure=NULL; + double* thickness=NULL; + double* accumulation=NULL; + double* melting=NULL; + + /*recover parameters : */ + parameters=*pparameters; + + count=parameters->Size(); + + + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!model->ismacayealpattyn)return; + + /*Get vx and vy: */ + ModelFetchData((void**)&vx,NULL,NULL,model_handle,"vx","Matrix","Mat"); + ModelFetchData((void**)&vy,NULL,NULL,model_handle,"vy","Matrix","Mat"); + ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat"); + + if(vx) for(i=0;inumberofnodes;i++)vx[i]=vx[i]/model->yts; + if(vy) for(i=0;inumberofnodes;i++)vy[i]=vy[i]/model->yts; + if(vz) for(i=0;inumberofnodes;i++)vz[i]=vz[i]/model->yts; + + count++; + param= new Param(count,"vx",DOUBLEVEC); + if(vx) param->SetDoubleVec(vx,model->numberofnodes); + else param->SetDoubleVec(vx,0); + parameters->AddObject(param); + + count++; + param= new Param(count,"vy",DOUBLEVEC); + if(vy) param->SetDoubleVec(vy,model->numberofnodes); + else param->SetDoubleVec(vy,0); + parameters->AddObject(param); + + count++; + param= new Param(count,"vz",DOUBLEVEC); + if(vz) param->SetDoubleVec(vz,model->numberofnodes); + else param->SetDoubleVec(vz,0); + parameters->AddObject(param); + + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + + /*Get thickness: */ + ModelFetchData((void**)&thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + + count++; + param= new Param(count,"thickness",DOUBLEVEC); + if(thickness) param->SetDoubleVec(thickness,model->numberofnodes); + else param->SetDoubleVec(thickness,0); + parameters->AddObject(param); + + /*Free thickness: */ + xfree((void**)&thickness); + + /*Get melting: */ + ModelFetchData((void**)&melting,NULL,NULL,model_handle,"melting","Matrix","Mat"); + + count++; + param= new Param(count,"melting",DOUBLEVEC); + if(melting) param->SetDoubleVec(melting,model->numberofnodes); + else param->SetDoubleVec(melting,0); + parameters->AddObject(param); + + /*Free melting: */ + xfree((void**)&melting); + + /*Get accumulation: */ + ModelFetchData((void**)&accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat"); + + count++; + param= new Param(count,"accumulation",DOUBLEVEC); + if(accumulation) param->SetDoubleVec(accumulation,model->numberofnodes); + else param->SetDoubleVec(accumulation,0); + parameters->AddObject(param); + + /*Free accumulation: */ + xfree((void**)&accumulation); + + + + + /*Assign output pointer: */ + *pparameters=parameters; +} Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 585) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 586) @@ -66,4 +66,5 @@ double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ double tria_viscosity_overshoot; + int tria_artdiff; /*penta constructor input: */ @@ -227,5 +228,5 @@ /*Create tria element using its constructor:*/ - tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff); /*Add tria element to elements dataset: */ @@ -479,10 +480,4 @@ xfree((void**)&model->uppernodes); - cleanup_and_return: - - /*Assign output pointer: */ - *pelements=elements; - *pnodes=nodes; - *pmaterials=materials; /*Keep partitioning information into model*/ @@ -497,3 +492,11 @@ VecFree(&gridborder); #endif + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + } Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 585) +++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 586) @@ -47,4 +47,9 @@ } + #ifdef _DEBUG_ + MatNorm(Kgg,NORM_INFINITY,&kmax2); + _printf_(" K_gg infinity norm after penalties: %g\n",kmax2); + #endif + /*Assign output pointers:*/ if(pkmax)*pkmax=kmax; Index: /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp =================================================================== --- /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 585) +++ /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 586) @@ -24,15 +24,7 @@ double* vy=NULL; double* vz=NULL; - + double* pressure=NULL; + double* u_g=NULL; - - /*thermal*/ - double* pressure=NULL; - double* temperature=NULL; - double* melting=NULL; - - double* p_g=NULL; - double* t_g=NULL; - double* m_g=NULL; /*control: */ @@ -46,13 +38,21 @@ double* u_g_obs=NULL; + double* p_g=NULL; + + /*prognostic: */ + double* a_g=NULL; + double* m_g=NULL; + double* h_g=NULL; + double* accumulation=NULL; + double* thickness=NULL; + double* melting=NULL; + /*diverse: */ int numberofnodes; int analysis_type; - int sub_analysis_type; int count; parameters->FindParam((void*)&analysis_type,"analysis_type"); - parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); count=parameters->Size(); @@ -63,5 +63,9 @@ - if ( (analysis_type==ControlAnalysisEnum()) || (analysis_type==DiagnosticAnalysisEnum()) || (analysis_type==ThermalAnalysisEnum())){ + if ( (analysis_type==ControlAnalysisEnum()) || + (analysis_type==DiagnosticAnalysisEnum()) || + (analysis_type==ThermalAnalysisEnum()) || + (analysis_type==PrognosticAnalysisEnum()) + ){ parameters->FindParam((void*)&vx,"vx"); @@ -170,39 +174,73 @@ param=(Param*)parameters->FindParamObject("pressure"); parameters->DeleteObject((Object*)param); - - if(sub_analysis_type==TransientAnalysisEnum()){ - - parameters->FindParam((void*)&temperature,"temperature"); - parameters->FindParam((void*)&melting,"melting"); - - /*Now, from temperature and melting, build t_g and m_g, correctly partitioned: */ - t_g=(double*)xcalloc(numberofnodes,sizeof(double)); - m_g=(double*)xcalloc(numberofnodes,sizeof(double)); - - for(i=0;iSetDoubleVec(t_g,numberofnodes); - parameters->AddObject(param); - - count++; - param= new Param(count,"m_g",DOUBLEVEC); - param->SetDoubleVec(m_g,numberofnodes); - parameters->AddObject(param); - - /*erase old parameter: */ - param=(Param*)parameters->FindParamObject("temperature"); - parameters->DeleteObject((Object*)param); - - param=(Param*)parameters->FindParamObject("melting"); - parameters->DeleteObject((Object*)param); - } - } - + } + + if(analysis_type==PrognosticAnalysisEnum()){ + + parameters->FindParam((void*)&melting,"melting"); + + /*Now, from melting, build m_g, correctly partitioned: */ + m_g=(double*)xcalloc(numberofnodes,sizeof(double)); + + for(i=0;iSetDoubleVec(m_g,numberofnodes); + parameters->AddObject(param); + + /*erase old parameter: */ + param=(Param*)parameters->FindParamObject("melting"); + parameters->DeleteObject((Object*)param); + + + + parameters->FindParam((void*)&thickness,"thickness"); + + /*Now, from thickness, build h_g, correctly partitioned: */ + h_g=(double*)xcalloc(numberofnodes,sizeof(double)); + + for(i=0;iSetDoubleVec(h_g,numberofnodes); + parameters->AddObject(param); + + /*erase old parameter: */ + param=(Param*)parameters->FindParamObject("thickness"); + parameters->DeleteObject((Object*)param); + + + + + parameters->FindParam((void*)&accumulation,"accumulation"); + + /*Now, from accumulation, build a_g, correctly partitioned: */ + a_g=(double*)xcalloc(numberofnodes,sizeof(double)); + + for(i=0;iSetDoubleVec(a_g,numberofnodes); + parameters->AddObject(param); + + /*erase old parameter: */ + param=(Param*)parameters->FindParamObject("accumulation"); + parameters->DeleteObject((Object*)param); + + + + } xfree((void**)&partition); @@ -215,10 +253,9 @@ xfree((void**)&vy_obs); xfree((void**)&pressure); - xfree((void**)&temperature); - xfree((void**)&melting); xfree((void**)&control_parameter); xfree((void**)&u_g_obs); xfree((void**)&p_g); - xfree((void**)&t_g); + xfree((void**)&a_g); + xfree((void**)&h_g); xfree((void**)&m_g); Index: /issm/trunk/src/c/ThicknessExtrudex/ThicknessExtrudex.cpp =================================================================== --- /issm/trunk/src/c/ThicknessExtrudex/ThicknessExtrudex.cpp (revision 586) +++ /issm/trunk/src/c/ThicknessExtrudex/ThicknessExtrudex.cpp (revision 586) @@ -0,0 +1,36 @@ +/*!\file ThicknessExtrudex + * \brief: vertical velocity extrusion + */ + +#include "./ThicknessExtrudex.h" + +#undef __FUNCT__ +#define __FUNCT__ "ThicknessExtrudex" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void ThicknessExtrudex( Vec tg, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials){ + + double* tg_serial=NULL; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes, materials); + nodes->Configure(elements,loads, nodes, materials); + + /*Serialize velcoity: */ + VecToMPISerial(&tg_serial,tg); + + /*Extrude velocity vertically: */ + elements->ThicknessExtrude(tg,tg_serial); + + /*Assemble vector: */ + VecAssemblyBegin(tg); + VecAssemblyEnd(tg); + + /*Free ressources:*/ + xfree((void**)&tg_serial); + +} Index: /issm/trunk/src/c/ThicknessExtrudex/ThicknessExtrudex.h =================================================================== --- /issm/trunk/src/c/ThicknessExtrudex/ThicknessExtrudex.h (revision 586) +++ /issm/trunk/src/c/ThicknessExtrudex/ThicknessExtrudex.h (revision 586) @@ -0,0 +1,14 @@ +/*!\file: ThicknessExtrudex.h + * \brief header file for velocity extrusion + */ + +#ifndef _THICKNESSEXTRUDEX_H +#define _THICKNESSEXTRUDEX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void ThicknessExtrudex( Vec tg, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials); + +#endif /* _THICKNESSEXTRUDEX_H */ + Index: /issm/trunk/src/c/VelocityDepthAveragex/VelocityDepthAveragex.cpp =================================================================== --- /issm/trunk/src/c/VelocityDepthAveragex/VelocityDepthAveragex.cpp (revision 586) +++ /issm/trunk/src/c/VelocityDepthAveragex/VelocityDepthAveragex.cpp (revision 586) @@ -0,0 +1,50 @@ +/*!\file VelocityDepthAveragex + * \brief: average velocity through thickness + */ + +#include "./VelocityDepthAveragex.h" + +#undef __FUNCT__ +#define __FUNCT__ "VelocityDepthAveragex" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void VelocityDepthAveragex( Vec ug, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials){ + + double* ug_serial=NULL; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes, materials); + nodes->Configure(elements,loads, nodes, materials); + + /*Serialize velcoity: */ + VecToMPISerial(&ug_serial,ug); + + /*Depth average velocity, onto base of mesh: */ + elements->VelocityDepthAverageAtBase(ug,ug_serial); + + /*Assemble vector: */ + VecAssemblyBegin(ug); + VecAssemblyEnd(ug); + + /*ok, ug is now correct at the base. We need to extrude the value of the base + * through thickness :*/ + + /*Serialize velcoity: */ + xfree((void**)&ug_serial); + VecToMPISerial(&ug_serial,ug); + + /*Extrude velocity, regardless of collapsed elements: */ + elements->VelocityExtrudeAllElements(ug,ug_serial); + + /*Assemble vector: */ + VecAssemblyBegin(ug); + VecAssemblyEnd(ug); + + /*Free ressources:*/ + xfree((void**)&ug_serial); + +} Index: /issm/trunk/src/c/VelocityDepthAveragex/VelocityDepthAveragex.h =================================================================== --- /issm/trunk/src/c/VelocityDepthAveragex/VelocityDepthAveragex.h (revision 586) +++ /issm/trunk/src/c/VelocityDepthAveragex/VelocityDepthAveragex.h (revision 586) @@ -0,0 +1,14 @@ +/*!\file: VelocityDepthAveragex.h + * \brief header file for velocity extrusion + */ + +#ifndef _VELOCITYDEPTHAVERAGEX_H +#define _VELOCITYDEPTHAVERAGEX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void VelocityDepthAveragex( Vec ug, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials); + +#endif /* _VELOCITYDEPTHAVERAGEX_H */ + Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 585) +++ /issm/trunk/src/c/issm.h (revision 586) @@ -45,7 +45,9 @@ #include "./ControlConstrainx/ControlConstrainx.h" #include "./VelocityExtrudex/VelocityExtrudex.h" +#include "./VelocityDepthAveragex/VelocityDepthAveragex.h" #include "./GriddataMeshToGridx/GriddataMeshToGridx.h" #include "./ComputePressurex/ComputePressurex.h" #include "./SlopeExtrudex/SlopeExtrudex.h" +#include "./ThicknessExtrudex/ThicknessExtrudex.h" Index: /issm/trunk/src/c/objects/DakotaPlugin.cpp =================================================================== --- /issm/trunk/src/c/objects/DakotaPlugin.cpp (revision 586) +++ /issm/trunk/src/c/objects/DakotaPlugin.cpp (revision 586) @@ -0,0 +1,92 @@ +/*!\file: DakotaPlugin.cpp (see DakotaPlugin.h for explanations) + */ + +#include "DakotaResponse.H" +#include "ParamResponsePair.H" +#include "DakotaPlugin.h" +#include "system_defs.h" +#include "ProblemDescDB.H" +#include "ParallelLibrary.H" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../objects/objects.h" +#include "../parallel/parallel.h" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include + + +namespace SIM { + +//constructor +DakotaPlugin::DakotaPlugin(const Dakota::ProblemDescDB& problem_db,FemModel* in_femmodels, ParameterInputs* in_inputs, int in_analysis_type, int in_sub_analysis_type):Dakota::DirectApplicInterface(problem_db){ + + femmodels=in_femmodels; + inputs=in_inputs; + analysis_type=in_analysis_type; + sub_analysis_type=in_sub_analysis_type; +} + +//destructor +DakotaPlugin::~DakotaPlugin(){ /* Virtual destructor handles referenceCount at Interface level. */ } + +int DakotaPlugin::derived_map_ac(const Dakota::String& driver) { + + int status=0; + int i; + double* variables=NULL; + char** variable_descriptors=NULL; + char* variable_descriptor=NULL; + double* responses=NULL; + + /*Before launching analysis, we need to transfer the dakota inputs into Issm + *readable variables: */ + + /*First, the variables: */ + variables=(double*)xmalloc(numACV*sizeof(double)); + for(i=0;iGetViscosity3d(&viscosity, &epsilon[0]); @@ -525,8 +525,8 @@ /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4); - + /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant onto this scalar matrix, so that we win some computational time: */ - + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=newviscosity*gauss_weight*Jdet; @@ -534,10 +534,10 @@ D[i][i]=D_scalar; } - + /* Do the triple product tB*D*Bprime: */ TripleMultiply( &B[0][0],5,numdof,1, - &D[0][0],5,5,0, - &Bprime[0][0],5,numdof,0, - &Ke_gg_gaussian[0][0],0); + &D[0][0],5,5,0, + &Bprime[0][0],5,numdof,0, + &Ke_gg_gaussian[0][0],0); /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ @@ -549,9 +549,9 @@ } //for (ig2=0; ig2GetRhoWater(); @@ -818,8 +818,8 @@ /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore - get tria gaussian points as well as segment gaussian points. For tria gaussian - points, order of integration is 2, because we need to integrate the product tB*D*B' - which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian - points, same deal, which yields 3 gaussian points.*/ + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ area_order=5; @@ -845,5 +845,5 @@ /*Compute strain rate: */ GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); - + /*Get viscosity: */ matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); @@ -883,9 +883,9 @@ /*Build alpha2_list used by drag stiffness matrix*/ Friction* friction=NewFriction(); - + /*Initialize all fields: */ friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); strcpy(friction->element_type,"2d"); - + friction->gravity=matpar->GetG(); friction->rho_ice=matpar->GetRhoIce(); @@ -910,5 +910,5 @@ } } - + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2); @@ -920,25 +920,25 @@ gauss_coord[2]=*(third_gauss_area_coord+igarea); gauss_coord[3]=-1; - + gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); gauss_coord_tria[1]=*(second_gauss_area_coord+igarea); gauss_coord_tria[2]=*(third_gauss_area_coord+igarea); - + /*Get the Jacobian determinant */ tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); - + /*Get L matrix if viscous basal drag present: */ GetLStokes(&LStokes[0][0], gauss_coord_tria); GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord); - + /*Compute strain rate: */ GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); - + /*Get viscosity at last iteration: */ matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - + /*Get normal vecyor to the bed */ SurfaceNormal(&surface_normal[0],xyz_list_tria); - + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result bed_normal[1]=-surface_normal[1]; @@ -974,5 +974,5 @@ } } //if ( (onbed==1) && (shelf==0)) - + /*Reduce the matrix */ ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]); @@ -986,6 +986,6 @@ /*Add Ke_gg to global matrix Kgg: */ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES); - -cleanup_and_return: + + cleanup_and_return: return; @@ -998,9 +998,9 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ if (analysis_type==ControlAnalysisEnum()){ - + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==DiagnosticAnalysisEnum()){ - + if (sub_analysis_type==HorizAnalysisEnum()){ @@ -1008,9 +1008,9 @@ } else if (sub_analysis_type==VertAnalysisEnum()){ - + CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type); } else if (sub_analysis_type==StokesAnalysisEnum()){ - + CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type); } @@ -1018,13 +1018,13 @@ } else if (analysis_type==SlopeComputeAnalysisEnum()){ - + CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==ThermalAnalysisEnum()){ - + CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==MeltingAnalysisEnum()){ - + CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type); } @@ -1056,5 +1056,5 @@ inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes); inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes); - + //Update material if necessary if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){ @@ -1063,5 +1063,5 @@ matice->SetB(B_average); } - + if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){ B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0; @@ -1088,5 +1088,5 @@ } } - + int Penta::GetOnBed(){ return onbed; @@ -1099,5 +1099,5 @@ } void Penta::GetBedList(double* bed_list){ - + int i; for(i=0;i<6;i++)bed_list[i]=b[i]; @@ -1115,5 +1115,5 @@ int i; Tria* tria=NULL; - + /*Bail out if this element if: * -> Non collapsed and not on the surface @@ -1132,5 +1132,5 @@ } else{ - + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). tria->Du(du_g,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type); @@ -1156,8 +1156,8 @@ #define __FUNCT__ "Penta::GradjDrag" void Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){ - - + + Tria* tria=NULL; - + /*Bail out if this element does not touch the bedrock: */ if (!onbed){ @@ -1165,5 +1165,5 @@ } else{ - + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type); @@ -1178,12 +1178,12 @@ throw ErrorException(__FUNCT__," not supported yet!"); } - + #undef __FUNCT__ #define __FUNCT__ "Penta::Misfit" double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type,int sub_analysis_type){ - + double J; Tria* tria=NULL; - + /*Bail out if this element if: * -> Non collapsed and not on the surface @@ -1209,5 +1209,5 @@ } } - + #undef __FUNCT__ #define __FUNCT__ "Penta::SpawnTria" @@ -1224,5 +1224,5 @@ double tria_accumulation[3]; double tria_geothermalflux[3]; - + /*configuration: */ int tria_node_ids[3]; @@ -1249,5 +1249,5 @@ tria_melting[1]=melting[g1]; tria_melting[2]=melting[g2]; - + tria_accumulation[0]=accumulation[g0]; tria_accumulation[1]=accumulation[g1]; @@ -1270,5 +1270,5 @@ tria_node_offsets[2]=node_offsets[g2]; - tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot); + tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot,artdiff); tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets); @@ -1286,5 +1286,5 @@ int doflist_per_node[MAXDOFSPERNODE]; int numberofdofspernode; - + for(i=0;i<6;i++){ nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); @@ -1320,5 +1320,5 @@ GetB(&B[0][0], xyz_list, gauss_l1l2l3l4); -#ifdef _DEBUG_ + #ifdef _DEBUG_ _printf_("B for grid1 : [ %lf %lf \n",B[0][0],B[0][1]); _printf_(" [ %lf %lf \n",B[1][0],B[1][1]); @@ -1326,5 +1326,5 @@ _printf_(" [ %lf %lf ]\n",B[3][0],B[3][1]); _printf_(" [ %lf %lf ]\n",B[4][0],B[4][1]); - + _printf_("B for grid2 : [ %lf %lf \n",B[0][2],B[0][3]); _printf_(" [ %lf %lf \n",B[1][2],B[1][3]); @@ -1332,5 +1332,5 @@ _printf_(" [ %lf %lf ]\n",B[3][2],B[3][3]); _printf_(" [ %lf %lf ]\n",B[4][2],B[4][3]); - + _printf_("B for grid3 : [ %lf %lf \n", B[0][4],B[0][5]); _printf_(" [ %lf %lf \n", B[1][4],B[1][5]); @@ -1338,5 +1338,5 @@ _printf_(" [ %lf %lf ]\n",B[3][4],B[3][5]); _printf_(" [ %lf %lf ]\n",B[4][4],B[4][5]); - + _printf_("B for grid4 : [ %lf %lf \n", B[0][6],B[0][7]); _printf_(" [ %lf %lf \n", B[1][6],B[1][7]); @@ -1344,5 +1344,5 @@ _printf_(" [ %lf %lf ]\n",B[3][6],B[3][7]); _printf_(" [ %lf %lf ]\n",B[4][6],B[4][7]); - + _printf_("B for grid5 : [ %lf %lf \n", B[0][8],B[0][9]); _printf_(" [ %lf %lf \n", B[1][8],B[1][9]); @@ -1360,10 +1360,10 @@ _printf_("Velocity for grid %i %lf %lf\n",i,*(vxvy_list+2*i+0),*(vxvy_list+2*i+1)); } -#endif + #endif /*Multiply B by velocity, to get strain rate: */ MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, - velocity,NDOF2*numgrids,1,0, - epsilon,0); + velocity,NDOF2*numgrids,1,0, + epsilon,0); } @@ -1385,5 +1385,5 @@ * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) */ - + int i; const int numgrids=6; @@ -1396,9 +1396,9 @@ GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4); -#ifdef _DEBUG_ + #ifdef _DEBUG_ for (i=0;iGetDofList(&nodedofs[0],&numberofdofspernode); + VecSetValues(ug,1,&nodedofs[0],&velocity[0],INSERT_VALUES); + VecSetValues(ug,1,&nodedofs[1],&velocity[1],INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + + } + +} + +#undef __FUNCT__ +#define __FUNCT__ "Penta::ThicknessExtrude" +void Penta::ThicknessExtrude(Vec tg,double* tg_serial){ + + /* node data: */ + const int numgrids=6; + const int numdof=1*numgrids; + int doflist[numdof]; + int nodedofs; + int numberofdofspernode; + + Node* node=NULL; + int i; + double thickness; + + + if(onbed==1){ + + GetDofList(&doflist[0],&numberofdofspernode); + + /*this penta is on the bed. For each node on the base of this penta, + * we grab the thickness. Once we know the thickness, we follow the upper nodes, + * inserting the same thickness value into tg, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + + /*get velocity for this base node: */ + thickness=tg_serial[doflist[numberofdofspernode*i+0]]; + + //go through all nodes which sit on top of this node, until we reach the surface, + //and pltg velocity in tg + for(;;){ + + node->GetDofList(&nodedofs,&numberofdofspernode); + VecSetValues(tg,1,&nodedofs,&thickness,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + + } + +} + +#undef __FUNCT__ +#define __FUNCT__ "Penta::VelocityDepthAverageAtBase" +void Penta::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){ + + /* node data: */ + const int numgrids=6; + const int numdof=2*numgrids; + int doflist[numdof]; + int nodedofs[2]; + int numberofdofspernode; + int dofx,dofy; + + Node* node=NULL; + Node* upper_node=NULL; + int i; + double base_velocity[2]; + double velocity2[2]; + double velocity1[2]; + double velocity_average[2]; + double sum[2]; + double z1,z2,dz; + double thickness; + + if(onbed==1){ + + GetDofList(&doflist[0],&numberofdofspernode); + + /*this penta is on the bed. For each node on the base of this penta, we are going to + * follow the upper nodes until we reach the surface. At each upper node, we'll grab the + * velocity for this node, and add it to ug. We'll finish by + * we grab the velocity. Once we know the velocity, we follow the upper nodes, + * inserting the same velocity value into ug, until we reach the surface: */ + for(i=0;i<3;i++){ + + //get thickness for this grid + thickness=h[i]; //pick up from the penta element itself. + + node=nodes[i]; //base nodes + + /*get dofs for this base node velocity: */ + dofx=doflist[numberofdofspernode*i+0]; + dofy=doflist[numberofdofspernode*i+1]; + + /*first thing, cancel velocity on the base nodes, so that we start with the value 0: */ + base_velocity[0]=-ug_serial[doflist[numberofdofspernode*i+0]]; //- to cancel + base_velocity[1]=-ug_serial[doflist[numberofdofspernode*i+1]]; + + VecSetValues(ug,1,&dofx,&base_velocity[0],ADD_VALUES); + VecSetValues(ug,1,&dofy,&base_velocity[1],ADD_VALUES); + + //go through all nodes which sit on top of this node, until we reach the surface, + //and plug deltaV*deltaH/H at base of ug: */ + + for(;;){ + + if (node->IsOnSurface())break; + + node->GetDofList(&nodedofs[0],&numberofdofspernode); + + velocity1[0]=ug_serial[nodedofs[numberofdofspernode*i+0]]; + velocity1[1]=ug_serial[nodedofs[numberofdofspernode*i+1]]; + z1=node->GetZ(); + + upper_node=node->GetUpperNode(); + upper_node->GetDofList(&nodedofs[0],&numberofdofspernode); + + velocity2[0]=ug_serial[nodedofs[numberofdofspernode*i+0]]; + velocity2[1]=ug_serial[nodedofs[numberofdofspernode*i+1]]; + z2=upper_node->GetZ(); + + dz=(z2-z1); + velocity_average[0]=(velocity1[0]+velocity2[0])/2.0; + velocity_average[1]=(velocity1[1]+velocity2[1])/2.0; + + sum[0]=velocity_average[0]*dz/thickness; + sum[1]=velocity_average[1]*dz/thickness; + + /*Plug velocity_average*deltaH/H into base of ug: */ + VecSetValues(ug,1,&dofx,&sum[0],ADD_VALUES); + VecSetValues(ug,1,&dofy,&sum[1],ADD_VALUES); + + /*get next node: */ + node=node->GetUpperNode(); + } + } + + } + +} + #undef __FUNCT__ @@ -1938,12 +2125,12 @@ int nodedof; int numberofdofspernode; - + Node* node=NULL; int i; double slope; - + if(onbed==1){ - + GetDofList(&doflist[0],&numberofdofspernode); @@ -1952,9 +2139,9 @@ /*For each node on the base of this penta, we grab the slope. Once we know the slope, we follow the upper nodes, * inserting the same slope value into sg, until we reach the surface: */ - + for(i=0;i<3;i++){ - + node=nodes[i]; //base nodes - + /*get velocity for this base node: */ slope=sg_serial[doflist[i]]; @@ -1984,5 +2171,5 @@ /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; - where hi is the interpolation function for grid i.*/ + where hi is the interpolation function for grid i.*/ int i; @@ -1994,9 +2181,9 @@ GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4); -#ifdef _DEBUG_ + #ifdef _DEBUG_ for (i=0;iRecover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!"); - inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); + inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); /*Get gaussian points and weights :*/ @@ -2101,5 +2288,5 @@ GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); -#ifdef _DEBUG_ + #ifdef _DEBUG_ for (i=0;iGetViscosity3dStokes(&viscosity,&epsilon[0]); - + /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); @@ -2984,5 +3171,5 @@ GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); - + /*Get bubble part of Bprime */ for(i=0;i<8;i++){ @@ -3033,5 +3220,5 @@ gauss_coord[2]=*(third_gauss_area_coord+igarea); gauss_coord[3]=-1; - + gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); gauss_coord_tria[1]=*(second_gauss_area_coord+igarea); @@ -3040,5 +3227,5 @@ /*Get the Jacobian determinant */ tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); - + /* Get bed at gaussian point */ GetParameterValue(&bed,&b[0],gauss_coord); @@ -3046,5 +3233,5 @@ /*Get L matrix : */ tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1); - + /*Get water_pressure at gaussian point */ water_pressure=gravity*rho_water*bed; @@ -3052,5 +3239,5 @@ /*Get normal vecyor to the bed */ SurfaceNormal(&surface_normal[0],xyz_list_tria); - + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result bed_normal[1]=-surface_normal[1]; @@ -3064,5 +3251,5 @@ } } //if ( (onbed==1) && (shelf==1)) - + /*Reduce the matrix */ ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]); @@ -3080,5 +3267,5 @@ #define __FUNCT__ "Penta::ReduceVectorStokes" void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){ - + int i,j; @@ -3123,5 +3310,5 @@ #define __FUNCT__ "Penta::GetNodalFunctionsStokes" void Penta::GetNodalFunctionsStokes(double* l1l7, double* gauss_coord){ - + /*This routine returns the values of the nodal functions at the gaussian point.*/ @@ -3181,5 +3368,4 @@ int dofs[3]={0,1,2}; double dt; - int dt_exists; double vxvyvz_list[numgrids][3]; @@ -3195,29 +3381,29 @@ /*matrices: */ - double K_terms[numdof][numdof]={0.0}; - double Ke_gaussian_conduct[numdof][numdof]; - double Ke_gaussian_advec[numdof][numdof]; - double Ke_gaussian_transient[numdof][numdof]; - double B[3][numdof]; - double Bprime[3][numdof]; - double B_conduct[3][numdof]; - double B_advec[3][numdof]; - double Bprime_advec[3][numdof]; - double L[numdof]; - double D_scalar; - double D[3][3]; - double l1l2l3[3]; - double tl1l2l3D[3]; - double tBD[3][numdof]; - double tBD_conduct[3][numdof]; - double tBD_advec[3][numdof]; - double tLD[numdof]; - - double Jdet; + double K_terms[numdof][numdof]={0.0}; + double Ke_gaussian_conduct[numdof][numdof]; + double Ke_gaussian_advec[numdof][numdof]; + double Ke_gaussian_transient[numdof][numdof]; + double B[3][numdof]; + double Bprime[3][numdof]; + double B_conduct[3][numdof]; + double B_advec[3][numdof]; + double Bprime_advec[3][numdof]; + double L[numdof]; + double D_scalar; + double D[3][3]; + double l1l2l3[3]; + double tl1l2l3D[3]; + double tBD[3][numdof]; + double tBD_conduct[3][numdof]; + double tBD_advec[3][numdof]; + double tLD[numdof]; + + double Jdet; /*Material properties: */ - double gravity,rho_ice,rho_water; - double heatcapacity,thermalconductivity; - double mixed_layer_capacity,thermal_exchange_velocity; + double gravity,rho_ice,rho_water; + double heatcapacity,thermalconductivity; + double mixed_layer_capacity,thermal_exchange_velocity; /*Collapsed formulation: */ @@ -3232,5 +3418,5 @@ GetDofList(&doflist[0],&numberofdofspernode); - /*recovre material parameters: */ + // /*recovre material parameters: */ rho_water=matpar->GetRhoWater(); rho_ice=matpar->GetRhoIce(); @@ -3238,4 +3424,5 @@ heatcapacity=matpar->GetHeatCapacity(); thermalconductivity=matpar->GetThermalConductivity(); + /*recover extra inputs from users, dt and velocity: */ @@ -3294,4 +3481,5 @@ D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar; + /* Do the triple product B'*D*B: */ MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0); @@ -3303,4 +3491,5 @@ GetB_advec(&B_advec[0][0],&xyz_list[0][0],gauss_coord); GetBprime_advec(&Bprime_advec[0][0],&xyz_list[0][0],gauss_coord); + //Build the D matrix @@ -3322,4 +3511,5 @@ MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0); MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0); + /*Transient: */ @@ -3361,5 +3551,5 @@ xfree((void**)&vert_gauss_weights); xfree((void**)&vert_gauss_coord); - + //Ice/ocean heat exchange flux on ice shelf base if(onbed && shelf){ @@ -3385,5 +3575,5 @@ * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) */ - + int i; int calculationdof=3; @@ -3419,5 +3609,5 @@ * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) */ - + int i; int calculationdof=3; @@ -3454,5 +3644,5 @@ * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) */ - + int i; int calculationdof=3; @@ -3477,5 +3667,5 @@ #define __FUNCT__ "Penta::CreateKMatrixMelting" void Penta::CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ - + Tria* tria=NULL; if (!onbed){ @@ -3483,5 +3673,5 @@ } else{ - + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type); @@ -3508,5 +3698,5 @@ /*Grid data: */ double xyz_list[numgrids][3]; - + /* gaussian points: */ int num_area_gauss,igarea,igvert; @@ -3521,5 +3711,5 @@ int area_order=2; int num_vert_gauss=3; - + double dt; double vx_list[numgrids]; @@ -3537,5 +3727,5 @@ /* element parameters: */ int friction_type; - + int dofs[3]={0,1,2}; int dofs1[1]={0}; @@ -3572,5 +3762,5 @@ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - + // /*recovre material parameters: */ rho_water=matpar->GetRhoWater(); @@ -3592,5 +3782,5 @@ if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); } - + for(i=0;iGetViscosity3dStokes(&viscosity,&epsilon[0]); + matice->GetViscosity3d(&viscosity,&epsilon[0]); /* Get Jacobian determinant: */ @@ -3661,4 +3851,5 @@ /* Ice/ocean heat exchange flux on ice shelf base */ if(onbed && shelf){ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 585) +++ /issm/trunk/src/c/objects/Penta.h (revision 586) @@ -110,5 +110,8 @@ void GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4); void GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4); + void VelocityExtrudeAllElements(Vec ug,double* ug_serial); void VelocityExtrude(Vec ug,double* ug_serial); + void ThicknessExtrude(Vec ug,double* ug_serial); + void VelocityDepthAverageAtBase(Vec ug,double* ug_serial); void SlopeExtrude(Vec sg,double* sg_serial); void ComputePressure(Vec p_g); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 585) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 586) @@ -17,4 +17,5 @@ #include "../include/typedefs.h" + /*For debugging purposes: */ #define ELID 141 //element for which to print debug statements @@ -33,5 +34,5 @@ Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3], double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel, - double tria_viscosity_overshoot){ + double tria_viscosity_overshoot,int tria_artdiff){ int i; @@ -50,6 +51,5 @@ melting[i]=tria_melting[i]; accumulation[i]=tria_accumulation[i]; - geothermalflux[i]=tria_geothermalflux[i]; - } + geothermalflux[i]=tria_geothermalflux[i]; } matice=NULL; matice_offset=UNDEF; @@ -64,4 +64,5 @@ onbed=1; viscosity_overshoot=tria_viscosity_overshoot; + artdiff=tria_artdiff; return; @@ -96,4 +97,5 @@ printf(" onbed: %i\n",onbed); printf(" viscosity_overshoot=%g\n",viscosity_overshoot); + printf(" artdiff=%g\n",artdiff); printf(" nodes: \n"); if(nodes[0])nodes[0]->Echo(); @@ -146,4 +148,5 @@ memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel); memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot); + memcpy(marshalled_dataset,&artdiff,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff); *pmarshalled_dataset=marshalled_dataset; @@ -177,4 +180,5 @@ +sizeof(epsvel) +sizeof(viscosity_overshoot) + +sizeof(artdiff) +sizeof(int); //sizeof(int) for enum type } @@ -220,4 +224,5 @@ memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel); memcpy(&viscosity_overshoot,marshalled_dataset,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot); + memcpy(&artdiff,marshalled_dataset,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff); /*nodes and materials are not pointing to correct objects anymore:*/ @@ -242,4 +247,5 @@ } + #undef __FUNCT__ #define __FUNCT__ "Tria::Configure" @@ -290,8 +296,15 @@ } + else if (analysis_type==PrognosticAnalysisEnum()){ + + CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type); + + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); } -} + +} + #undef __FUNCT__ @@ -299,4 +312,5 @@ void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ + /* local declarations */ @@ -509,4 +523,208 @@ } +#undef __FUNCT__ +#define __FUNCT__ "Tria::CreateKMatrixPrognostic" +void Tria::CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ + + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + /* matrices: */ + double L[numgrids]; + double DL[2][2]={0.0}; + double DLprime[2][2]={0.0}; + double DL_scalar; + double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_thickness[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + + double Jdettria; + + /*input parameters for structural analysis (diagnostic): */ + double vxvy_list[numgrids][2]={0,0}; + double vx_list[numgrids]={0,0}; + double vy_list[numgrids]={0,0}; + double dvx[2]; + double dvy[2]; + double vx,vy; + double v_gauss[2]={0.0}; + double K[2][2]={0.0}; + double dt; + int dofs[2]={0,1}; + int found=0; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*recover extra inputs from users, at current convergence iteration: */ + found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity_average in inputs!"); + + for(i=0;iRecover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + //Create Artificial diffusivity once for all if requested + if(artdiff){ + //Get the Jacobian determinant + gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0; + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + + //Build K matrix (artificial diffusivity matrix) + v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); + v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); + + K[0][0]=pow(Jdettria,.5)/2.0*fabs(v_gauss[0]); + K[1][1]=pow(Jdettria,.5)/2.0*fabs(v_gauss[1]); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igpartition; + + /*serialize outputs: */ + VecDuplicate(partition,&partition_shifted); + VecCopy(partition,partition_shifted); + VecShift(partition_shifted,1.0); //matlab indexing starts at 1 + + VecToMPISerial(&serial_partition,partition_shifted); + VecGetSize(partition,&nods); + + VecToMPISerial(&serial_hg,h_g); + VecGetSize(h_g,&h_g_size); + + /* Open output file to write raw binary data: */ + if(my_rank==0){ + fid=pfopen(filename,"wb"); + + /*Write solution type: */ + WriteDataToDisk(analysis_type,NULL,NULL,"String",fid); + + /*Write partition: */ + WriteDataToDisk(serial_partition,&nods,&one,"Mat",fid); + + /*Write solution to disk: */ + WriteDataToDisk(serial_hg,&h_g_size,&one,"Mat",fid); + + /*Close file: */ + pfclose(fid,filename); + } + + /*Free ressources: */ + VecFree(&partition_shifted); + xfree((void**)&serial_partition); + xfree((void**)&serial_hg); +} Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 585) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 586) @@ -23,5 +23,7 @@ char* outputfilename=NULL; char* lockname=NULL; + char* qmuname=NULL; int numberofnodes; + int qmu_analysis=0; /*Fem models : */ @@ -52,4 +54,5 @@ outputfilename=argv[3]; lockname=argv[4]; + qmuname=argv[5]; /*Open handle to data on disk: */ @@ -68,4 +71,5 @@ CreateFemModel(&femmodels[4],fid,"slope_compute",""); + _printf_("initialize inputs:\n"); femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g"); @@ -79,9 +83,21 @@ femmodels[0].parameters->DeleteObject((Object*)param); - _printf_("call computational core:\n"); - diagnostic_core(&u_g,&p_g,femmodels,inputs); + /*are we running the solutoin sequence, or a qmu wrapper around it? : */ + femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis"); + if(!qmu_analysis){ + /*run diagnostic analysis: */ + _printf_("call computational core:\n"); + diagnostic_core(&u_g,&p_g,femmodels,inputs); - _printf_("write results to disk:\n"); - OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename); + _printf_("write results to disk:\n"); + OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename); + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on diagnostic core:\n"); + + qmu(qmuname,&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum()); + } _printf_("write lock file:\n"); Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 585) +++ /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 586) @@ -14,4 +14,6 @@ void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs){ + + extern int my_rank; /*fem models: */ Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 585) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 586) @@ -46,4 +46,6 @@ kflag=1; pflag=1; + + fem->parameters->FindParam((void*)&connectivity,"connectivity"); fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode"); @@ -57,5 +59,5 @@ /*Copy loads for backup: */ loads=fem->loads->Copy(); - + /*Initialize ug and ug_old */ if (numberofdofspernode>=3)dofs[2]=1;//only keep vz if running with more than 3 dofs per node Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 585) +++ /issm/trunk/src/c/parallel/parallel.h (revision 586) @@ -31,8 +31,12 @@ void OutputThermal(Vec* t_g,Vec* m_g, double* time,FemModel* femmodels,char* filename); void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets); +void OutputPrognostic(Vec h_g,FemModel* femmodel,char* filename); void WriteLockFile(char* filename); void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type); //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename); +void qmu(const char* dakota_input_file,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void SpawnCore(double* responses,double* variables,char** variable_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Vec u_g,Vec p_g,int analysis_type,int sub_analysis_type); #endif Index: /issm/trunk/src/c/parallel/prognostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic.cpp (revision 586) +++ /issm/trunk/src/c/parallel/prognostic.cpp (revision 586) @@ -0,0 +1,109 @@ +/*!\file: prognostic.cpp + * \brief: prognostic solution + */ + +#include "../issm.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "prognostic" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + int waitonlock=0; + + FemModel fem; + Vec h_g=NULL; + Vec u_g=NULL; + double* u_g_serial=NULL; + double* h_g_initial=NULL; + double* melting_g=NULL; + double* accumulation_g=NULL; + double dt; + + + ParameterInputs* inputs=NULL; + Param* param=NULL; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!"); + #endif + + PetscInitialize(&argc,&argv,(char *)0,""); + + /*Size and rank: */ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); + + _printf_("recover , input file name and output file name:\n"); + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + _printf_("read and create finite element model:\n"); + CreateFemModel(&fem,fid,"prognostic",""); + + //retrieve parameters used to fill inputs + fem.parameters->FindParam((void*)&u_g_serial,"u_g"); + fem.parameters->FindParam((void*)&h_g_initial,"h_g"); + fem.parameters->FindParam((void*)&melting_g,"m_g"); + fem.parameters->FindParam((void*)&accumulation_g,"a_g"); + fem.parameters->FindParam((void*)&dt,"dt"); + fem.parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + + _printf_("depth averaging velocity..."); + u_g=SerialToVec(u_g_serial,numberofnodes*3); xfree((void**)&u_g_serial);//vx,vy and vz should be present at this point. + VelocityDepthAveragex( u_g, fem.elements,fem.nodes, fem.loads, fem.materials); + + _printf_("initialize inputs:\n"); + inputs=new ParameterInputs; + inputs->Add("velocity_average",u_g,3,numberofnodes); + inputs->Add("thickness",h_g_initial,1,numberofnodes); + inputs->Add("melting",melting_g,1,numberofnodes); + inputs->Add("accumulation",accumulation_g,1,numberofnodes); + inputs->Add("dt",dt); + + /*lighten up on parameters : to be done */ + + _printf_("call computational core:\n"); + diagnostic_core_linear(&h_g,&fem,inputs,PrognosticAnalysisEnum(),NoneAnalysisEnum()); + + _printf_("write results to disk:\n"); + OutputPrognostic(h_g,&fem,outputfilename); + + _printf_("write lock file:\n"); + fem.parameters->FindParam((void*)&waitonlock,"waitonlock"); + if (waitonlock){ + WriteLockFile(lockname); + } + + _printf_("closing MPI and Petsc\n"); + MPI_Barrier(MPI_COMM_WORLD); + + /*Close MPI libraries: */ + PetscFinalize(); + + + /*end module: */ + MODULEEND(); + + return 0; //unix success return; +} Index: /issm/trunk/src/c/parallel/qmu.cpp =================================================================== --- /issm/trunk/src/c/parallel/qmu.cpp (revision 586) +++ /issm/trunk/src/c/parallel/qmu.cpp (revision 586) @@ -0,0 +1,94 @@ +/*!\file: qmu.cpp + * \brief: run qmu capabilities using Dakota library + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "qmu" + +#include "ParallelLibrary.H" +#include "ProblemDescDB.H" +#include "DakotaStrategy.H" +#include "DakotaModel.H" +#include "DakotaInterface.H" + +#include "../objects/DakotaPlugin.h" +#include "../shared/shared.h" +#include "../include/macros.h" +#include "parallel.h" + +void qmu(const char* dakota_input_file,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + + extern int my_rank; + int status=0; + Dakota::ModelLIter ml_iter; + + if(my_rank==0){ + + // Instantiate/initialize the parallel library and problem description + // database objects. + Dakota::ParallelLibrary parallel_lib("serial"); + 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(NULL, NULL, 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,femmodels,inputs,analysis_type,sub_analysis_type), false); + } + + // Execute the strategy + problem_db.lock(); // prevent run-time DB queries + selected_strategy.run_strategy(); + + //Warn other cpus that we are done running the dakota iterator + MPI_Bcast(&status,1,MPI_INT,0,MPI_COMM_WORLD); + + } + else{ + + for(;;){ + + SpawnCore(NULL,NULL,NULL,0,femmodels,inputs,analysis_type,sub_analysis_type); + + //Figure out if cpu 0 is done iterating + MPI_Bcast(&status,1,MPI_INT,0,MPI_COMM_WORLD); + + //yes! + if(status){ + break; //yes, we are done + } + else{ + /*We are not done: */ + + + } + } + + + } +} Index: /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 585) +++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 586) @@ -42,4 +42,7 @@ numdofs=1; } + else if (analysis_type==PrognosticAnalysisEnum()){ + numdofs=1; + } else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis type: ",analysis_type," not implemented yet!")); Index: /issm/trunk/src/c/shared/String/root.cpp =================================================================== --- /issm/trunk/src/c/shared/String/root.cpp (revision 586) +++ /issm/trunk/src/c/shared/String/root.cpp (revision 586) @@ -0,0 +1,35 @@ +/*!\file: root.cpp + * \brief get root of a string: ex, "thickness"=root("thickness1"); + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "root" + +char* root(char* string){ + + int outlength; + char* outstring=NULL; + int length; + + length=strlen(string); + + outlength=length; + for(i=0;i +extern "C" void METIS_PartMeshNodal(int *, int *, idxtype *, int *, int *, int *, int *, idxtype *, idxtype *); + #endif Index: /issm/trunk/src/c/toolkits/petsc/patches/SerialToVec.cpp =================================================================== --- /issm/trunk/src/c/toolkits/petsc/patches/SerialToVec.cpp (revision 586) +++ /issm/trunk/src/c/toolkits/petsc/patches/SerialToVec.cpp (revision 586) @@ -0,0 +1,63 @@ +/* \file SerialToVec.cpp + * \brief: convert a serial vector on all cpus into a parallel vector + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +/*Petsc includes: */ +#include "petscmat.h" +#include "petscvec.h" +#include "petscksp.h" + +/*Matlab includes: */ +#include "mex.h" + +#include "../../../shared/shared.h" + +Vec SerialToVec(double* vector,int vector_size){ + + int i; + + /*output: */ + Vec outvector=NULL; + + /*petsc indices: */ + int* idxn=NULL; + double* values=NULL; + int lower_row,upper_row,range; + + + /*Create parallel vector: */ + outvector=NewVec(vector_size); + + /*plug values from local vector into new parallel vector: */ + VecGetOwnershipRange(outvector,&lower_row,&upper_row); + upper_row--; + range=upper_row-lower_row+1; + + if (range){ + idxn=(int*)xmalloc(range*sizeof(int)); + values=(double*)xmalloc(range*sizeof(double)); + for (i=0;i %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name); +fprintf(fid,' %s %s.bin %s.outbin %s.lock qmu.in 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name); fclose(fid); Index: /issm/trunk/src/m/classes/public/LaunchQueueJobGeneric.m =================================================================== --- /issm/trunk/src/m/classes/public/LaunchQueueJobGeneric.m (revision 585) +++ /issm/trunk/src/m/classes/public/LaunchQueueJobGeneric.m (revision 586) @@ -20,7 +20,15 @@ disp('uploading input file and queueing script'); if strcmpi(hostname,md.cluster), - system(['cp ' md.name '.bin' ' ' md.name '.queue' ' ' executionpath]); + if md.qmu_analysis, + system(['cp ' md.name '.bin' ' ' md.name '.queue qmu/qmu.in ' ' ' executionpath]); + else + system(['cp ' md.name '.bin' ' ' md.name '.queue' ' ' executionpath]); + end else - system(['scp ' md.name '.bin' ' ' md.name '.queue' ' ' md.cluster ':' executionpath]); + if md.qmu_analysis, + system(['scp ' md.name '.bin' ' ' md.name '.queue qmu/qmu.in ' ' ' md.cluster ':' executionpath]); + else + system(['scp ' md.name '.bin' ' ' md.name '.queue' ' ' md.cluster ':' executionpath]); + end end Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 585) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 586) @@ -153,4 +153,14 @@ WriteData(fid,md.rifts(i).friction,'Scalar',['friction' num2str(i)]); end + +%Qmu +WriteData(fid,md.qmu_analysis,'Integer','qmu_analysis'); +if md.qmu_analysis, + responses=md.responses(md.iresp).rf; + WriteData(fid,numel(responses),'Integer','numberofresponses'); + for i=1:numel(responses), + WriteData(fid,responses(i).descriptor,'String',['descriptor' num2str(i-1)]); + end +end %Get penalties to connect collapsed and non-collapsed 3d meshes: Index: /issm/trunk/src/m/classes/public/process_solve_options.m =================================================================== --- /issm/trunk/src/m/classes/public/process_solve_options.m (revision 585) +++ /issm/trunk/src/m/classes/public/process_solve_options.m (revision 586) @@ -46,5 +46,5 @@ end if ~found, - sub_analysis_type=''; + sub_analysis_type='none'; end @@ -56,9 +56,8 @@ strcmpi(analysis_type,'mesh') | ... strcmpi(analysis_type,'mesh2grid') | ... - strcmpi(analysis_type,'qmu') | ... strcmpi(analysis_type,'transient') ), error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']); end -if ~(strcmpi(sub_analysis_type,'') | ... +if ~(strcmpi(sub_analysis_type,'none') | ... strcmpi(sub_analysis_type,'steady') | ... strcmpi(sub_analysis_type,'horiz') | ... Index: /issm/trunk/src/m/classes/public/qmu.m =================================================================== --- /issm/trunk/src/m/classes/public/qmu.m (revision 586) +++ /issm/trunk/src/m/classes/public/qmu.m (revision 586) @@ -0,0 +1,65 @@ +function md=qmu(md,varargin) +%QMU - apply Quantification of Margins and Uncertainties techniques +% to a solution sequence (like diagnostic.m, progonstic.m, etc ...), +% using the Dakota software from Sandia. +% Usage: +% md=qmu(md,varargin) +% where varargin is a lit of paired arguments. +% arguments can be: 'analysis_type': 'diagnostic','thermal','prognostic','transient' +% arguments can be: 'sub_analysis_type': 'transient','steady','' (default if empty = 'steady') +% arguments can be: 'package': 'macayeal','ice','cielo' (default if not specified = 'cielo') +% +% +% Examples: +% md=qmu(md,'analysis_type','diagnostic','package','cielo'); +% md=qmu(md,'analysis_type','control','package','ice'); +% md=qmu(md,'analysis_type','thermal','sub_analysis_type','transient','package','ice'); +% md=qmu(md,'analysis_type','thermal','sub_analysis_type','steady','package','cielo'); +% md=qmu(md,'analysis_type','thermal','package','cielo'); +% +% On top of these solution arguments (found also in solve.m), user can specify QMU +% specific arguments: +% +% qmudir: any directory where to run the qmu analysis +% qmufile: input file for Dakota +% ivar: selection number for variables input (if several are specified in variables) +% iresp: same thing for response functions +% imethod: same thing for methods +% iparams: same thing for params +% overwrite: overwrite qmudir +% outfiles: (John?) +% rstfile: backup file name +% rundakota: (John?) +% runmpi: (John?) + +%some checks on list of arguments +global ISSM_DIR; + +%recover options +options=recover_solve_options(md,varargin{:}); + +%add default options +options=process_solve_options(options); + +%recover some fields +md.analysis_type=options.analysis_type; +md.sub_analysis_type=options.sub_analysis_type; +package=options.package; + +%Use package to set solution namespace +usenamespace(package); + +if ~ismodelselfconsistent(md,package), + error(' '); %previous error messages should explain what is going on. +end + +displaystring(md.debug,'\n%s\n','launching solution sequence'); + +%Ok, are we running serially or in parallel? +if ~strcmpi(md.cluster,'none'), + %ok, we are running parallel. + md=qmuparallel(md,varargin{:}); +else + %we are running serially. + md=qmuserial(md,varargin{:}); +end Index: /issm/trunk/src/m/classes/public/qmuparallel.m =================================================================== --- /issm/trunk/src/m/classes/public/qmuparallel.m (revision 586) +++ /issm/trunk/src/m/classes/public/qmuparallel.m (revision 586) @@ -0,0 +1,38 @@ +function md=qmuparallel(md,varargin) +%QMUPARALLEL - drive qmu analysis in parallel +% +% Usage: +% md=qmuparallel(md); + +%Get cluster.rc location +cluster_rc_location=which('cluster.rc'); + +%let solution sequences know that we are running in qmu mode +md.qmu_analysis=1; + +%create dakota input file +md=qmuin(md,varargin{:}); + +%Figure out parameters for this particular cluster +[codepath,executionpath]=ClusterParameters(md.cluster,cluster_rc_location); + +%Marshall model data into a binary file. +marshall(md); + +%Create dakota input file. + +%Now, we need to build the queuing script, used by the cluster to launch the job. +BuildQueueingScript(md,executionpath,codepath); + +%Now, launch the queueing script +LaunchQueueJob(md,executionpath); + +%Do we return, or just wait for results? +if md.waitonlock, + %we wait for the done file + waitonlock([executionpath '/' md.name '.lock']); + %load results + md=loadresultsfromcluster(md); +else + return; +end Index: /issm/trunk/src/m/solutions/cielo/prognostic.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/prognostic.m (revision 586) +++ /issm/trunk/src/m/solutions/cielo/prognostic.m (revision 586) @@ -0,0 +1,40 @@ +function md=prognostic(md) +%THERMAL - prognostic solution sequence. +% +% Usage: +% md=prognostic(md) + + %timing + t1=clock; + + %Build all models requested for diagnostic simulation + displaystring(md.debug,'%s',['reading prognostic model data']); + md.analysis_type='prognostic'; m_p=CreateFemModel(md); + + displaystring(debug,'\n%s',['depth averaging velocity...']); + u_g=VelocityDepthAverage(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,m_p.parameters.u_g); + + % figure out number of dof: just for information purposes. + md.dof=m_p.nodesets.fsize; + + %initialize inputs + displaystring(md.debug,'\n%s',['setup inputs...']); + inputs=inputlist; + inputs=add(inputs,'velocity_average',u_g,'doublevec',3,m_p.parameters.numberofnodes); + inputs=add(inputs,'thickness',m_p.parameters.h_g,'doublevec',1,m_p.parameters.numberofnodes); + inputs=add(inputs,'melting',m_p.parameters.m_g,'doublevec',1,m_p.parameters.numberofnodes); + inputs=add(inputs,'accumulation',m_p.parameters.a_g,'doublevec',1,m_p.parameters.numberofnodes); + inputs=add(inputs,'dt',m_p.parameters.dt,'double'); + + displaystring(md.debug,'\n%s',['computing new thickness...']); + h_g=diagnostic_core_linear(m_p,inputs,'prognostic',''); + + displaystring(md.debug,'\n%s',['extrude new thickness...']); + h_g=ThicknessExtrude(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,h_g); + + displaystring(md.debug,'\n%s',['load results...']); + md.new_thickness=h_g; + + %stop timing + t2=clock; + displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/dakota/dakota_m_write.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_m_write.m (revision 585) +++ /issm/trunk/src/m/solutions/dakota/dakota_m_write.m (revision 586) @@ -208,5 +208,4 @@ disp(sprintf(' Writing %d %s responses.',length(dresp),class(dresp))); - for i=1:length(dresp) ifnvals=ifnvals+1; Index: /issm/trunk/src/m/solutions/dakota/qmuin.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/qmuin.m (revision 586) +++ /issm/trunk/src/m/solutions/dakota/qmuin.m (revision 586) @@ -0,0 +1,83 @@ +function md=qmuin(md,package,varargin) +%INPUT function md=qmuin(md,package,varargin) +%Create dakota input file for parallel runs + +global ISSM_DIR; + +% qmudir =['qmu_' datestr(now,'yyyymmdd_HHMMSS')]; +qmudir ='qmu'; +% qmufile can not be changed unless cielo_ice_script.sh is also changed +qmufile='qmu'; +ivar =1; +iresp =1; +imethod=1; +iparams=1; +runmpi =false; + +% process any extra input arguments + +for i=1:2:nargin-2 + switch varargin{i} + case 'qmudir' + qmudir =varargin{i+1}; + disp(sprintf('qmudir =%s',qmudir)); + case 'qmufile' + qmufile=varargin{i+1}; + disp(sprintf('qmufile=%s',qmufile)); + case 'ivar' + ivar =varargin{i+1}; + disp(sprintf('ivar =%d',ivara)); + case 'iresp' + iresp =varargin{i+1}; + disp(sprintf('iresp =%d',iresp)); + case 'imethod' + imethod=varargin{i+1}; + disp(sprintf('imethod=%d',imethod)); + case 'iparams' + iparams=varargin{i+1}; + disp(sprintf('iparams=%d',iparams)); + case 'overwrite' + overwrite=varargin{i+1}; + disp(sprintf('overwrite=%s',overwrite)); + case 'outfiles' + outfiles =varargin{i+1}; + disp(sprintf('outfiles =%s',outfiles)); + case 'rstfile' + rstfile =varargin{i+1}; + disp(sprintf('rstfile =%s',rstfile)); + case 'rundakota' + rundakota=varargin{i+1}; + disp(sprintf('rundakota=%s',rundakota)); + case 'runmpi' + runmpi =varargin{i+1}; + disp(sprintf('runmpi =%d',runmpi)); + end +end + +if exist(qmudir,'dir') + if ~exist('overwrite','var') + overwrite=input(['Overwrite existing ''' qmudir ''' directory? Y/N [N]: '], 's'); + end + if strncmpi(overwrite,'y',1) + system(['rm -rf ' qmudir]); +% else +% error('Existing ''%s'' directory not overwritten.',qmudir); + end +end +mkdir(qmudir) +cd(qmudir) + +%create m and in files for dakota +if (~isfield(md.qmu_params(iparams),'direct') || ... + ~md.qmu_params(iparams).direct) && ... + (~isfield(md.qmu_params(iparams),'analysis_driver') || ... + isempty(md.qmu_params(iparams).analysis_driver)) + md.qmu_params(iparams).analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh']; +end +dakota_in_data(md.qmu_method(imethod),md.variables(ivar),md.responses(iresp),md.qmu_params(iparams),qmufile,package,md); + +system(['rm -rf qmu.m']); +cd ../ + +%keep iresp for response function handling in marshall routine. +md.iresp=iresp; Index: /issm/trunk/src/m/solutions/dakota/qmuserial.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/qmuserial.m (revision 586) +++ /issm/trunk/src/m/solutions/dakota/qmuserial.m (revision 586) @@ -0,0 +1,183 @@ +function md=qmuserial(md,package,varargin) +%INPUT function md=qmu(md,package) +%Deal with coupled ISSM or Cielo/ Dakota runs, to do sensitivity analyses. + +global ISSM_DIR; + +% qmudir =['qmu_' datestr(now,'yyyymmdd_HHMMSS')]; +qmudir ='qmu'; +% qmufile can not be changed unless cielo_ice_script.sh is also changed +qmufile='qmu'; +ivar =1; +iresp =1; +imethod=1; +iparams=1; +runmpi =false; + +% process any extra input arguments + +for i=1:2:nargin-2 + switch varargin{i} + case 'qmudir' + qmudir =varargin{i+1}; + disp(sprintf('qmudir =%s',qmudir)); + case 'qmufile' + qmufile=varargin{i+1}; + disp(sprintf('qmufile=%s',qmufile)); + case 'ivar' + ivar =varargin{i+1}; + disp(sprintf('ivar =%d',ivara)); + case 'iresp' + iresp =varargin{i+1}; + disp(sprintf('iresp =%d',iresp)); + case 'imethod' + imethod=varargin{i+1}; + disp(sprintf('imethod=%d',imethod)); + case 'iparams' + iparams=varargin{i+1}; + disp(sprintf('iparams=%d',iparams)); + case 'overwrite' + overwrite=varargin{i+1}; + disp(sprintf('overwrite=%s',overwrite)); + case 'outfiles' + outfiles =varargin{i+1}; + disp(sprintf('outfiles =%s',outfiles)); + case 'rstfile' + rstfile =varargin{i+1}; + disp(sprintf('rstfile =%s',rstfile)); + case 'rundakota' + rundakota=varargin{i+1}; + disp(sprintf('rundakota=%s',rundakota)); + case 'runmpi' + runmpi =varargin{i+1}; + disp(sprintf('runmpi =%d',runmpi)); + end +end + +%first create temporary directory in which we will work +if exist(qmudir,'dir') + if ~exist('overwrite','var') + overwrite=input(['Overwrite existing ''' qmudir ''' directory? Y/N [N]: '], 's'); + end + if strncmpi(overwrite,'y',1) + system(['rm -rf ' qmudir]); +% else +% error('Existing ''%s'' directory not overwritten.',qmudir); + end +end +mkdir(qmudir) +cd(qmudir) +system('cp $ISSM_DIR/startup.m .'); + +%save our model in qmu so that it can be repeatedly used by Dakota. +save('Qmu.model','md') + +%create m and in files for dakota +if (~isfield(md.qmu_params(iparams),'direct') || ... + ~md.qmu_params(iparams).direct) && ... + (~isfield(md.qmu_params(iparams),'analysis_driver') || ... + isempty(md.qmu_params(iparams).analysis_driver)) + md.qmu_params(iparams).analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh']; +end +dakota_in_data(md.qmu_method(imethod),md.variables(ivar),md.responses(iresp),md.qmu_params(iparams),qmufile,package,md); + +% check for existence of results.out files to use +if exist('results.out.1','file') || exist('results.out.zip','file') + if ~exist('outfiles','var') + outfiles=input(['Use existing ''results.out'' files? Y/N [N]: '], 's'); + end + if ~strncmpi(outfiles,'y',1) + system('rm -f results.out.[1-9]*'); + else + if exist('results.out.zip','file') && ~exist('results.out.1','file') + display('Inflating ''results.out.zip'' file.'); + system('unzip -q results.out.zip'); + end + end +end + +% check for existence of dakota.rst file to use +rstflag=''; +if (~exist('outfiles','var') || ~strncmpi(outfiles,'y',1)) && ... + exist('dakota.rst','file') + if ~exist('rstfile','var') + rstfiles=input(['Use existing ''dakota.rst'' file? Y/N [N]: '], 's'); + end + if strncmpi(rstfiles,'y',1) + system('rm -f results.out.[1-9]*'); + system('dakota_restart_util print dakota.rst | grep completed'); + rstflag=' -read_restart dakota.rst'; + end +end + +%call dakota +if ~exist('rundakota','var') + rundakota=input(['Run Dakota analysis ''' qmufile '''? Y/N [N]: '], 's'); +end +if ~strncmpi(rundakota,'y',1) + cd .. + return +end + +if ~runmpi + system(['dakota -i ' qmufile '.in -o ' qmufile '.out -e ' qmufile '.err' rstflag]); +else +% use 'mpd --ncpus=8 &' to initialize mpi and 'mpdringtest' to verify. +% exporting MPIRUN_NPROCS sets mpi in dakota. +% system('mpd --ncpus=8 &'); +% system('mpdringtest'); + system(['export MPIRUN_NPROCS=8;mpirun -np 4 dakota -i ' qmufile '.in -o ' qmufile '.out -e ' qmufile '.err' rstflag]); +end + +% check to see if dakota returned errors in the err file +if exist([qmufile '.err'],'file') + fide=fopen([qmufile '.err'],'r'); + fline=fgetl(fide); + if ischar(fline) + while ischar(fline) + disp(sprintf('%s',fline)); + fline=fgetl(fide); + end + status=fclose(fide); + cd ../ + error(['Dakota returned error in ''' qmufile '.err'' file. ''' qmudir ''' directory retained.']) + end + status=fclose(fide); +else + cd ../ + error(['Dakota did not generate ''' qmufile '.err'' file. ''' qmudir ''' directory retained.']) +end + +%parse inputs and results from dakota +[method,dvar,dresp_in]=dakota_in_parse([qmufile '.in']); +md.dakotaresults.method =method; +md.dakotaresults.dvar =dvar; +md.dakotaresults.dresp_in =dresp_in; + +[method,dresp_out,scm,pcm,srcm,prcm]=dakota_out_parse([qmufile '.out']); +md.dakotaresults.dresp_out=dresp_out; +md.dakotaresults.scm =scm; +md.dakotaresults.pcm =pcm; +md.dakotaresults.srcm =srcm; +md.dakotaresults.prcm =prcm; + +if exist('dakota_tabular.dat','file') + [method,dresp_dat ]=dakota_out_parse('dakota_tabular.dat'); + md.dakotaresults.dresp_dat=dresp_dat; +end + +%save input and output files into model +%md.dakotain =readfile([qmufile '.in']); +%md.dakotaout=readfile([qmufile '.out']); +%if exist('dakota_tabular.dat','file') +% md.dakotadat=readfile('dakota_tabular.dat'); +%end + +% move all the individual function evalutations into zip files +system('zip -mq params.in.zip params.in.[1-9]*'); +system('zip -mq results.out.zip results.out.[1-9]*'); +system('zip -mq matlab.out.zip matlab*.out.[1-9]*'); + +%get out of local directory and erase +cd ../ +% system(['rm -rf ' qmudir]); Index: /issm/trunk/src/m/utils/Nightly/testsgetpackage.m =================================================================== --- /issm/trunk/src/m/utils/Nightly/testsgetpackage.m (revision 585) +++ /issm/trunk/src/m/utils/Nightly/testsgetpackage.m (revision 586) @@ -25,5 +25,5 @@ elseif strcmpi(string,'cielo_parallel'), package='cielo'; - md.cluster='wilkes'; + md.cluster='astrid'; else Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 585) +++ /issm/trunk/src/mex/Makefile.am (revision 586) @@ -38,9 +38,11 @@ SystemMatrices\ Test\ + ThicknessExtrude\ TriMesh\ TriMeshProcessRifts\ TriMeshRefine\ UpdateFromInputs\ - VelocityExtrude + VelocityExtrude\ + VelocityDepthAverage endif @@ -167,4 +169,7 @@ SystemMatrices/SystemMatrices.h +ThicknessExtrude_SOURCES = ThicknessExtrude/ThicknessExtrude.cpp\ + ThicknessExtrude/ThicknessExtrude.h + TriMesh_SOURCES = TriMesh/TriMesh.cpp\ TriMesh/TriMesh.h @@ -182,2 +187,4 @@ VelocityExtrude/VelocityExtrude.h +VelocityDepthAverage_SOURCES = VelocityDepthAverage/VelocityDepthAverage.cpp\ + VelocityDepthAverage/VelocityDepthAverage.h Index: /issm/trunk/src/mex/ThicknessExtrude/ThicknessExtrude.cpp =================================================================== --- /issm/trunk/src/mex/ThicknessExtrude/ThicknessExtrude.cpp (revision 586) +++ /issm/trunk/src/mex/ThicknessExtrude/ThicknessExtrude.cpp (revision 586) @@ -0,0 +1,54 @@ +/*\file ThicknessExtrude.c + *\brief: extrude thickness vertically + */ + +#include "./ThicknessExtrude.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*diverse: */ + int noerr=1; + + /*input datasets: */ + DataSet* elements=NULL; + DataSet* nodes=NULL; + DataSet* loads=NULL; + DataSet* materials=NULL; + Vec tg=NULL; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&ThicknessExtrudeUsage); + + /*Input datasets: */ + FetchData((void**)&elements,NULL,NULL,ELEMENTS,"DataSet",NULL); + FetchData((void**)&nodes,NULL,NULL,NODES,"DataSet",NULL); + FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL); + FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL); + FetchData((void**)&tg,NULL,NULL,TG,"Vector",NULL); + + /*!Call core code: */ + ThicknessExtrudex(tg,elements,nodes,loads,materials); + + /*write output : */ + WriteData(TGOUT,tg,0,0,"Vector",NULL); + /*Free ressources: */ + delete elements; + delete nodes; + delete loads; + delete materials; + VecFree(&tg); + + /*end module: */ + MODULEEND(); + +} + +void ThicknessExtrudeUsage(void) +{ + _printf_("\n"); + _printf_(" usage: [tg] = %s(elements, nodes,loads, materials, tg);\n",__FUNCT__); + _printf_("\n"); +} Index: /issm/trunk/src/mex/ThicknessExtrude/ThicknessExtrude.h =================================================================== --- /issm/trunk/src/mex/ThicknessExtrude/ThicknessExtrude.h (revision 586) +++ /issm/trunk/src/mex/ThicknessExtrude/ThicknessExtrude.h (revision 586) @@ -0,0 +1,35 @@ + +/* + ThicknessExtrude.h +*/ + + +#ifndef _THICKNESSEXTRUDE_H +#define _THICKNESSEXTRUDE_H + +/* local prototypes: */ +void ThicknessExtrudeUsage(void); + +#include "../../c/issm.h" + +#undef __FUNCT__ +#define __FUNCT__ "ThicknessExtrude" + +/* serial input macros: */ +#define ELEMENTS (mxArray*)prhs[0] +#define NODES (mxArray*)prhs[1] +#define LOADS (mxArray*)prhs[2] +#define MATERIALS (mxArray*)prhs[3] +#define TG (mxArray*)prhs[4] + +/* serial output macros: */ +#define TGOUT (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 5 + + +#endif /* _THICKNESSEXTRUDE_H */ Index: /issm/trunk/src/mex/VelocityDepthAverage/VelocityDepthAverage.cpp =================================================================== --- /issm/trunk/src/mex/VelocityDepthAverage/VelocityDepthAverage.cpp (revision 586) +++ /issm/trunk/src/mex/VelocityDepthAverage/VelocityDepthAverage.cpp (revision 586) @@ -0,0 +1,54 @@ +/*\file VelocityDepthAverage.c + *\brief: average velocity through thickness + */ + +#include "./VelocityDepthAverage.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*diverse: */ + int noerr=1; + + /*input datasets: */ + DataSet* elements=NULL; + DataSet* nodes=NULL; + DataSet* loads=NULL; + DataSet* materials=NULL; + Vec ug=NULL; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&VelocityDepthAverageUsage); + + /*Input datasets: */ + FetchData((void**)&elements,NULL,NULL,ELEMENTS,"DataSet",NULL); + FetchData((void**)&nodes,NULL,NULL,NODES,"DataSet",NULL); + FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL); + FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL); + FetchData((void**)&ug,NULL,NULL,UG,"Vector",NULL); + + /*!Call core code: */ + VelocityDepthAveragex(ug,elements,nodes,loads,materials); + + /*write output : */ + WriteData(UGOUT,ug,0,0,"Vector",NULL); + /*Free ressources: */ + delete elements; + delete nodes; + delete loads; + delete materials; + VecFree(&ug); + + /*end module: */ + MODULEEND(); + +} + +void VelocityDepthAverageUsage(void) +{ + _printf_("\n"); + _printf_(" usage: [ug] = %s(elements, nodes,loads, materials, ug);\n",__FUNCT__); + _printf_("\n"); +} Index: /issm/trunk/src/mex/VelocityDepthAverage/VelocityDepthAverage.h =================================================================== --- /issm/trunk/src/mex/VelocityDepthAverage/VelocityDepthAverage.h (revision 586) +++ /issm/trunk/src/mex/VelocityDepthAverage/VelocityDepthAverage.h (revision 586) @@ -0,0 +1,35 @@ + +/* + VelocityDepthAverage.h +*/ + + +#ifndef _VELOCITYDEPTHAVERAGE_H +#define _VELOCITYDEPTHAVERAGE_H + +/* local prototypes: */ +void VelocityDepthAverageUsage(void); + +#include "../../c/issm.h" + +#undef __FUNCT__ +#define __FUNCT__ "VelocityDepthAverage" + +/* serial input macros: */ +#define ELEMENTS (mxArray*)prhs[0] +#define NODES (mxArray*)prhs[1] +#define LOADS (mxArray*)prhs[2] +#define MATERIALS (mxArray*)prhs[3] +#define UG (mxArray*)prhs[4] + +/* serial output macros: */ +#define UGOUT (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 5 + + +#endif /* _VELOCITYDEPTHAVERAGE_H */