Skip to content

Commit

Permalink
CHG: adding factors in analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
hseroussi committed Jun 25, 2024
1 parent 44136bd commit 3a66e7c
Show file tree
Hide file tree
Showing 14 changed files with 114 additions and 72 deletions.
11 changes: 6 additions & 5 deletions src/c/analyses/EnthalpyAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -489,7 +489,7 @@ ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
int stabilization;
IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
IssmDouble h,hx,hy,hz,vx,vy,vz;
IssmDouble tau_parameter,diameter;
IssmDouble tau_parameter,diameter,factor;
IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
IssmDouble D_scalar;
IssmDouble* xyz_list = NULL;
Expand Down Expand Up @@ -534,11 +534,11 @@ ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
if(dt!=0.) D_scalar=D_scalar*dt;

/*Conduction: */
factor = D_scalar*kappa/rho_ice;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += D_scalar*kappa/rho_ice*(
dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
);
Ke->values[i*numnodes+j] += factor*(
dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]);
}
}

Expand Down Expand Up @@ -599,9 +599,10 @@ ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
}
if(dt!=0.){
D_scalar=gauss->weight*Jdet;
factor = tau_parameter*D_scalar;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
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]);
Ke->values[i*numnodes+j]+=factor*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/c/analyses/ExtrudeFromBaseAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/

for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j]);
Ke->values[i*numnodes+j] += D*(dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j]);
}
}
}
Expand Down
3 changes: 2 additions & 1 deletion src/c/analyses/ExtrudeFromTopAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -77,9 +77,10 @@ ElementMatrix* ExtrudeFromTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
D=gauss->weight*Jdet;
element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);

IssmDouble factor = gauss->weight*Jdet;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j]);
Ke->values[i*numnodes+j] += factor*(dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j]);
}
}
}
Expand Down
17 changes: 10 additions & 7 deletions src/c/analyses/FreeSurfaceBaseAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -324,18 +324,18 @@ ElementMatrix* FreeSurfaceBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
}
}
else if(stabilization==5){
D_scalar=gauss->weight*Jdet*dt;
D_scalar=gauss->weight*Jdet*dt*tau;
if(dim==2){
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j]+=tau*D_scalar*
Ke->values[i*numnodes+j]+=D_scalar*
(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*
(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
}
}
}
else{
for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
}
}
}
Expand All @@ -353,7 +353,7 @@ ElementVector* FreeSurfaceBaseAnalysis::CreatePVector(Element* element){/*{{{*/

/*Intermediaries*/
int domaintype,dim,stabilization;
IssmDouble Jdet,dt,intrusiondist;
IssmDouble Jdet,dt,intrusiondist,factor;
IssmDouble gmb,fmb,mb,bed,vx,vy,vz,tau,gldistance;
Element* basalelement = NULL;
IssmDouble *xyz_list = NULL;
Expand Down Expand Up @@ -485,21 +485,24 @@ ElementVector* FreeSurfaceBaseAnalysis::CreatePVector(Element* element){/*{{{*/
_error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
}

for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb) + dt*vz)*basis[i];
factor = Jdet*gauss->weight*(bed+dt*(mb) + dt*vz);
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];

if(stabilization==5){
/*SUPG*/
basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
if(dim==1){
vx_input->GetInputAverage(&vx);
tau=h/(2.*fabs(vx)+1e-10);
for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*mb+dt*vz)*tau*(vx*dbasis[0*numnodes+i]);
factor = Jdet*gauss->weight*(dt*mb+dt*vz)*tau;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*(vx*dbasis[0*numnodes+i]);
}
else{
vx_input->GetInputAverage(&vx);
vy_input->GetInputAverage(&vy);
tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed*0.+dt*mb+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
factor = Jdet*gauss->weight*(bed*0.+dt*mb+dt*vz)*tau;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
}
}
}
Expand Down
18 changes: 11 additions & 7 deletions src/c/analyses/FreeSurfaceTopAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -148,7 +148,7 @@ ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
int domaintype,dim,stabilization;
Element* topelement = NULL;
IssmDouble *xyz_list = NULL;
IssmDouble Jdet,D_scalar,dt,h;
IssmDouble Jdet,D_scalar,dt,h,factor;
IssmDouble vel,vx,vy,tau;

