Changeset 8429


Ignore:
Timestamp:
05/25/11 17:19:53 (14 years ago)
Author:
Mathieu Morlighem
Message:

merged the transients 2D and 3D

Location:
issm/trunk/src/c
Files:
1 added
2 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r8399 r8429  
    4949        HydrologyAnalysisEnum,
    5050        //transient
    51         Transient2DSolutionEnum,
    52         Transient2DAnalysisEnum,
    53         Transient3DSolutionEnum,
    54         Transient3DAnalysisEnum,
     51        TransientSolutionEnum,
    5552        GroundingLineMigration2DSolutionEnum,
    5653        SteadyAnalysisEnum,

  • issm/trunk/src/c/Makefile.am

    r8410 r8429  
    13121312                                        ./solutions/hydrology_core.cpp\
    13131313                                        ./solutions/hydrology_core_step.cpp\
    1314                                         ./solutions/transient2d_core.cpp\
    1315                                         ./solutions/transient3d_core.cpp\
     1314                                        ./solutions/transient_core.cpp\
    13161315                                        ./solutions/groundinglinemigration2d_core.cpp\
    13171316                                        ./solutions/steadystate_core.cpp\

  • issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r8399 r8429  
    4444                case HydrologySolutionEnum : return "HydrologySolution";
    4545                case HydrologyAnalysisEnum : return "HydrologyAnalysis";
    46                 case Transient2DSolutionEnum : return "Transient2DSolution";
    47                 case Transient2DAnalysisEnum : return "Transient2DAnalysis";
    48                 case Transient3DSolutionEnum : return "Transient3DSolution";
    49                 case Transient3DAnalysisEnum : return "Transient3DAnalysis";
     46                case TransientSolutionEnum : return "TransientSolution";
    5047                case GroundingLineMigration2DSolutionEnum : return "GroundingLineMigration2DSolution";
    5148                case SteadyAnalysisEnum : return "SteadyAnalysis";

  • issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp

    r8365 r8429  
    1818void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,FILE* iomodel_handle,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){
    1919
    20         bool       continuous = true;
    21         Elements  *elements   = NULL;
    22         Materials *materials  = NULL;
    23         Parameters *parameters  = NULL;
     20        bool        continuous = true;
     21        Elements   *elements   = NULL;
     22        Materials  *materials  = NULL;
     23        Parameters *parameters = NULL;
    2424                       
    2525        /*Create elements, vertices and materials, independent of analysis_type: */
    … …  
    3232        /*Now, branch onto analysis dependent model generation: */
    3333        switch(analysis_type){
    34                 case  DiagnosticHorizAnalysisEnum:
     34                case DiagnosticHorizAnalysisEnum:
    3535                        CreateNodesDiagnosticHoriz(pnodes, iomodel,iomodel_handle);
    3636                        CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle);
    … …  
    3939                        break;
    4040               
    41                 case  DiagnosticVertAnalysisEnum:
     41                case DiagnosticVertAnalysisEnum:
    4242                        CreateNodesDiagnosticVert(pnodes, iomodel,iomodel_handle);
    4343                        CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle);
    … …  
    120120        /*Update Elements in case we are running a transient solution: */
    121121        parameters=*pparameters;
    122         if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){
     122        if(analysis_counter==(nummodels-1)&& solution_type==TransientSolutionEnum){
    123123                UpdateElementsTransient(elements,parameters,iomodel,iomodel_handle,analysis_counter,analysis_type);
    124124        }
    … …  
    130130         * will need to start at the end of the updated counters: */
    131131        UpdateCounters(iomodel,pnodes,ploads,pconstraints);
    132 
    133132}

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

    r8330 r8429  
    2222       
    2323        /*output: */
    24         Elements* elements=NULL;
    25         Nodes* nodes=NULL;
    26         Vertices* vertices=NULL;
    27         Materials* materials=NULL;
    28         Constraints* constraints=NULL;
    29         Loads* loads=NULL;
    30         Parameters* parameters=NULL;
     24        Elements    *elements    = NULL;
     25        Nodes       *nodes       = NULL;
     26        Vertices    *vertices    = NULL;
     27        Materials   *materials   = NULL;
     28        Constraints *constraints = NULL;
     29        Loads       *loads       = NULL;
     30        Parameters  *parameters  = NULL;
    3131
    32         /*intermediary: */
    33         IoModel* iomodel=NULL;
    34 
    35         iomodel = new IoModel(IOMODEL);
     32        /*Initialize IoModel from input file*/
     33        IoModel* iomodel = new IoModel(IOMODEL);
    3634        SetVerbosityLevel(iomodel->verbose);
    3735
    … …  
    3937
    4038                analysis_type=analysis_type_list[i];
     39
     40                /*Hack for trasient runs (to be improved)*/
     41                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && iomodel->dim==2) continue;
     42                if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && iomodel->dim==2) continue;
    4143       
    4244                _printf_(VerboseMProcessor(),"   create datasets for analysis %s\n",EnumToStringx(analysis_type));

  • issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r8399 r8429  
    4242        else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    4343        else if (strcmp(name,"HydrologyAnalysis")==0) return HydrologyAnalysisEnum;
    44         else if (strcmp(name,"Transient2DSolution")==0) return Transient2DSolutionEnum;
    45         else if (strcmp(name,"Transient2DAnalysis")==0) return Transient2DAnalysisEnum;
    46         else if (strcmp(name,"Transient3DSolution")==0) return Transient3DSolutionEnum;
    47         else if (strcmp(name,"Transient3DAnalysis")==0) return Transient3DAnalysisEnum;
     44        else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    4845        else if (strcmp(name,"GroundingLineMigration2DSolution")==0) return GroundingLineMigration2DSolutionEnum;
    4946        else if (strcmp(name,"SteadyAnalysis")==0) return SteadyAnalysisEnum;

  • issm/trunk/src/c/shared/Matlab/PrintfFunction.cpp

    r7737 r8429  
    66#include "stdarg.h"
    77#include "stdio.h"
     8#include "../../shared/shared.h"
     9#include "../../include/include.h"
    810
    911#ifdef _SERIAL_
    … …  
    1214
    1315int PrintfFunction(char* format,...){
     16        /*http://linux.die.net/man/3/vsnprintf*/
    1417
    1518        /*string to be printed: */
    16         char string[1024]; //assume that nobody will print more than 1024 characters at once.
    17         int  string_size;
    18         extern int my_rank;
    19         extern int num_procs;
     19        char *buffer = NULL;
     20        int   n,size = 100;
     21        int         string_size;
     22        extern int  my_rank;
     23        extern int  num_procs;
    2024                                                                                                                                                                                                     
    2125        //variable list of arguments
    22         va_list ap;
     26        va_list args;
    2327
    24         //First use vsprintf to get the whole input string.
    25         va_start(ap,format);
    26         #ifndef WIN32
    27         string_size=vsprintf(string,format,ap); //printf style coding
    28         #else
    29         string_size=vsprintf_s(string,format,ap); //printf style coding
    30         #endif
    31         va_end(ap);
     28        while(true){
    3229
    33         //Ok, if we are running in paralle, get node 0 to print
    34         #if defined(_PARALLEL_)
    35                 if(my_rank==0)printf(string);
    36         #else
    37                 mexPrintf(string);
    38         #endif
     30                /*allocate buffer for given string size*/
     31                buffer=(char*)xmalloc(size*sizeof(char));
    3932
     33                /* Try to print in the allocated space. */
     34                va_start(args, format);
     35#ifndef WIN32
     36                n=vsnprintf(buffer,size,format,args);
     37#else
     38                n=vsnprintf_s(buffer,size,format,args);
     39#endif
     40                va_end(args);
     41
     42                /* If that worked, return the string. */
     43                if(n>-1 && n<size) break;
     44
     45                /* Else try again with more space. */
     46                if(n>-1)   /* glibc 2.1 */
     47                 size=n+1; /* precisely what is needed */
     48                else       /* glibc 2.0 */
     49                 size*=2;  /* twice the old size */
     50
     51                xfree((void**)&buffer);
     52        }
     53
     54        /*Ok, if we are running in paralle, get node 0 to print*/
     55#if defined(_PARALLEL_)
     56        if(my_rank==0)printf(buffer);
     57#else
     58        mexPrintf(buffer);
     59#endif
     60
     61        /*Clean up and return*/
     62        xfree((void**)&buffer);
    4063        return 1;
    4164}

  • issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp

    r8288 r8429  
    4949                        solutioncore=&bedslope_core;
    5050                        break;
    51                 case Transient2DSolutionEnum:
    52                         solutioncore=&transient2d_core;
    53                         break;
    54                 case Transient3DSolutionEnum:
    55                         solutioncore=&transient3d_core;
     51                case TransientSolutionEnum:
     52                        solutioncore=&transient_core;
    5653                        break;
    5754                case GroundingLineMigration2DSolutionEnum:

  • issm/trunk/src/c/solutions/SolutionConfiguration.cpp

    r8287 r8429  
    1717#include "../include/include.h"
    1818#include "../solvers/solvers.h"
    19 
    2019
    2120void SolutionConfiguration(int** panalyses,int* pnumanalyses, void (**psolutioncore)(FemModel*),int solutiontype){
    … …  
    9998                        break;
    10099
    101                 case Transient2DSolutionEnum:
    102                         numanalyses=6;
    103                         analyses=(int*)xmalloc(numanalyses*sizeof(int));
    104                         analyses[0]=DiagnosticHorizAnalysisEnum;
    105                         analyses[1]=DiagnosticVertAnalysisEnum;
    106                         analyses[2]=DiagnosticHutterAnalysisEnum;
    107                         analyses[3]=SurfaceSlopeAnalysisEnum;
    108                         analyses[4]=BedSlopeAnalysisEnum;
    109                         analyses[5]=PrognosticAnalysisEnum;
    110                         break;
    111                
    112                 case Transient3DSolutionEnum:
     100                case TransientSolutionEnum:
    113101                        numanalyses=8;
    114102                        analyses=(int*)xmalloc(numanalyses*sizeof(int));
    … …  
    118106                        analyses[3]=SurfaceSlopeAnalysisEnum;
    119107                        analyses[4]=BedSlopeAnalysisEnum;
    120                         analyses[5]=PrognosticAnalysisEnum;
    121                         analyses[6]=ThermalAnalysisEnum;
    122                         analyses[7]=MeltingAnalysisEnum;
     108                        analyses[5]=ThermalAnalysisEnum;
     109                        analyses[6]=MeltingAnalysisEnum;
     110                        analyses[7]=PrognosticAnalysisEnum;
    123111                        break;
    124112               
    … …  
    140128        else              xfree((void**)&analyses);
    141129        if(psolutioncore)*psolutioncore=solutioncore;
    142 
    143130}

  • issm/trunk/src/c/solutions/controltao_core.cpp

    r8390 r8429  
    1717        /*TAO*/
    1818        int                i,n,info;
    19         TaoMethod          method     = "tao_blmvm";
    20         //TaoMethod          method     = "tao_lmvm";
    21         //TaoMethod          method     = "tao_cg";
    22         //TaoMethod          method     = "tao_gpcg"; -> Hessian
     19        //TaoMethod          method = "tao_blmvm";
     20        //TaoMethod        method = "tao_lmvm";
     21        TaoMethod        method = "tao_cg";
    2322        TaoTerminateReason reason;
    2423        TAO_SOLVER         tao;
    2524        TAO_APPLICATION    controlapp;
    26         Vec                initial_solution=NULL;
    27         AppCtx          user;                /* user-defined work context */
     25        Vec                initial_solution = NULL;
     26        AppCtx             user;
     27        PetscInt           iter;
     28        double             ff,gnorm;
    2829
    29 #ifdef _HAVE_TAO_
     30        /*Initialize TAO*/
    3031        int argc; char **args; PetscGetArgs(&argc,&args);
    3132        int ierr=TaoInitialize(&argc,&args,(char*)0,"");
    3233        if(ierr) _error_("Could not initialize Tao");
    33 #endif
     34
     35        /*Line search options*/
     36        info = PetscOptionsSetValue("-tao_ls_stepmax","10e11"); if(info) _error_("STOP"); //does not work
     37        info = PetscOptionsSetValue("-tao_ls_stepmin","10e5"); if(info) _error_("STOP");    //does not work
     38        info = PetscOptionsSetValue("-tao_ls_maxfev","8"); if(info) _error_("STOP");
     39
     40        /*TAO options: http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoSetFromOptions.html*/
     41        info = PetscOptionsSetValue("-tao_monitor",""); if(info) _error_("STOP");
     42        info = PetscOptionsSetValue("-tao_gatol","10e-18"); if(info) _error_("STOP");
     43        info = PetscOptionsSetValue("-tao_max_its","10"); if(info) _error_("STOP");
     44        info = PetscOptionsSetValue("-tao_max_funcs","20"); if(info) _error_("STOP");
     45
     46        /*Additional options*/
     47        //info = PetscOptionsSetValue("-info","/u/astrid-r1b/morlighe/svn/issm/trunk/test/NightlyRun/taolog.txt"); if(info) _error_("STOP");
    3448
    3549        /*Initialize argument*/
    3650        user.femmodel=femmodel;
    3751
     52        /*Set up and solve TAO*/
     53        GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,RheologyBbarEnum,VertexEnum);
    3854        info = TaoCreate(PETSC_COMM_WORLD,method,&tao); if(info) _error_("STOP");
    3955        info = TaoApplicationCreate(PETSC_COMM_WORLD,&controlapp); if(info) _error_("STOP");
    40         GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,RheologyBbarEnum,VertexEnum);
    4156        info = TaoAppSetInitialSolutionVec(controlapp,initial_solution);  if(info) _error_("STOP");
    4257        info = TaoAppSetObjectiveAndGradientRoutine(controlapp,FormFunctionGradient,(void*)&user);  if(info) _error_("STOP");
    4358        info = TaoSetOptions(controlapp,tao);  if(info) _error_("STOP");
    44         info = TaoSetTolerances(tao,1e-24,1e-24,1e-24,1e-24); if(info) _error_("STOP");
    45         //info = TaoSetFunctionLowerBound(tao,0.77); if(info) _error_("STOP");
    4659        info = TaoSolveApplication(controlapp,tao); //if(info) _error_("STOP");
    47         PetscInt          iter;
    48         double          ff,gnorm;
     60
     61        /*Get solution status*/
    4962        info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); //CHKERRQ(info);
    50 
    51         switch(reason){
    52                 /*http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoGetTerminationReason.html*/
     63        switch(reason){ /*http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoGetTerminationReason.html*/
    5364                case TAO_CONVERGED_ATOL:       _printf_(true,"TAO_CONVERGED_ATOL (res <= atol)\n"); break;
    5465                case TAO_CONVERGED_RTOL:       _printf_(true,"TAO_CONVERGED_RTOL (res/res0 <= rtol)\n"); break;
    … …  
    7283        }
    7384        info = TaoView(tao);  if(info) _error_("STOP");
     85
     86        /*Clean up*/
    7487        info = TaoDestroy(tao);  if(info) _error_("STOP");
    7588        info = TaoAppDestroy(controlapp);  if(info) _error_("STOP");
    76 
    77         /*Clean up*/
    7889        VecFree(&initial_solution);
    7990
    … …  
    89100        Vec       gradient = NULL;
    90101
     102        /*Temp*/
     103//      double*   Xserial=NULL;
     104//      VecToMPISerial(&Xserial,X);
     105//      printf("X= [%20.20g %20.20g %20.20g]\n",Xserial[0],Xserial[1],Xserial[2]);
     106        /*End Temp*/
     107
    91108        InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum,CmResponseEnum);
    92109        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,RheologyBbarEnum,VertexEnum);
    93110        adjointdiagnostic_core(user->femmodel);
    94111        Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, RheologyBbarEnum);
    95         //VecView(gradient,PETSC_VIEWER_STDOUT_SELF);
    96         //VecScale(gradient,-1.0);
     112        VecScale(gradient,10e7);
     113        //VecScale(gradient,-1.);
    97114        VecCopy(gradient,G);
     115        CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum);
     116
     117        //printf("X\n");
     118        //VecView(X,PETSC_VIEWER_STDOUT_SELF);
     119
     120        //printf("Gradient\n");
    98121        //VecView(G,PETSC_VIEWER_STDOUT_SELF);
    99         CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum);
     122
    100123        printf("f(x) = %g\n",*fcn);
    101124        return 0;

  • issm/trunk/src/c/solutions/issm.cpp

    r8330 r8429  
    1212
    1313        /*I/O: */
    14         FILE     *output_fid       = NULL;
    15         FILE     *petscoptionsfid = NULL;
    16         char     *lockname         = NULL;
    17         bool      qmu_analysis     = false;
    18         bool      control_analysis = false;
    19         bool      waitonlock       = false;
     14        FILE *output_fid       = NULL;
     15        FILE *petscoptionsfid = NULL;
     16        char *lockname         = NULL;
     17        bool  qmu_analysis     = false;
     18        bool  control_analysis = false;
     19        bool  waitonlock       = false;
    2020
    2121        /*FemModel: */
    22         FemModel *femmodel         = NULL;
     22        FemModel *femmodel = NULL;
    2323
    2424        /*configuration: */

  • issm/trunk/src/c/solutions/solutions.h

    r8288 r8429  
    3030void slopecompute_core(FemModel* femmodel);
    3131void steadystate_core(FemModel* femmodel);
    32 void transient2d_core(FemModel* femmodel);
    33 void transient3d_core(FemModel* femmodel);
     32void transient_core(FemModel* femmodel);
    3433void groundinglinemigration2d_core(FemModel* femmodel);
    3534double objectivefunctionC(double search_scalar,OptArgs* optargs);
Note: See TracChangeset for help on using the changeset viewer.