Changeset 24136


Ignore:
Timestamp:
09/09/19 02:54:07 (6 years ago)
Author:
rueckamp
Message:

NEW: anisotropic SUPG for Enthalpy/Thermal Solver

Location:
issm/trunk-jpl
Files:
2 added
11 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

    r24060 r24136  
    567567        IssmDouble  h,hx,hy,hz,vx,vy,vz;
    568568        IssmDouble  tau_parameter,diameter;
     569        IssmDouble  tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;   
    569570        IssmDouble  D_scalar;
    570571        IssmDouble* xyz_list = NULL;
    … …  
    595596        Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
    596597        Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
    597         if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
    598598
    599599        /*Enthalpy diffusion parameter*/
    … …  
    660660                                                &Ke->values[0],1);
    661661                }
     662                /*SUPG*/
    662663                else if(stabilization==2){
    663664                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     665                        diameter=element->MinEdgeLength(xyz_list);                     
    664666                        tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
    665667                        for(int i=0;i<numnodes;i++){
    666668                                for(int j=0;j<numnodes;j++){
    667669                                        Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
    668                                           ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
     670                                          ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*
     671                                          ((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
    669672                                }
    670673                        }
    … …  
    678681                        }
    679682                }
     683                /*anisotropic SUPG*/
     684                else if(stabilization==3){
     685                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     686                        element->ElementSizes(&hx,&hy,&hz);                     
     687                        element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u-um,v-vm,w-wm,hx,hy,hz,kappa/rho_ice);
     688                        tau_parameter_hor=tau_parameter_anisotropic[0];
     689                        tau_parameter_ver=tau_parameter_anisotropic[1];
     690                        for(int i=0;i<numnodes;i++){
     691                                for(int j=0;j<numnodes;j++){
     692                                        Ke->values[i*numnodes+j]+=D_scalar*
     693                                                (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+i]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+i]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+i])*
     694                                                (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+j]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+j]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+j]);
     695                                }
     696                        }
     697                }
    680698        }
    681699
    … …  
    769787        IssmDouble  pressure, d1pressure[3], d2pressure;
    770788        IssmDouble  waterfractionpicard;
    771         IssmDouble  kappa,tau_parameter,diameter,kappa_w;
     789        IssmDouble  kappa,tau_parameter,diameter,hx,hy,hz,kappa_w;
     790        IssmDouble  tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
    772791        IssmDouble  u,v,w;
    773792        IssmDouble  scalar_def, scalar_sens ,scalar_transient;
    … …  
    801820        Input* enthalpy_input=NULL;
    802821        if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
    803         if(stabilization==2){
    804                 diameter=element->MinEdgeLength(xyz_list);
    805                 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
    806         }
    807822
    808823        /* Start  looping on the number of gaussian points: */
    … …  
    848863                        for(i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
    849864                }
    850 
     865               
     866                /* SUPG */
    851867                if(stabilization==2){
    852868                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    853 
     869                        diameter=element->MinEdgeLength(xyz_list);
     870                        kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
    854871                        vx_input->GetInputValue(&u,gauss);
    855872                        vy_input->GetInputValue(&v,gauss);
    … …  
    862879                                for(i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
    863880                        }
     881                }
     882                /* anisotropic SUPG */
     883                else if(stabilization==3){
     884                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     885                        element->ElementSizes(&hx,&hy,&hz);
     886                        kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
     887                        vx_input->GetInputValue(&u,gauss);
     888                        vy_input->GetInputValue(&v,gauss);
     889                        vz_input->GetInputValue(&w,gauss);
     890                        element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa/rho_ice);
     891                        tau_parameter_hor=tau_parameter_anisotropic[0];
     892                        tau_parameter_ver=tau_parameter_anisotropic[1];
     893                       
     894                        for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*(tau_parameter_hor*u*dbasis[0*numnodes+i]+tau_parameter_hor*v*dbasis[1*numnodes+i]+tau_parameter_ver*w*dbasis[2*numnodes+i]);
    864895                }
    865896        }

  • TabularUnified issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp

    r24093 r24136  
    378378        IssmDouble  h,hx,hy,hz,vx,vy,vz,D_scalar;
    379379        IssmDouble  tau_parameter,diameter;
     380        IssmDouble  tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;   
    380381        IssmDouble* xyz_list = NULL;
    381382
    … …  
    406407        Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
    407408        Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
    408         if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
    409409
    410410        /* Start  looping on the number of gaussian points: */
    … …  
    469469                }
    470470                else if(stabilization==2){
     471                        diameter=element->MinEdgeLength(xyz_list);
    471472                        tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
    472473                        for(int i=0;i<numnodes;i++){
    473474                                for(int j=0;j<numnodes;j++){
    474475                                        Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
    475                                           ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
     476                                          ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*
     477                                          ((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
    476478                                }
    477479                        }
    … …  
    482484                                                Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
    483485                                        }
     486                                }
     487                        }
     488                }
     489                /*anisotropic SUPG*/
     490                else if(stabilization==3){
     491                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     492                        element->ElementSizes(&hx,&hy,&hz);
     493                        element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u-um,v-vm,w-wm,hx,hy,hz,kappa);
     494                        tau_parameter_hor=tau_parameter_anisotropic[0];
     495                        tau_parameter_ver=tau_parameter_anisotropic[1];
     496                        for(int i=0;i<numnodes;i++){
     497                                for(int j=0;j<numnodes;j++){
     498                                        Ke->values[i*numnodes+j]+=D_scalar*
     499                                                (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+i]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+i]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+i])*
     500                                                (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+j]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+j]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+j]);
    484501                                }
    485502                        }
    … …  
    635652        IssmDouble  Jdet,phi,dt;
    636653        IssmDouble  temperature;
    637         IssmDouble  tau_parameter,diameter;
     654        IssmDouble  tau_parameter,diameter,hx,hy,hz;
     655        IssmDouble  tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
    638656        IssmDouble  u,v,w;
    639657        IssmDouble  scalar_def,scalar_transient;
    … …  
    661679        Input* temperature_input = NULL;
    662680        if(reCast<bool,IssmDouble>(dt)){temperature_input = element->GetInput(TemperatureEnum); _assert_(temperature_input);}
    663         if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
    664681
    665682        /* Start  looping on the number of gaussian points: */
    … …  
    686703                if(stabilization==2){
    687704                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    688 
     705                        diameter=element->MinEdgeLength(xyz_list);
    689706                        vx_input->GetInputValue(&u,gauss);
    690707                        vy_input->GetInputValue(&v,gauss);
    … …  
    697714                                for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
    698715                        }
     716                }
     717                /* anisotropic SUPG */
     718                else if(stabilization==3){
     719                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     720                        element->ElementSizes(&hx,&hy,&hz);
     721                        vx_input->GetInputValue(&u,gauss);
     722                        vy_input->GetInputValue(&v,gauss);
     723                        vz_input->GetInputValue(&w,gauss);
     724                        element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa);
     725                        tau_parameter_hor=tau_parameter_anisotropic[0];
     726                        tau_parameter_ver=tau_parameter_anisotropic[1];
     727
     728                        for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*(tau_parameter_hor*u*dbasis[0*numnodes+i]+tau_parameter_hor*v*dbasis[1*numnodes+i]+tau_parameter_ver*w*dbasis[2*numnodes+i]);
    699729                }
    700730        }

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

    r24089 r24136  
    319319                virtual Element*   SpawnTopElement(void)=0;
    320320                virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
     321                virtual IssmDouble StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0;         
    321322                virtual void        StrainRateparallel(void)=0;
    322323                virtual void        StrainRateperpendicular(void)=0;

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

    r24119 r24136  
    31423142IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){/*{{{*/
    31433143        /*Compute stabilization parameter*/
    3144         /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
     3144        /*kappa=thermalconductivity/(rho_ice*heatcapacity) for thermal model*/
    31453145        /*kappa=enthalpydiffusionparameter for enthalpy model*/
    31463146
    … …  
    31553155
    31563156        return tau_parameter;
     3157}
     3158IssmDouble Penta::StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){/*{{{*/
     3159        /*Compute stabilization parameter*/
     3160        /*kappa=thermalconductivity/(rho_ice*heatcapacity) for thermal model*/
     3161        /*kappa=enthalpydiffusionparameter for enthalpy model*/
     3162
     3163        IssmDouble normu,hk,C,area,p;
     3164
     3165        /* compute tau for the horizontal direction */
     3166        p=2.; C=3.;
     3167        normu=pow(pow(u,p)+pow(v,p)+pow(w,p),1./p);
     3168        hk=sqrt(pow(hx,2)+pow(hy,2));
     3169       
     3170        if(normu*hk/(C*2*kappa)<1){
     3171                tau_parameter_anisotropic[0]=pow(hk,2)/(C*2*2*kappa);
     3172        }
     3173        else tau_parameter_anisotropic[0]=hk/(2*normu);
     3174       
     3175        /* compute tau for the vertical direction */
     3176        hk=hz;
     3177        if(normu*hk/(C*2*kappa)<1){
     3178                tau_parameter_anisotropic[1]=pow(hk,2)/(C*2*2*kappa);
     3179        }
     3180        else tau_parameter_anisotropic[1]=hk/(2*normu);
    31573181}
    31583182/*}}}*/

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

    r24119 r24136  
    162162                Tria*            SpawnTria(int index1,int index2,int index3);
    163163                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
     164                IssmDouble     StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa);
    164165                void           StressIntensityFactor();
    165166                void           StrainRateparallel();

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

    r24089 r24136  
    138138                Element*    SpawnTopElement(void){_error_("not implemented yet");};
    139139                IssmDouble  StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
     140                IssmDouble  StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy,  IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");};
    140141                void        StrainRateparallel(void){_error_("not implemented yet");};
    141142                void        StrainRateperpendicular(void){_error_("not implemented yet");};

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

    r24089 r24136  
    147147                Tria*       SpawnTria(int index1,int index2,int index3);
    148148                IssmDouble  StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
     149                IssmDouble  StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");};
    149150                void        StrainRateparallel(void){_error_("not implemented yet");};
    150151                void        StrainRateperpendicular(void){_error_("not implemented yet");};

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

    r24089 r24136  
    215215                Seg*             SpawnSeg(int index1,int index2);
    216216                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
     217                IssmDouble     StabilizationParameterAnisotropic(IssmDouble* tau_parameter_ansiotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");};
    217218                void           UpdateConstraintsExtrudeFromBase(void);
    218219                void           UpdateConstraintsExtrudeFromTop(void);

  • TabularUnified issm/trunk-jpl/src/m/classes/thermal.js

    r22677 r24136  
    4141
    4242                fielddisplay(this,'spctemperature','temperature constraints (NaN means no constraint) [K]');
    43                 fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
     43                fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG');
    4444                fielddisplay(this,'reltol','relative tolerance convergence criterion for enthalpy');
    4545                fielddisplay(this,'maxiter','maximum number of non linear iterations');
    … …  
    9090                if(!ArrayAnyEqual(ArrayIsMember('ThermalAnalysis',analyses),1) & !ArrayAnyEqual(ArrayIsMember('EnthalpyAnalysis',analyses),1)  | (solution == 'TransientSolution' & md.trans.isthermal==0)) return;
    9191
    92                 checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 ,1, 2]);
     92                checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 ,1, 2, 3]);
    9393                checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1);
    9494                checkfield(md,'fieldname','thermal.fe','values',['P1','P1xP2','P1xP3']);

  • TabularUnified issm/trunk-jpl/src/m/classes/thermal.m

    r21670 r24136  
    7878                        if (~ismember('ThermalAnalysis',analyses) & ~ismember('EnthalpyAnalysis',analyses)) | (strcmp(solution,'TransientSolution') & md.transient.isthermal==0), return; end
    7979
    80                         md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2]);
     80                        md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2 3]);
    8181                        md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1,'>=',0);
    8282                        md = checkfield(md,'fieldname','thermal.fe','values',{'P1','P1xP2','P1xP3'});
    … …  
    107107
    108108                        fielddisplay(self,'spctemperature','temperature constraints (NaN means no constraint) [K]');
    109                         fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
     109                        fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG');
    110110                        fielddisplay(self,'reltol','relative tolerance convergence criterion for enthalpy');
    111111                        fielddisplay(self,'maxiter','maximum number of non linear iterations');

  • TabularUnified issm/trunk-jpl/src/m/classes/thermal.py

    r23170 r24136  
    3333                string='   Thermal solution parameters:'
    3434                string="%s\n%s"%(string,fielddisplay(self,'spctemperature','temperature constraints (NaN means no constraint) [K]'))
    35                 string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'))
     35                string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG'))
    3636                string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of non linear iterations'))
    3737                string="%s\n%s"%(string,fielddisplay(self,'reltol','relative tolerance criterion'))
    … …  
    9696                        return md
    9797
    98                 md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2])
     98                md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2,3])
    9999                md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1)
    100100                md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1)
Note: See TracChangeset for help on using the changeset viewer.