/*Get top element*/
Expand Down Expand Up @@ -282,17 +282,18 @@ ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
}
else if(stabilization==5){
D_scalar=gauss->weight*Jdet*dt;
factor = tau*D_scalar;
if(dim==2){
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j]+=tau*D_scalar*
Ke->values[i*numnodes+j]+=factor*
(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*
(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
}
}
}
else{
for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=factor*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
}
}
}
Expand All @@ -310,7 +311,7 @@ ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/

/*Intermediaries*/
int domaintype,dim,stabilization;
IssmDouble Jdet,dt;
IssmDouble Jdet,dt,factor;
IssmDouble ms,surface,vx,vy,vz,tau;
Element* topelement = NULL;
IssmDouble *xyz_list = NULL;
Expand Down Expand Up @@ -377,7 +378,8 @@ ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/
vz_input->GetInputValue(&vz,gauss);
surface_input->GetInputValue(&surface,gauss);

for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vz)*basis[i];
factor = Jdet*gauss->weight*(surface + dt*ms + dt*vz);
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}

if(stabilization==5){
Expand All @@ -386,13 +388,15 @@ ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/
if(dim==1){
vx_input->GetInputAverage(&vx);
tau=h/(2.*fabs(vx)+1e-10);
for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]);
factor = Jdet*gauss->weight*(dt*ms+dt*vz)*tau;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*(vx*dbasis[0*numnodes+i]);
}
else{
vx_input->GetInputAverage(&vx);
vy_input->GetInputAverage(&vy);
tau=h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
factor = Jdet*gauss->weight*(dt*ms+dt*vz)*tau;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
}
}

Expand Down
13 changes: 8 additions & 5 deletions src/c/analyses/GLheightadvectionAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ ElementMatrix* GLheightadvectionAnalysis::CreateKMatrix(Element* element){/*{{{*
/*Intermediaries */
const IssmPDouble yts = 365*24*3600.;
int domaintype,dim;
IssmDouble Jdet,D_scalar,onboundary;
IssmDouble Jdet,D_scalar,onboundary,factor;
IssmDouble vel,vx,vy;
IssmDouble* xyz_list = NULL;
Input* vx_input = NULL;
Expand Down Expand Up @@ -122,9 +122,10 @@ ElementMatrix* GLheightadvectionAnalysis::CreateKMatrix(Element* element){/*{{{*
vel = 30./yts*500000.;
}

factor = gauss->weight*Jdet;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
Ke->values[i*numnodes+j] += factor*(
(vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j])
+ vel/500000.*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]));
}
Expand All @@ -142,9 +143,10 @@ ElementMatrix* GLheightadvectionAnalysis::CreateKMatrix(Element* element){/*{{{*
/*Diffusion */
if(sqrt(vx*vx+vy*vy)<1000./31536000.){
IssmPDouble kappa = -10.;
factor = D_scalar*kappa;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
Ke->values[i*numnodes+j] += factor*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
}
}
}
Expand All @@ -158,8 +160,9 @@ ElementMatrix* GLheightadvectionAnalysis::CreateKMatrix(Element* element){/*{{{*

/*Artificial diffusivity*/
vel=sqrt(vx*vx + vy*vy)+1.e-14;
D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx); D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx); D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
factor = D_scalar*h/(2.*vel);
D[0][0]=factor*fabs(vx*vx); D[0][1]=factor*fabs(vx*vy);
D[1][0]=factor*fabs(vy*vx); D[1][1]=factor*fabs(vy*vy);
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += (
Expand Down
17 changes: 11 additions & 6 deletions src/c/analyses/HydrologyGlaDSAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -301,9 +301,10 @@ ElementMatrix* HydrologyGlaDSAnalysis::CreateKMatrix(Element* element){/*{{{*/
}

/*Diffusive term*/
IssmDouble factor = gauss->weight*Jdet;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
Ke->values[i*numnodes+j] += factor*(
coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j]
+ coeff*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
}
Expand All @@ -319,18 +320,20 @@ ElementMatrix* HydrologyGlaDSAnalysis::CreateKMatrix(Element* element){/*{{{*/
v1 = 0;
}
}
factor = gauss->weight*Jdet*(-v1);
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
Ke->values[i*numnodes+j] += factor*basis[i]*basis[j];
}
}

/*Transient term if dt>0*/
if(dt>0.){
/*Diffusive term*/
factor = gauss->weight*Jdet*e_v/(rho_water*g*dt);
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*basis[i]*basis[j];
Ke->values[i*numnodes+j] += factor*basis[i]*basis[j];
}
}
}
Expand Down Expand Up @@ -447,12 +450,14 @@ ElementVector* HydrologyGlaDSAnalysis::CreatePVector(Element* element){/*{{{*/
}
}

