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

CHG: moved MassFluxx from modules to FemModel method

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.