Skip to content

Commit

Permalink
Add run-time options for removing certain terms. Fix definition of vo…
Browse files Browse the repository at this point in the history 
…lumetric heating terms.
  • Loading branch information
@jedbrown
jedbrown committed May 28, 2011
1 parent 5a30c92 commit c1a7aac
Show file tree
Hide file tree
Showing 3 changed files with 34 additions and 28 deletions.
25 changes: 13 additions & 12 deletions src/fs/tests/vht.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,7 +186,8 @@ static dErr VHTCaseSetFromOptions(VHTCase scase)
err = dOptionsRealUnits("-rheo_gravity","Gravity in z-direction (probably negative)","",u->Acceleration,rheo->gravity,&rheo->gravity,NULL);dCHK(err);
err = dOptionsRealUnits("-rheo_kinetic","Coefficient of kinetic energy used to determine internal energy","",NULL,rheo->kinetic,&rheo->kinetic,NULL);dCHK(err);
err = dOptionsRealUnits("-rheo_Kstab","Stabilization diffusivity, use to prevent negative diffusivity due to non-monotone splice","",NULL,rheo->Kstab,&rheo->Kstab,NULL);dCHK(err);
err = dOptionsRealUnits("-rheo_momentum_transport","Multiplier for the transport term in momentum balance","",NULL,rheo->momentum_transport,&rheo->momentum_transport,NULL);dCHK(err);
err = dOptionsRealUnits("-rheo_mask_momtrans","Multiplier for the transport term in momentum balance","",NULL,rheo->mask_momtrans,&rheo->mask_momtrans,NULL);dCHK(err);
err = dOptionsRealUnits("-rheo_mask_rho","Multiplier for density variation due to omega","",NULL,rheo->mask_rho,&rheo->mask_rho,NULL);dCHK(err);
if (scase->setfromoptions) {err = (*scase->setfromoptions)(scase);dCHK(err);}
} err = PetscOptionsEnd();dCHK(err);
dFunctionReturn(0);
Expand Down Expand Up @@ -1047,9 +1048,10 @@ static dErr VHTRheoViscosity(struct VHTRheology *rheo,dScalar p,VHTScalarD *T,VH
}
static void VHTStashGetRho(const struct VHTStash *st,const struct VHTRheology *rheo,VHTScalarD *rho)
{
rho->x = (1 - st->omega.x) * rheo->rhoi + st->omega.x * rheo->rhow;
rho->dp = (rheo->rhow - rheo->rhoi) * st->omega.dp;
rho->dE = (rheo->rhow - rheo->rhoi) * st->omega.dE;
dScalar mask = rheo->mask_rho;
rho->x = (1 - mask*st->omega.x) * rheo->rhoi + mask*st->omega.x * rheo->rhow;
rho->dp = (rheo->rhow - rheo->rhoi) * mask*st->omega.dp;
rho->dE = (rheo->rhow - rheo->rhoi) * mask*st->omega.dE;
}
static void VHTStashGetDui(const struct VHTStash *st,const struct VHTRheology *rheo,const dScalar drhou[9],dScalar Dui[6])
{
Expand Down Expand Up @@ -1106,8 +1108,7 @@ static dErr VHTPointwiseComputeStash(struct VHTRheology *rheo,const dScalar rhou

static dErr VHTPointwiseFunction(VHTCase scase,const dReal x[3],dReal weight,
const dScalar rhou[3],const dScalar drhou[9],const dScalar p[1],const dScalar dp[3],const dScalar E[1],const dScalar dE[3],
struct VHTStash *st,
dScalar rhou_[3],dScalar drhou_[6],dScalar p_[1],dScalar E_[1],dScalar dE_[3])
struct VHTStash *st, dScalar rhou_[3],dScalar drhou_[6],dScalar p_[1],dScalar E_[1],dScalar dE_[3])
{
struct VHTRheology *rheo = &scase->rheo;
dScalar frhou[3],fp[1],fE[1],ui[3],heatflux[3],Sigma,symstress[6],stress[9];
Expand All @@ -1122,16 +1123,16 @@ static dErr VHTPointwiseFunction(VHTCase scase,const dReal x[3],dReal weight,
for (dInt i=0; i<3; i++) heatflux[i] = -st->K[0]*dp[i] - st->K[1]*dE[i];
for (dInt i=0; i<6; i++) symstress[i] = st->eta.x * st->Dui[i] - (i<3)*p[0]; // eta Du - p I
dTensorSymUncompress3(symstress,stress);
Sigma = dColonSymScalar3(st->Dui,symstress); // Strain heating
for (dInt i=0; i<3; i++) rhou_[i] = -weight * frhou[i]; // Momentum forcing term
Sigma = st->eta.x*dColonSymScalar3(st->Dui,st->Dui); // Strain heating
for (dInt i=0; i<3; i++) rhou_[i] = -weight * frhou[i]; // Momentum forcing term
rhou_[2] -= weight * rho.x * rheo->gravity; // Gravitational source term
for (dInt i=0; i<3; i++)
for (dInt j=0; j<3; j++)
drhou_[i*3+j] = -weight * (rheo->momentum_transport*rhou[i]*st->u[j] - stress[i*3+j]);
drhou_[i*3+j] = -weight * (rheo->mask_momtrans*rhou[i]*st->u[j] - stress[i*3+j]);
p_[0] = -weight * (drhou[0]+drhou[4]+drhou[8] + fp[0]); // -q tr(drhou) - forcing, note tr(drhou) = div(rhou)
E_[0] = -weight * (Sigma + fE[0]); // Strain heating and thermal forcing
E_[0] -= weight * rhou[2] * rheo->gravity; // Gravitational source term for energy
for (dInt i=0; i<3; i++) dE_[i] = -weight * (ui[i]*E[0] + heatflux[i]); // Transport and diffusion
for (dInt i=0; i<3; i++) dE_[i] = -weight * (ui[i]*(E[0]+p[0]) + heatflux[i]); // Transport and diffusion
dFunctionReturn(0);
}
static dErr VHTPointwiseJacobian(const struct VHTStash *st,struct VHTRheology *rheo,dReal weight,
Expand All @@ -1151,7 +1152,7 @@ static dErr VHTPointwiseJacobian(const struct VHTStash *st,struct VHTRheology *r
for (dInt i=0; i<3; i++) ui1[i] = u1[i] - rheo->kappa_w*domega1[i] + st->wmom[i]/dSqr(rho.x)*rho1;

for (dInt i=0; i<3; i++) rhou_[i] = -weight * rho1 * (i==2) * rheo->gravity;
for (dInt i=0; i<3; i++) for (dInt j=0; j<3; j++) drhou_[i*3+j] = -weight * (rheo->momentum_transport*(rhou1[i]*st->u[j] + rho.x*st->u[i]*u1[j]) - Stress1[i*3+j]);
for (dInt i=0; i<3; i++) for (dInt j=0; j<3; j++) drhou_[i*3+j] = -weight * (rheo->mask_momtrans*(rhou1[i]*st->u[j] + rho.x*st->u[i]*u1[j]) - Stress1[i*3+j]);
p_[0] = -weight*(drhou1[0]+drhou1[4]+drhou1[8]); // -q tr(Du)
E_[0] = -weight*(Sigma1 - rhou1[2]*rheo->gravity);
for (dInt i=0; i<3; i++) dE_[i] = -weight * (ui[i] * (E1[0]+p1[0]) + ui1[i] * (st->E + st->p)
Expand All @@ -1168,7 +1169,7 @@ static void VHTPointwiseJacobian_uu(const struct VHTStash *st,const struct VHTRh
VHTStashGetRho(st,rheo,&rho);
VHTStashGetStress1(st,rheo,drhou1,p1,E1,Stress1,NULL);
for (dInt i=0; i<3; i++) u1[i] = rhou1[i] / rho.x; // perturb u = rhou/rho
for (dInt i=0; i<3; i++) for (dInt j=0; j<3; j++) drhou_[i*3+j] = -weight * (rheo->momentum_transport*(rhou1[i]*st->u[j] + rho.x*st->u[i]*u1[j]) - Stress1[i*3+j]);
for (dInt i=0; i<3; i++) for (dInt j=0; j<3; j++) drhou_[i*3+j] = -weight * (rheo->mask_momtrans*(rhou1[i]*st->u[j] + rho.x*st->u[i]*u1[j]) - Stress1[i*3+j]);
}

static void VHTPointwiseJacobian_pu(dReal weight,const dScalar drhou1[9],dScalar p_[1])
Expand Down
34 changes: 19 additions & 15 deletions src/fs/tests/vhtexact.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
from __future__ import division

import sympy
from sympy import Matrix, Symbol, ccode
from sympy import Matrix, Symbol, ccode, S
from dohpexact import *

# Some coupling terms make the symbolic expressions explode in size, which takes too long and crashes the compiler. Such
Expand All @@ -11,6 +11,7 @@
# numerical rootfinding problem. Since we cannot symbolically differentiate an implicit solve, we are shorting it out in
# a physically plausible way, but one that is inconsistent for large amplitude moisture.
MASK = 0
MASK2 = 0

def SecondInvariant(Du):
return 0.5*Du.dot(Du)
Expand All @@ -37,21 +38,20 @@ def separate(x,w):
neg = integrate(0.5+integrate(normal(x,a),x),x)
pos = integrate(0.5-integrate(normal(x,a),x),x)
"""
from sympy import erf,exp,pi
from sympy import erf,exp,pi,sqrt
# This explicit CSE does not seem to make a performance difference with sympy, but it probably would if we use AD
theerf = erf(x*2**(1/2)*(w**(-2))**(1/2)/2)
xtheerf = x*theerf
theerf = erf(x/sqrt(2*w**2))
theexp = exp(-x**2/(2*w**2))
neg = 0.5*x + (xtheerf + 2**(1/2)*w**2*(w**(-2))**(1/2)*theexp/pi**(1/2))/(2*(w**(-2))**(1/2)*(w**2)**(1/2))
pos = 0.5*x - (xtheerf + 2**(1/2)*w**2*(w**(-2))**(1/2)*theexp/pi**(1/2))/(2*(w**(-2))**(1/2)*(w**2)**(1/2))
neg1 = 0.5 + theerf/(2*(w**(-2))**(1/2)*(w**2)**(1/2)) # first derivwtives
pos1 = 0.5 - theerf/(2*(w**(-2))**(1/2)*(w**2)**(1/2))
neg = 0.5*x*(1 - theerf) - w*theexp/sqrt(2*pi)
pos = 0.5*x*(1 + theerf) + w*theexp/sqrt(2*pi)
neg1 = 0.5*(1 - theerf) # first derivwtives
pos1 = 0.5*(1 + theerf)
return neg,pos,neg1,pos1

class VHTExact(Exact):
def __init__(self, name=None, model=None, param='a b c'):
if model is None:
model = 'B0 Bomega R Q V T0 eps gamma0 pe beta_CC rhoi rhow T3 c_i Latent splice_delta k_T kappa_w gravity Kstab momentum_transport'
model = 'B0 Bomega R Q V T0 eps gamma0 pe beta_CC rhoi rhow T3 c_i Latent splice_delta k_T kappa_w gravity Kstab mask_momtrans'
Exact.__init__(self, name=name, model=model, param=param, fieldspec=[('rhou',3), ('p',1), ('E',1)])
def unpack(self, U, dU):
rhou, p, E = U[:3,:], U[3], U[4]
Expand All @@ -63,10 +63,10 @@ def solve_eqstate(self, rhou, p, E, drhou, dp, dE):
T_m = T3 - beta_CC * p # melting temperature at current pressure
e_m = c_i * (T_m - T0) # melting energy at current pressure
# At this point, we should solve an implicit system for (e,rho,omega). Unfortunately, doing so would require
# solving an inhomogeneous system involving splice(), which I don't know how to do symbolically. Putting a
# solving an inhomogeneous system involving separate(), is problematic to do symbolically. Putting a
# numeric rootfinder into this symbolic description would make manufactured solutions much more complex.
# Therefore, we cheat by simply using ice density to remove kinetic energy and convert energy/volume to
# energy/mass.
# energy/mass. Note that with algorithmic differentiation, the implicit solve would not be a problem.
rhotmp = rhoi # Cheat
e = (E - MASK*1/(2*rhotmp) * rhou.dot(rhou)) / rhotmp
de = (dE - MASK*1/(rhotmp) * (rhou.T * drhou).T) / rhotmp
Expand All @@ -85,8 +85,12 @@ def meltfraction(e):
return omega.subs(x,e), omega1.subs(x,e)
omega, omega1 = meltfraction(e4omega1)
domega = omega1 * e4omega1 * (G1p*dp + G1E*dE)
rho = (1-omega)*rhoi + omega*rhow
drho = (rhow-rhoi) * domega
if MASK2 == 0: # Shameful, but SymPy chokes
rho = rhoi
drho = 0*domega
else:
rho = (1-omega)*rhoi + omega*rhow
drho = (rhow-rhoi) * domega
return (e,T,omega,rho), (de,dT,domega,drho)
def eta(self, p, T, omega, gamma): # Power law with Arrhenius relation
from sympy import exp
Expand All @@ -98,7 +102,7 @@ def eta(self, p, T, omega, gamma): # Power law with Arrhenius relation
def weak_homogeneous(self, x, U, dU, V, dV):
(rhou,p,E), (drhou,dp,dE) = self.unpack(U,dU)
(e,T,omega,rho), (de,dT,domega,drho) = self.solve_eqstate(rhou,p,E,drhou,dp,dE)
k_T, kappa_w, Kstab, L, grav, momtrans = self.model_get('k_T kappa_w Kstab Latent gravity momentum_transport')
k_T, kappa_w, Kstab, L, grav, momtrans = self.model_get('k_T kappa_w Kstab Latent gravity mask_momtrans')
gravvec = grav*Matrix([0,0,1])
u = rhou / rho # total velocity
du = (1/rho) * drhou - MASK*u * drho.T # total velocity gradient
Expand All @@ -112,7 +116,7 @@ def weak_homogeneous(self, x, U, dU, V, dV):
(rhou_,p_,E_), (drhou_,dp_,dE_) = self.unpack(V,dV)
conserve_momentum = -drhou_.dot(momtrans*rhou*u.T - eta*Dui + p*I) - rhou_.dot(rho*gravvec)
conserve_mass = -p_ * drhou.trace()
conserve_energy = -dE_.dot((E+p)*ui + heatflux - Kstab*dE) - Dui.dot(eta*Dui) - E_ * rhou.dot(gravvec)
conserve_energy = -dE_.dot((E+p)*ui + heatflux - Kstab*dE) -E_*(eta*Dui.dot(Dui) + rhou.dot(gravvec))
return conserve_momentum + conserve_mass + conserve_energy
def create_prototype(self, definition=False):
return 'dErr VHTCaseCreate_%(name)s(VHTCase case)%(term)s' % dict(name=self.name, term=('' if definition else ';'))
Expand Down
3 changes: 2 additions & 1 deletion src/fs/tests/vhtimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,8 @@ struct VHTRheology {
dReal gravity; /* Strength of gravity in z direction (probably negative) */
dReal kinetic; /* Parameter to turn on the use of kinetic energy when computing velocity */
dReal Kstab; /* Stabilization for energy diffusion */
dReal momentum_transport; /* Multiplier for the transport term in momentum balance */
dReal mask_momtrans; /* Multiplier for the transport term in momentum balance */
dReal mask_rho; /* Multiplier for the true rho */
};
struct VHTUnitTable {
dUnit Length;
Expand Down

0 comments on commit c1a7aac

Please sign in to comment.