Changeset 25942


Ignore:
Timestamp:
01/18/21 10:20:53 (4 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixing Inwoo's commit

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

    r25941 r25942  
    55#include "../modules/modules.h"
    66#include "../solutionsequences/solutionsequences.h"
     7
     8//#define INWOOVZ 1 //For Inwoo Park 0: incompressible assumption(IA), 1: internal deformation(ID), 2: IA+ID
    79
    810/*Model processing*/
    … …  
    112114        }
    113115        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
     116        #ifdef INWOOVZ
    114117        iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
     118        #endif
    115119
    116120
    … …  
    154158
    155159        /* specific parameters*/
    156         parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvzsolution",StressbalanceIsVzSolutionEnum));
    157         parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvelbc",StressbalanceIsVelBCEnum));
    158160
    159161}/*}}}*/
    … …  
    518520        IssmDouble   rho_ice,g;
    519521
    520         int          isvzsolution;
    521         IssmDouble   isvelbc, rheology_n, eta;
    522 
    523522        /*Get the approximation and do nothing if the element in FS or None*/
    524523        element->GetInputValue(&approximation,ApproximationEnum);
    … …  
    526525                return;
    527526        }
    528 
    529         /*Get parameters*/
    530         element->FindParam(&isvzsolution,StressbalanceIsVzSolutionEnum);
    531         element->FindParam(&isvelbc,StressbalanceIsVelBCEnum);
    532527
    533528        /*Get dof list and vertices coordinates: */
    … …  
    544539        IssmDouble*  pressure  = xNew<IssmDouble>(numnodes);
    545540        IssmDouble*  surface   = xNew<IssmDouble>(numnodes);
    546         IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
    547         IssmDouble*  base            = xNew<IssmDouble>(numnodes);
    548         IssmDouble*  smb       = xNew<IssmDouble>(numnodes);
    549541
    550542        /*Use the dof list to index into the solution vector vz: */
    … …  
    588580        }
    589581
     582
     583        #ifdef INWOOVZ
     584        IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
     585        IssmDouble*  base            = xNew<IssmDouble>(numnodes);
     586        IssmDouble*  smb       = xNew<IssmDouble>(numnodes);
     587
     588        IssmDouble rheology_n=element->material->GetN();
     589        IssmDouble isvelbc   = 10./(365.25*24*3600); /*10 m/yr*/
    590590        /*Set analytical vertical velocity field at slow flow region.*/
    591         if (isvzsolution == 1){
    592                 rheology_n=element->material->GetN();
     591        if(INWOOVZ == 1){
    593592                element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
    594593                element->GetInputListOnNodes(&base[0],BaseEnum,0.);
    595594                element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
    596595                for(i=0;i<numnodes;i++){
    597                         eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
     596                        IssmDouble eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
    598597                        vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
    599598                }       
    600599        }
    601         else if (isvzsolution == 2) {
    602                 rheology_n=element->material->GetN();
     600        else if(INWOOVZ== 2) {
    603601                element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
    604602                element->GetInputListOnNodes(&base[0],BaseEnum,0.);
    … …  
    606604                for(i=0;i<numnodes;i++){
    607605                        if ( vx[i]*vx[i]+vy[i]*vy[i] < isvelbc*isvelbc ){
    608                                 eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
     606                                IssmDouble eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
    609607                                vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
    610608                        }
    611609                }       
    612610        }
     611        else if(INWOOVZ==0){
     612                /*nothing to add*/
     613        }
     614        else{
     615                _error_("not supported yet");
     616        }
     617
     618        /*Cleanup*/
     619        xDelete<IssmDouble>(smb);
     620        xDelete<IssmDouble>(base);
     621        xDelete<IssmDouble>(thickness);
     622        #endif
    613623
    614624
    … …  
    637647
    638648        /*Free ressources:*/
    639         xDelete<IssmDouble>(smb);
    640         xDelete<IssmDouble>(base);
    641         xDelete<IssmDouble>(thickness);
    642649        xDelete<IssmDouble>(surface);
    643650        xDelete<IssmDouble>(pressure);

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

    r25934 r25942  
    30153015        IssmDouble* values    = xNew<IssmDouble>(numnodes);
    30163016
     3017                if(this->Id()==104 || this->Id()==266 || this->Id()==264 || this->Id()==247 || this->Id()==263 || this->Id()==249 || this->Id()==250){
     3018                        if(enum_type==MaskIceLevelsetEnum) printf("-----------------\n");
     3019                }
     3020
    30173021        /*Use the dof list to index into the solution vector: */
    30183022        for(int i=0;i<numnodes;i++){
    30193023                values[i]=solution[doflist[i]];
     3024                if(this->Id()==104 || this->Id()==266 || this->Id()==264 || this->Id()==247 || this->Id()==263 || this->Id()==249 || this->Id()==250){
     3025                        if(enum_type==MaskIceLevelsetEnum) printf("%g\n",values[i]);
     3026                        if(enum_type==MaskIceLevelsetEnum && values[i]>0.5 && values[i]<1.5 )this->nodes[i]->DeepEcho();
     3027                }
    30203028                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    30213029                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector, SID = " << this->Sid());
    30223030        }
     3031
    30233032
    30243033        /*Add input to the element: */

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r25839 r25942  
    13991399syn keyword cType Cfsurfacesquare
    14001400syn keyword cType Channel
     1401syn keyword cType classes
    14011402syn keyword cType Constraint
    14021403syn keyword cType Constraints
    … …  
    14051406syn keyword cType ControlInput
    14061407syn keyword cType Covertree
     1408syn keyword cType DatasetInput
    14071409syn keyword cType DataSetParam
    1408 syn keyword cType DatasetInput
    14091410syn keyword cType Definition
    14101411syn keyword cType DependentObject
    … …  
    14191420syn keyword cType ElementInput
    14201421syn keyword cType ElementMatrix
     1422syn keyword cType Elements
    14211423syn keyword cType ElementVector
    1422 syn keyword cType Elements
    14231424syn keyword cType ExponentialVariogram
    14241425syn keyword cType ExternalResult
    … …  
    14271428syn keyword cType Friction
    14281429syn keyword cType Gauss
     1430syn keyword cType GaussianVariogram
     1431syn keyword cType gaussobjects
    14291432syn keyword cType GaussPenta
    14301433syn keyword cType GaussSeg
    14311434syn keyword cType GaussTetra
    14321435syn keyword cType GaussTria
    1433 syn keyword cType GaussianVariogram
    14341436syn keyword cType GenericExternalResult
    14351437syn keyword cType GenericOption
    … …  
    14461448syn keyword cType IssmDirectApplicInterface
    14471449syn keyword cType IssmParallelDirectApplicInterface
     1450syn keyword cType krigingobjects
    14481451syn keyword cType Load
    14491452syn keyword cType Loads
    … …  
    14561459syn keyword cType Matice
    14571460syn keyword cType Matlitho
     1461syn keyword cType matrixobjects
    14581462syn keyword cType MatrixParam
    14591463syn keyword cType Misfit
    … …  
    14681472syn keyword cType Observations
    14691473syn keyword cType Option
     1474syn keyword cType Options
    14701475syn keyword cType OptionUtilities
    1471 syn keyword cType Options
    14721476syn keyword cType Param
    14731477syn keyword cType Parameters
    … …  
    14831487syn keyword cType Regionaloutput
    14841488syn keyword cType Results
     1489syn keyword cType Riftfront
    14851490syn keyword cType RiftStruct
    1486 syn keyword cType Riftfront
    14871491syn keyword cType SealevelMasks
    14881492syn keyword cType Seg
    14891493syn keyword cType SegInput
     1494syn keyword cType Segment
    14901495syn keyword cType SegRef
    1491 syn keyword cType Segment
    14921496syn keyword cType SpcDynamic
    14931497syn keyword cType SpcStatic
    … …  
    15081512syn keyword cType Vertex
    15091513syn keyword cType Vertices
    1510 syn keyword cType classes
    1511 syn keyword cType gaussobjects
    1512 syn keyword cType krigingobjects
    1513 syn keyword cType matrixobjects
    15141514syn keyword cType AdjointBalancethickness2Analysis
    15151515syn keyword cType AdjointBalancethicknessAnalysis
    … …  
    15301530syn keyword cType FreeSurfaceBaseAnalysis
    15311531syn keyword cType FreeSurfaceTopAnalysis
     1532syn keyword cType GiaAnalysis
    15321533syn keyword cType GLheightadvectionAnalysis
    1533 syn keyword cType GiaAnalysis
    15341534syn keyword cType HydrologyDCEfficientAnalysis
    15351535syn keyword cType HydrologyDCInefficientAnalysis

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r25941 r25942  
    441441        StressbalanceRiftPenaltyThresholdEnum,
    442442        StressbalanceShelfDampeningEnum,
    443         StressbalanceIsVzSolutionEnum,
    444         StressbalanceIsVelBCEnum,
    445443        ThermalIsdrainicecolumnEnum,
    446444        ThermalIsdynamicbasalspcEnum,

  • issm/trunk-jpl/src/m/classes/stressbalance.m

    r25941 r25942  
    2323                loadingforce           = NaN;
    2424                requested_outputs      = {};
    25                 isvzsolution           = NaN;
    26                 isvelbc                = NaN;   
    2725        end
    2826        methods
    … …  
    7775                         self.rift_penalty_lock=10;
    7876
    79                          %Set vertical velocity option
    80                          self.isvzsolution = 0;
    81                          self.isvelbc      = 10.;
    82 
    8377                         %output default:
    8478                         self.requested_outputs={'default'};
    … …  
    10296                        md = checkfield(md,'fieldname','stressbalance.referential','size',[md.mesh.numberofvertices 6]);
    10397                        md = checkfield(md,'fieldname','stressbalance.loadingforce','size',[md.mesh.numberofvertices 3]);
    104                         md = checkfield(md,'fieldname','stressbalance.isvzsolution','numel',[1],'values',[0 1 2]);
    105                         md = checkfield(md,'fieldname','stressbalance.isvelbc','numel',[1],'>=',0.0);
    10698                        md = checkfield(md,'fieldname','stressbalance.requested_outputs','stringrow',1);
    10799                        if ~any(isnan(md.stressbalance.vertex_pairing)),
    … …  
    166158                        fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
    167159                        fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized');
    168 
    169                         disp(sprintf('\n      %s','Estimation of vertical velocity field options:'));
    170                         fielddisplay(self,'isvzsolution','0: incompressible assumption(IA), 1: internal deformation(ID), 2: IA+ID');
    171                         fielddisplay(self,'isvelbc','If turned on isvzsolution as 2, region (vel < isvelbc) calculates for ID.');
    172160
    173161                        disp(sprintf('\n      %s','Other:'));
    … …  
    199187                        WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer');
    200188                        WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1);
    201                         WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvzsolution','format','Integer','numel',1,'values',[0 1 2]);
    202                         WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvelbc','format','Double','numel',1,'scale',1./yts);
    203189
    204190                        if size(self.loadingforce,2)==3,
Note: See TracChangeset for help on using the changeset viewer.