Changeset 18348


Ignore:
Timestamp:
08/08/14 13:18:01 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: trying to push to subfunction

Location:
issm/trunk-jpl/src/c
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r18286 r18348  
    722722        }
    723723        delete gauss;
     724}
     725/*}}}*/
     726void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,int input_enum){/*{{{*/
     727
     728        Input* input = GetInput(input_enum);
     729        if(!input){
     730                _error_("Input "<<EnumToStringx(input_enum)<<" not found");
     731        }
     732
     733        int numnodes = this->GetNumberOfNodes();
     734        IssmDouble input_min = input->Min();
     735        IssmDouble input_max = input->Max();
     736        for(int i=0;i<numnodes;i++){
     737                if(min[nodes[i]->Sid()]>input_min) min[nodes[i]->Sid()] = input_min;
     738                if(max[nodes[i]->Sid()]<input_max) max[nodes[i]->Sid()] = input_max;
     739        }
    724740}
    725741/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r18293 r18348  
    8585                void       GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    8686                void       GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
     87                void       GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,int input_enum);
    8788                void       GetInputValue(bool* pvalue,int enum_type);
    8889                void       GetInputValue(int* pvalue,int enum_type);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r18259 r18348  
    11251125
    11261126}/*}}}*/
     1127void FemModel::GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,int input_enum){/*{{{*/
     1128
     1129        /*Get vector sizes for current configuration*/
     1130        int configuration_type;
     1131        this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     1132        int numnodes = this->nodes->NumberOfNodes(configuration_type);
     1133
     1134        /*Initialize output vectors*/
     1135        IssmDouble* uLmin_local = xNew<IssmDouble>(numnodes);
     1136        IssmDouble* uLmax_local = xNew<IssmDouble>(numnodes);
     1137        IssmDouble* uLmin = xNew<IssmDouble>(numnodes);
     1138        IssmDouble* uLmax = xNew<IssmDouble>(numnodes);
     1139        for(int i=0;i<numnodes;i++){
     1140                uLmin_local[i] = +1.e+50;
     1141                uLmax_local[i] = -1.e+50;
     1142        }
     1143
     1144        for(int i=0;i<this->elements->Size();i++){
     1145                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     1146                element->GetInputLocalMinMaxOnNodes(uLmin_local,uLmax_local,input_enum);
     1147        }
     1148
     1149        /*Synchronize all CPUs*/
     1150        ISSM_MPI_Allreduce((void*)uLmin_local,(void*)uLmin,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
     1151        ISSM_MPI_Allreduce((void*)uLmax_local,(void*)uLmax,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
     1152        xDelete<IssmDouble>(uLmin_local);
     1153        xDelete<IssmDouble>(uLmax_local);
     1154
     1155        /*Assign output pointers: */
     1156        *pmin=uLmin;
     1157        *pmax=uLmax;
     1158
     1159}/*}}}*/
    11271160void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/
    11281161

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r18237 r18348  
    6161
    6262                /*Modules*/
     63                void GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,int input_enum);
    6364                void MassFluxx(IssmDouble* presponse);
    6465                void MaxAbsVxx(IssmDouble* presponse);

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

    r18347 r18348  
    88#include "../modules/modules.h"
    99
     10#ifdef _HAVE_PETSC_
     11void CreateDMatrix(Mat* pD,Mat K){/*{{{*/
     12        /*Create D matrix such that:
     13         *
     14         * d_ij = max( -k_ij,0,-k_ji) off diagonal
     15         *
     16         * d_ii = - sum_{i!=j} d_ij for the diagonal
     17         *
     18         */
     19
     20        /*Intermediaries*/
     21        int        ncols,ncols2,rstart,rend;
     22        double     d,diagD;
     23        Mat        D        = NULL;
     24        Mat        K_transp = NULL;
     25        int*       cols  = NULL;
     26        int*       cols2 = NULL;
     27        double*    vals  = NULL;
     28        double*    vals2 = NULL;
     29
     30        /*First, we need to transpose K so that we access both k_ij and k_ji*/
     31        MatTranspose(K,MAT_INITIAL_MATRIX,&K_transp);
     32
     33        /*Initialize output (D has the same non zero pattern as K)*/
     34        MatDuplicate(K,MAT_SHARE_NONZERO_PATTERN,&D);
     35
     36        /*Go through the rows of K an K' and build D*/
     37        MatGetOwnershipRange(K,&rstart,&rend);
     38        for(int row=rstart; row<rend; row++){
     39                diagD = 0.;
     40                MatGetRow(K       ,row,&ncols, (const int**)&cols, (const double**)&vals);
     41                MatGetRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     42                _assert_(ncols==ncols2);
     43                for(int j=0; j<ncols; j++) {
     44                        _assert_(cols[j]==cols2[j]);
     45                        d = max(max(-vals[j],-vals2[j]),0.);
     46                        MatSetValue(D,row,cols[j],(const double)d,INSERT_VALUES);
     47                        if(cols[j]!=row) diagD -= d;
     48                }
     49                MatSetValue(D,row,row,(const double)diagD,INSERT_VALUES);
     50                MatRestoreRow(K       ,row,&ncols, (const int**)&cols, (const double**)&vals);
     51                MatRestoreRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     52        }
     53        MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);
     54        MatAssemblyEnd(  D,MAT_FINAL_ASSEMBLY);
     55
     56        /*Clean up and assign output pointer*/
     57        MatFree(&K_transp);
     58        *pD = D;
     59}/*}}}*/
     60void CreateLHS(Mat* pLHS,IssmDouble* pdmax,Mat K,Mat D,Vec Ml,IssmDouble theta,IssmDouble deltat,FemModel* femmodel,int configuration_type){/*{{{*/
     61        /*Create Left Hand side of Lower order solution
     62         *
     63         * LHS = [ML − theta*detlat *(K+D)^n+1]
     64         *
     65         */
     66
     67        /*Intermediaries*/
     68        int        dof,ncols,ncols2,rstart,rend;
     69        double     d,mi,dmax = 0.;
     70        Mat        LHS   = NULL;
     71        int*       cols  = NULL;
     72        int*       cols2 = NULL;
     73        double*    vals  = NULL;
     74        double*    vals2 = NULL;
     75
     76        MatDuplicate(K,MAT_SHARE_NONZERO_PATTERN,&LHS);
     77        MatGetOwnershipRange(K,&rstart,&rend);
     78        for(int row=rstart; row<rend; row++){
     79                MatGetRow(K,row,&ncols, (const int**)&cols, (const double**)&vals);
     80                MatGetRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     81                _assert_(ncols==ncols2);
     82                for(int j=0; j<ncols; j++) {
     83                        _assert_(cols[j]==cols2[j]);
     84                        d = -theta*deltat*(vals[j] + vals2[j]);
     85                        if(cols[j]==row){
     86                                VecGetValues(Ml,1,(const int*)&cols[j],&mi);
     87                                d += mi;
     88                        }
     89                        if(fabs(d)>dmax) dmax = fabs(d);
     90                        MatSetValue(LHS,row,cols[j],(const double)d,INSERT_VALUES);
     91                }
     92                MatRestoreRow(K,row,&ncols, (const int**)&cols, (const double**)&vals);
     93                MatRestoreRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     94        }
     95
     96        /*Penalize Dirichlet boundary*/
     97        dmax = dmax * 1.e+3;
     98        for(int i=0;i<femmodel->constraints->Size();i++){
     99                Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
     100                if(constraint->InAnalysis(configuration_type)){
     101                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
     102                        if(dof!=-1){
     103                                MatSetValue(LHS,dof,dof,(const double)dmax,INSERT_VALUES);
     104                        }
     105                }
     106        }
     107        MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
     108        MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
     109
     110        /*Clean up and assign output pointer*/
     111        *pdmax = dmax;
     112        *pLHS  = LHS;
     113}/*}}}*/
     114#endif
    10115void solutionsequence_fct(FemModel* femmodel){
    11116
    … …  
    18123        Vector<IssmDouble>*  ys = NULL;
    19124
    20         IssmDouble theta,deltat;
     125        IssmDouble theta,deltat,dmax=0.;
    21126        int        dof,ncols,ncols2,ncols3,rstart,rend;
    22127        int        configuration_type,analysis_type;
    23         double     d,diagD,mi;
    24         double     dmax = 0.;
     128        double     d,mi;
    25129        int*       cols  = NULL;
    26130        int*       cols2 = NULL;
    … …  
    45149        analysis->FctKMatrix(&K,NULL,femmodel);
    46150
    47         /*Create D Matrix*/
    48         #ifdef _HAVE_PETSC_
    49         Mat K_transp = NULL;
     151        /*Convert matrices to PETSc matrices*/
    50152        Mat D_petsc  = NULL;
    51153        Mat LHS      = NULL;
    … …  
    53155        Vec Ml_petsc = Ml->pvector->vector;
    54156        Mat Mc_petsc = Mc->pmatrix->matrix;
    55         MatTranspose(K_petsc,MAT_INITIAL_MATRIX,&K_transp);
    56         MatDuplicate(K_petsc,MAT_SHARE_NONZERO_PATTERN,&D_petsc);
    57         MatGetOwnershipRange(K_transp,&rstart,&rend);
    58         for(int row=rstart; row<rend; row++){
    59                 diagD = 0.;
    60                 MatGetRow(K_petsc ,row,&ncols, (const int**)&cols, (const double**)&vals);
    61                 MatGetRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    62                 _assert_(ncols==ncols2);
    63                 for(int j=0; j<ncols; j++) {
    64                         _assert_(cols[j]==cols2[j]);
    65                         d = max(max(-vals[j],-vals2[j]),0.);
    66                         MatSetValues(D_petsc,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
    67                         if(cols[j]!=row) diagD -= d;
    68                 }
    69                 MatSetValues(D_petsc,1,&row,1,&row,(const double*)&diagD,INSERT_VALUES);
    70                 MatRestoreRow(K_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
    71                 MatRestoreRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    72         }
    73         MatAssemblyBegin(D_petsc,MAT_FINAL_ASSEMBLY);
    74         MatAssemblyEnd(  D_petsc,MAT_FINAL_ASSEMBLY);
    75         MatFree(&K_transp);
     157
     158        /*Create D Matrix*/
     159        #ifdef _HAVE_PETSC_
     160        CreateDMatrix(&D_petsc,K_petsc);
    76161
    77162        /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
    78163        MatDuplicate(K_petsc,MAT_SHARE_NONZERO_PATTERN,&LHS);
     164        MatGetOwnershipRange(K_petsc,&rstart,&rend);
    79165        for(int row=rstart; row<rend; row++){
    80166                MatGetRow(K_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
Note: See TracChangeset for help on using the changeset viewer.