Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24135) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24136) @@ -567,4 +567,5 @@ IssmDouble h,hx,hy,hz,vx,vy,vz; IssmDouble tau_parameter,diameter; + IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver; IssmDouble D_scalar; IssmDouble* xyz_list = NULL; @@ -595,5 +596,4 @@ Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input); Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input); - if(stabilization==2) diameter=element->MinEdgeLength(xyz_list); /*Enthalpy diffusion parameter*/ @@ -660,11 +660,14 @@ &Ke->values[0],1); } + /*SUPG*/ else if(stabilization==2){ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + diameter=element->MinEdgeLength(xyz_list); tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice); for(int i=0;ivalues[i*numnodes+j]+=tau_parameter*D_scalar* - ((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]); + ((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]); } } @@ -678,4 +681,19 @@ } } + /*anisotropic SUPG*/ + else if(stabilization==3){ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->ElementSizes(&hx,&hy,&hz); + element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u-um,v-vm,w-wm,hx,hy,hz,kappa/rho_ice); + tau_parameter_hor=tau_parameter_anisotropic[0]; + tau_parameter_ver=tau_parameter_anisotropic[1]; + for(int i=0;ivalues[i*numnodes+j]+=D_scalar* + (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])* + (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]); + } + } + } } @@ -769,5 +787,6 @@ IssmDouble pressure, d1pressure[3], d2pressure; IssmDouble waterfractionpicard; - IssmDouble kappa,tau_parameter,diameter,kappa_w; + IssmDouble kappa,tau_parameter,diameter,hx,hy,hz,kappa_w; + IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver; IssmDouble u,v,w; IssmDouble scalar_def, scalar_sens ,scalar_transient; @@ -801,8 +820,4 @@ Input* enthalpy_input=NULL; if(reCast(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);} - if(stabilization==2){ - diameter=element->MinEdgeLength(xyz_list); - kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.); - } /* Start looping on the number of gaussian points: */ @@ -848,8 +863,10 @@ for(i=0;ivalues[i]+=scalar_transient*basis[i]; } - + + /* SUPG */ if(stabilization==2){ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - + diameter=element->MinEdgeLength(xyz_list); + kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.); vx_input->GetInputValue(&u,gauss); vy_input->GetInputValue(&v,gauss); @@ -862,4 +879,18 @@ for(i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]); } + } + /* anisotropic SUPG */ + else if(stabilization==3){ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->ElementSizes(&hx,&hy,&hz); + kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.); + vx_input->GetInputValue(&u,gauss); + vy_input->GetInputValue(&v,gauss); + vz_input->GetInputValue(&w,gauss); + element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa/rho_ice); + tau_parameter_hor=tau_parameter_anisotropic[0]; + tau_parameter_ver=tau_parameter_anisotropic[1]; + + for(i=0;ivalues[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]); } } Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 24135) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 24136) @@ -378,4 +378,5 @@ IssmDouble h,hx,hy,hz,vx,vy,vz,D_scalar; IssmDouble tau_parameter,diameter; + IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver; IssmDouble* xyz_list = NULL; @@ -406,5 +407,4 @@ Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input); Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input); - if(stabilization==2) diameter=element->MinEdgeLength(xyz_list); /* Start looping on the number of gaussian points: */ @@ -469,9 +469,11 @@ } else if(stabilization==2){ + diameter=element->MinEdgeLength(xyz_list); tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); for(int i=0;ivalues[i*numnodes+j]+=tau_parameter*D_scalar* - ((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]); + ((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]); } } @@ -482,4 +484,19 @@ 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]); } + } + } + } + /*anisotropic SUPG*/ + else if(stabilization==3){ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->ElementSizes(&hx,&hy,&hz); + element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u-um,v-vm,w-wm,hx,hy,hz,kappa); + tau_parameter_hor=tau_parameter_anisotropic[0]; + tau_parameter_ver=tau_parameter_anisotropic[1]; + for(int i=0;ivalues[i*numnodes+j]+=D_scalar* + (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])* + (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]); } } @@ -635,5 +652,6 @@ IssmDouble Jdet,phi,dt; IssmDouble temperature; - IssmDouble tau_parameter,diameter; + IssmDouble tau_parameter,diameter,hx,hy,hz; + IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver; IssmDouble u,v,w; IssmDouble scalar_def,scalar_transient; @@ -661,5 +679,4 @@ Input* temperature_input = NULL; if(reCast(dt)){temperature_input = element->GetInput(TemperatureEnum); _assert_(temperature_input);} - if(stabilization==2) diameter=element->MinEdgeLength(xyz_list); /* Start looping on the number of gaussian points: */ @@ -686,5 +703,5 @@ if(stabilization==2){ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - + diameter=element->MinEdgeLength(xyz_list); vx_input->GetInputValue(&u,gauss); vy_input->GetInputValue(&v,gauss); @@ -697,4 +714,17 @@ for(int i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]); } + } + /* anisotropic SUPG */ + else if(stabilization==3){ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->ElementSizes(&hx,&hy,&hz); + vx_input->GetInputValue(&u,gauss); + vy_input->GetInputValue(&v,gauss); + vz_input->GetInputValue(&w,gauss); + element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa); + tau_parameter_hor=tau_parameter_anisotropic[0]; + tau_parameter_ver=tau_parameter_anisotropic[1]; + + for(int i=0;ivalues[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]); } } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24135) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24136) @@ -319,4 +319,5 @@ virtual Element* SpawnTopElement(void)=0; virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0; + virtual IssmDouble StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0; virtual void StrainRateparallel(void)=0; virtual void StrainRateperpendicular(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24135) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24136) @@ -3142,5 +3142,5 @@ IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){/*{{{*/ /*Compute stabilization parameter*/ - /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/ + /*kappa=thermalconductivity/(rho_ice*heatcapacity) for thermal model*/ /*kappa=enthalpydiffusionparameter for enthalpy model*/ @@ -3155,4 +3155,28 @@ return tau_parameter; +} +IssmDouble Penta::StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){/*{{{*/ + /*Compute stabilization parameter*/ + /*kappa=thermalconductivity/(rho_ice*heatcapacity) for thermal model*/ + /*kappa=enthalpydiffusionparameter for enthalpy model*/ + + IssmDouble normu,hk,C,area,p; + + /* compute tau for the horizontal direction */ + p=2.; C=3.; + normu=pow(pow(u,p)+pow(v,p)+pow(w,p),1./p); + hk=sqrt(pow(hx,2)+pow(hy,2)); + + if(normu*hk/(C*2*kappa)<1){ + tau_parameter_anisotropic[0]=pow(hk,2)/(C*2*2*kappa); + } + else tau_parameter_anisotropic[0]=hk/(2*normu); + + /* compute tau for the vertical direction */ + hk=hz; + if(normu*hk/(C*2*kappa)<1){ + tau_parameter_anisotropic[1]=pow(hk,2)/(C*2*2*kappa); + } + else tau_parameter_anisotropic[1]=hk/(2*normu); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24135) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24136) @@ -162,4 +162,5 @@ Tria* SpawnTria(int index1,int index2,int index3); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); + IssmDouble StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa); void StressIntensityFactor(); void StrainRateparallel(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24135) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24136) @@ -138,4 +138,5 @@ Element* SpawnTopElement(void){_error_("not implemented yet");}; IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; + IssmDouble StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");}; void StrainRateparallel(void){_error_("not implemented yet");}; void StrainRateperpendicular(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24135) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24136) @@ -147,4 +147,5 @@ Tria* SpawnTria(int index1,int index2,int index3); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; + IssmDouble StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");}; void StrainRateparallel(void){_error_("not implemented yet");}; void StrainRateperpendicular(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24135) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24136) @@ -215,4 +215,5 @@ Seg* SpawnSeg(int index1,int index2); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; + IssmDouble StabilizationParameterAnisotropic(IssmDouble* tau_parameter_ansiotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");}; void UpdateConstraintsExtrudeFromBase(void); void UpdateConstraintsExtrudeFromTop(void); Index: /issm/trunk-jpl/src/m/classes/thermal.js =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.js (revision 24135) +++ /issm/trunk-jpl/src/m/classes/thermal.js (revision 24136) @@ -41,5 +41,5 @@ fielddisplay(this,'spctemperature','temperature constraints (NaN means no constraint) [K]'); - fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'); + fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG'); fielddisplay(this,'reltol','relative tolerance convergence criterion for enthalpy'); fielddisplay(this,'maxiter','maximum number of non linear iterations'); @@ -90,5 +90,5 @@ if(!ArrayAnyEqual(ArrayIsMember('ThermalAnalysis',analyses),1) & !ArrayAnyEqual(ArrayIsMember('EnthalpyAnalysis',analyses),1) | (solution == 'TransientSolution' & md.trans.isthermal==0)) return; - checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 ,1, 2]); + checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 ,1, 2, 3]); checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1); checkfield(md,'fieldname','thermal.fe','values',['P1','P1xP2','P1xP3']); Index: /issm/trunk-jpl/src/m/classes/thermal.m =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.m (revision 24135) +++ /issm/trunk-jpl/src/m/classes/thermal.m (revision 24136) @@ -78,5 +78,5 @@ if (~ismember('ThermalAnalysis',analyses) & ~ismember('EnthalpyAnalysis',analyses)) | (strcmp(solution,'TransientSolution') & md.transient.isthermal==0), return; end - md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2]); + md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2 3]); md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1,'>=',0); md = checkfield(md,'fieldname','thermal.fe','values',{'P1','P1xP2','P1xP3'}); @@ -107,5 +107,5 @@ fielddisplay(self,'spctemperature','temperature constraints (NaN means no constraint) [K]'); - fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'); + fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG'); fielddisplay(self,'reltol','relative tolerance convergence criterion for enthalpy'); fielddisplay(self,'maxiter','maximum number of non linear iterations'); Index: /issm/trunk-jpl/src/m/classes/thermal.py =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.py (revision 24135) +++ /issm/trunk-jpl/src/m/classes/thermal.py (revision 24136) @@ -33,5 +33,5 @@ string=' Thermal solution parameters:' string="%s\n%s"%(string,fielddisplay(self,'spctemperature','temperature constraints (NaN means no constraint) [K]')) - string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG')) + string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG')) string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of non linear iterations')) string="%s\n%s"%(string,fielddisplay(self,'reltol','relative tolerance criterion')) @@ -96,5 +96,5 @@ return md - md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2]) + md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2,3]) md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1) md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1) Index: /issm/trunk-jpl/test/NightlyRun/test331.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test331.m (revision 24136) +++ /issm/trunk-jpl/test/NightlyRun/test331.m (revision 24136) @@ -0,0 +1,29 @@ +%Test Name: SquareSheetConstrainedAnisotropicSUPG +md=triangle(model(),'../Exp/Square.exp',180000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/SquareSheetConstrained.par'); +md=extrude(md,3,1.); +md=setflowequation(md,'SSA','all'); +md.timestepping.time_step=0.; +md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1); +md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1); +md.thermal.isenthalpy=0; +md.thermal.isdynamicbasalspc=0; +md.thermal.stabilization=3; + +md.cluster=generic('name',oshostname(),'np',3); + +field_names={}; +field_tolerances={}; +field_values={}; + +for i=[0 1] +md.thermal.isenthalpy=i; + disp(' '); + disp(['====== Testing Thermal model with anisotropic SUPG and isenthalpy=',num2str(i),' =====']); + md=solve(md,'Thermal'); + %Fields and tolerances to track changes + field_names ={field_names{:},['Temperature' i]}; + field_tolerances={field_tolerances{:},1e-13}; + field_values={field_values{:},(md.results.ThermalSolution.Temperature)}; +end