for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v2-m)*basis[i];
IssmDouble factor = - Jdet*gauss->weight*(w-v2-m);
for(int i=0;i<numnodes;i++) pe->values[i]+= factor*basis[i];

/*Transient term if dt>0*/
if(dt>0.){
factor = gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old;
phiold_input->GetInputValue(&phi_old,gauss);
for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
for(int i=0;i<numnodes;i++) pe->values[i] += factor*basis[i];
}
}

Expand Down Expand Up @@ -646,7 +651,7 @@ void HydrologyGlaDSAnalysis::UpdateSheetThickness(FemModel* femmodel){/*{{{*/
void HydrologyGlaDSAnalysis::UpdateSheetThickness(Element* element){/*{{{*/

/*Intermediaries */
IssmDouble Jdet,vx,vy,ub,h_old,N,h_r,H,b;
IssmDouble vx,vy,ub,h_old,N,h_r,H,b;
IssmDouble A,B,n,phi,phi_0;
IssmDouble alpha,beta;
IssmDouble oceanLS,iceLS;
Expand Down
16 changes: 8 additions & 8 deletions src/c/analyses/HydrologyShaktiAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -240,11 +240,12 @@ ElementMatrix* HydrologyShaktiAnalysis::CreateKMatrix(Element* element){/*{{{*/
IssmDouble pressure_water = rho_water*g*(head-bed);
if(pressure_water>pressure_ice) pressure_water = pressure_ice;

IssmDouble factor = conductivity*gauss->weight*Jdet;
IssmDouble factor2 = gauss->weight*Jdet*storage/dt + gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
+ gauss->weight*Jdet*storage/dt*basis[i]*basis[j]
+gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j];
Ke->values[i*numnodes+j] += factor*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
+ factor2*basis[i]*basis[j];
}
}
}
Expand Down Expand Up @@ -357,15 +358,14 @@ ElementVector* HydrologyShaktiAnalysis::CreatePVector(Element* element){/*{{{*/

meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));

for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
(
meltrate*(1/rho_water-1/rho_ice)
IssmDouble factor = Jdet*gauss->weight*
(meltrate*(1/rho_water-1/rho_ice)
+A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
+(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
-beta*sqrt(vx*vx+vy*vy)
+ieb
+storage*head_old/dt
)*basis[i];
+storage*head_old/dt);
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];

}

Expand Down
8 changes: 5 additions & 3 deletions src/c/analyses/HydrologyShreveAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -237,7 +237,7 @@ ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/

/*Intermediaries */
IssmDouble Jdet,dt;
IssmDouble mb,oldw;
IssmDouble mb,oldw,factor;
IssmDouble* xyz_list = NULL;

/*Fetch number of nodes and dof for this finite element*/
Expand Down Expand Up @@ -265,10 +265,12 @@ ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
oldw_input->GetInputValue(&oldw,gauss);

if(dt!=0.){
for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i];
factor = Jdet*gauss->weight*(oldw+dt*mb);
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}
else{
for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*mb*basis[i];
factor = Jdet*gauss->weight*mb;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}
}

Expand Down
3 changes: 2 additions & 1 deletion src/c/analyses/L2ProjectionBaseAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -196,7 +196,8 @@ ElementVector* L2ProjectionBaseAnalysis::CreatePVector(Element* element){/*{{{*/
input->GetInputValue(&value,gauss);
}

for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
IssmDouble factor = Jdet*gauss->weight*value;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}

/*Clean up and return*/
Expand Down
5 changes: 3 additions & 2 deletions src/c/analyses/L2ProjectionEPLAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,7 @@ ElementVector* L2ProjectionEPLAnalysis::CreatePVector(Element* element){/*{{{*/

/*Intermediaries */
int input_enum,index;
IssmDouble Jdet,slopes[2];
IssmDouble Jdet,slopes[2],factor;
Input *input = NULL;
IssmDouble *xyz_list = NULL;

Expand Down Expand Up @@ -204,7 +204,8 @@ ElementVector* L2ProjectionEPLAnalysis::CreatePVector(Element* element){/*{{{*/
basalelement->NodalFunctions(basis,gauss);

input->GetInputDerivativeValue(&slopes[0],xyz_list,gauss);
for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*slopes[index]*basis[i];
factor = Jdet*gauss->weight*slopes[index];
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}

/*Clean up and return*/
Expand Down
Loading

0 comments on commit 3a66e7c

Please sign in to comment.