Changeset 13721


Ignore:
Timestamp:
10/17/12 15:59:20 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moved MassFluxx from modules to FemModel method

Location:
issm/trunk-jpl/src/c
Files:
1 deleted
4 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/Makefile.am

    r13709 r13721  
    521521                                             ./modules/IceVolumex/IceVolumex.cpp\
    522522                                             ./modules/ElementResponsex/ElementResponsex.h\
    523                                              ./modules/ElementResponsex/ElementResponsex.cpp\
    524                                              ./modules/MassFluxx/MassFluxx.cpp\
    525                                              ./modules/MassFluxx/MassFluxx.h
     523                                             ./modules/ElementResponsex/ElementResponsex.cpp
    526524#}}}
    527525#Slope sources  {{{
  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r13700 r13721  
    228228}
    229229/*}}}*/
     230/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/
     231void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){
     232
     233        /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use
     234         * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several
     235         * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the
     236         * Slope configuration.*/
     237
     238        int found=-1;
     239        for(int i=0;i<nummodels;i++){
     240                if (analysis_type_list[i]==configuration_type){
     241                        found=i;
     242                        break;
     243                }
     244        }
     245        if(found!=-1) analysis_counter=found;
     246        else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");
     247
     248        /*Now, plug analysis_counter and analysis_type inside the parameters: */
     249        this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
     250        this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
     251        this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);
     252
     253        /*configure elements, loads and nodes, for this new analysis: */
     254        this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
     255        this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
     256        this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
     257
     258        #ifdef _HAVE_PETSC_
     259        /*take care of petsc options, that depend on this analysis type (present only after model processor)*/
     260        if(this->parameters->Exist(PetscOptionsStringsEnum)){
     261                PetscOptionsFromAnalysis(this->parameters,analysis_type);
     262                if(VerboseSolver()) _pprintLine_("      petsc Options set for analysis type: " << EnumToStringx(analysis_type));
     263        }
     264        #endif
     265
     266}
     267/*}}}*/
     268/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/
     269void FemModel::SetCurrentConfiguration(int configuration_type){
     270        this->SetCurrentConfiguration(configuration_type,configuration_type);
     271}
     272/*}}}*/
    230273/*FUNCTION FemModel::Solve {{{*/
    231274void FemModel::Solve(void){
     
    280323/*}}}*/
    281324
    282 /*Numerics: */
    283 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/
    284 void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){
    285 
    286         /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use
    287          * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several
    288          * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the
    289          * Slope configuration.*/
    290 
    291         int found=-1;
    292         for(int i=0;i<nummodels;i++){
    293                 if (analysis_type_list[i]==configuration_type){
    294                         found=i;
    295                         break;
    296                 }
    297         }
    298         if(found!=-1) analysis_counter=found;
    299         else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");
    300 
    301         /*Now, plug analysis_counter and analysis_type inside the parameters: */
    302         this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
    303         this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
    304         this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);
    305 
    306         /*configure elements, loads and nodes, for this new analysis: */
    307         this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
    308         this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
    309         this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
    310 
    311         #ifdef _HAVE_PETSC_
    312         /*take care of petsc options, that depend on this analysis type (present only after model processor)*/
    313         if(this->parameters->Exist(PetscOptionsStringsEnum)){
    314                 PetscOptionsFromAnalysis(this->parameters,analysis_type);
    315                 if(VerboseSolver()) _pprintLine_("      petsc Options set for analysis type: " << EnumToStringx(analysis_type));
    316         }
    317         #endif
    318 
    319 }
    320 /*}}}*/
    321 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/
    322 void FemModel::SetCurrentConfiguration(int configuration_type){
    323         this->SetCurrentConfiguration(configuration_type,configuration_type);
    324 }
    325 /*}}}*/
    326 
    327325/*Modules:*/
    328 int FemModel::UpdateVertexPositionsx(void){ /*{{{*/
     326int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
    329327
    330328        int     i;
     
    378376        NodesDofx(nodes,parameters,analysis_type);
    379377
    380 }
    381 /*}}}*/
    382 void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
    383 
    384         int      i;
    385 
    386         /*output: */
    387         IssmDouble   dt;
    388 
    389         /*intermediary: */
    390         Element *element     = NULL;
    391         IssmDouble   min_dt      = 0;
    392         IssmDouble   node_min_dt = 0;
    393 
    394         /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */
    395         element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();
    396 
    397         for (i=1;i<elements->Size();i++){
    398                 element=(Element*)elements->GetObjectByOffset(i);
    399                 dt=element->TimeAdapt();
    400                 if(dt<min_dt)min_dt=dt;
    401         }
    402 
    403         /*Figure out minimum across the cluster: */
    404         #ifdef _HAVE_MPI_
    405         MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );
    406         MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());
    407         min_dt=node_min_dt;
    408         #endif
    409 
    410         /*Assign output pointers:*/
    411         *pdt=min_dt;
    412378}
    413379/*}}}*/
     
    440406                case MaxVzEnum:                  MaxVzx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    441407                case MaxAbsVzEnum:               MaxAbsVzx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    442                 case MassFluxEnum:               MassFluxx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     408                case MassFluxEnum:               this->MassFluxx(          responses,process_units); break;
    443409                case SurfaceAbsVelMisfitEnum:    SurfaceAbsVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
    444410                case SurfaceRelVelMisfitEnum:    SurfaceRelVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
     
    543509}
    544510/*}}}*/
     511void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
     512
     513        int      i;
     514
     515        /*output: */
     516        IssmDouble   dt;
     517
     518        /*intermediary: */
     519        Element *element     = NULL;
     520        IssmDouble   min_dt      = 0;
     521        IssmDouble   node_min_dt = 0;
     522
     523        /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */
     524        element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();
     525
     526        for (i=1;i<elements->Size();i++){
     527                element=(Element*)elements->GetObjectByOffset(i);
     528                dt=element->TimeAdapt();
     529                if(dt<min_dt)min_dt=dt;
     530        }
     531
     532        /*Figure out minimum across the cluster: */
     533#ifdef _HAVE_MPI_
     534        MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );
     535        MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());
     536        min_dt=node_min_dt;
     537#endif
     538
     539        /*Assign output pointers:*/
     540        *pdt=min_dt;
     541}
     542/*}}}*/
     543#ifdef _HAVE_RESPONSES_
     544void FemModel::MassFluxx(IssmDouble* pmass_flux,bool process_units){/*{{{*/
     545
     546        int          i,j;
     547        Element     *element       = NULL;
     548        int          element_id;
     549        bool         ispresent     = false;
     550        IssmDouble   mass_flux     = 0;
     551        IssmDouble   all_mass_flux = 0;
     552        int          counter;
     553        IssmDouble **array         = NULL;
     554        int          M;
     555        int         *mdims_array   = NULL;
     556        int         *ndims_array   = NULL;
     557        IssmDouble  *segments      = NULL;
     558        int          num_segments;
     559
     560        /*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */
     561        this->parameters->FindParam(&ispresent,MassFluxSegmentsPresentEnum);
     562        if(!ispresent)_error_("no mass flux segments available!");
     563        this->parameters->FindParam(&array,&M,&mdims_array,&ndims_array,MassFluxSegmentsEnum);
     564
     565        /*Retrieve index of segments being used for MassFlux computation: */
     566        parameters->FindParam(&counter,IndexEnum);
     567
     568        /*retrieve segments from array: */
     569        segments     = array[counter-1]; //matlab to "C" indexing
     570        num_segments = mdims_array[counter-1];
     571
     572        /*Go through segments, and then elements, and figure out which elements belong to a segment.
     573         * When we find one, use the element to compute the mass flux on the segment: */
     574        for(i=0;i<num_segments;i++){
     575                element_id=reCast<int,IssmDouble>(*(segments+5*i+4));
     576                for(j=0;j<elements->Size();j++){
     577                        element=(Element*)this->elements->GetObjectByOffset(j);
     578                        if (element->Id()==element_id){
     579                                /*We found the element which owns this segment, use it to compute the mass flux: */
     580                                mass_flux+=element->MassFlux(segments+5*i+0,process_units);
     581                                break;
     582                        }
     583                }
     584        }
     585
     586#ifdef _HAVE_MPI_
     587        MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
     588        mass_flux=all_mass_flux;
     589#endif
     590
     591        /*Free ressources:*/
     592        for(i=0;j<M;i++){
     593                IssmDouble* matrix=array[j];
     594                xDelete<IssmDouble>(matrix);
     595        }
     596        xDelete<int>(mdims_array);
     597        xDelete<int>(ndims_array);
     598        xDelete<IssmDouble*>(array);
     599
     600        /*Assign output pointers: */
     601        *pmass_flux=mass_flux;
     602
     603}/*}}}*/
     604#endif
    545605#ifdef _HAVE_CONTROL_
    546606void FemModel::ThicknessAbsGradientx( IssmDouble* pJ, bool process_units, int weight_index){/*{{{*/
     
    601661/*}}}*/
    602662#endif
    603 
    604663#ifdef  _HAVE_DAKOTA_
    605664void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
  • issm/trunk-jpl/src/c/classes/FemModel.h

    r13700 r13721  
    5050                FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
    5151                ~FemModel();
     52
     53                /*Methods:*/
     54                void Echo();
    5255                void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
    53 
    54                 /*Methods: {{{*/
    55                 void Echo();
    5656                void Solve(void);
    5757                void OutputResults(void);
    5858                void SetStaticComm();
    59                 /*}}}*/
    60                 /*Fem: {{{*/
    61                 void  SetCurrentConfiguration(int configuration_type);
    62                 void  SetCurrentConfiguration(int configuration_type,int analysis_type);
    63                 /*}}}*/
    64                 /*Modules: {{{*/
    65                  #ifdef  _HAVE_DAKOTA_
     59                void SetCurrentConfiguration(int configuration_type);
     60                void SetCurrentConfiguration(int configuration_type,int analysis_type);
     61
     62                /*Modules*/
     63                #ifdef _HAVE_RESPONSES_
     64                void MassFluxx(IssmDouble* pmass_flux,bool process_units);
     65                #endif
     66                #ifdef  _HAVE_DAKOTA_
    6667                void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
    6768                #endif
    68 
    6969                void RequestedOutputsx(int* requested_outputs, int numoutputs);
    7070                void RequestedDependentsx(void);
     
    7878                void UpdateConstraintsx(void);
    7979                int  UpdateVertexPositionsx(void);
    80 
    81 
    82                 /*}}}*/
    83 
    8480};
    8581
  • issm/trunk-jpl/src/c/modules/modules.h

    r13700 r13721  
    6363#include "./Krigingx/Krigingx.h"
    6464#include "./Shp2Kmlx/Shp2Kmlx.h"
    65 #include "./MassFluxx/MassFluxx.h"
    6665#include "./MaxAbsVxx/MaxAbsVxx.h"
    6766#include "./MaxAbsVyx/MaxAbsVyx.h"
Note: See TracChangeset for help on using the changeset viewer.