[cig-commits] commit: updated MADDs-4 & 5 P1 forms
Mercurial
hg at geodynamics.org
Fri Dec 11 11:36:33 PST 2009
changeset: 112:7feb2fbe6d2f
tag: tip
user: Marc Spiegelman <mspieg at ldeo.columbia.edu>
date: Fri Dec 11 14:31:48 2009 -0500
files: MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl MADDs-5/MADDs-5a/cpp/forms/MagmaEquations_P1.ufl MADDs-5/MADDs-5a/cpp/forms/MeltFlux_P1.ufl MADDs-5/MADDs-5a/cpp/forms/SemiLagrangianRHS_P1.ufl
description:
updated MADDs-4 & 5 P1 forms
diff -r 93fed2bfd487 -r 7feb2fbe6d2f MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl
--- a/MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl Fri Dec 11 13:33:52 2009 -0500
+++ b/MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl Fri Dec 11 14:31:48 2009 -0500
@@ -31,7 +31,7 @@ hsquared = Constant(triangle) # (h/delta
hsquared = Constant(triangle) # (h/delta)^2
# facet normal for forcing term
-n = P2.cell().n
+n = P1.cell().n
#Permeability
K = f*f*f
diff -r 93fed2bfd487 -r 7feb2fbe6d2f MADDs-5/MADDs-5a/cpp/forms/MagmaEquations_P1.ufl
--- a/MADDs-5/MADDs-5a/cpp/forms/MagmaEquations_P1.ufl Fri Dec 11 13:33:52 2009 -0500
+++ b/MADDs-5/MADDs-5a/cpp/forms/MagmaEquations_P1.ufl Fri Dec 11 14:31:48 2009 -0500
@@ -13,15 +13,19 @@
# c = (h/delta)^2
# Gamma = melting rate function (assumed W dF/dz)
#
-# Compile this form with FFC: ffc -l dolfin MagmaEquation.ufl
+# Compile this form with FFC: ffc -l dolfin -O MagmaEquation.ufl
P1 = FiniteElement("Lagrange", triangle, 1)
P2 = FiniteElement("Lagrange", triangle, 2)
-M = P1 + P1
+ME = P1 + P1
-(v, u) = TestFunctions(M)
-(dp, df) = TrialFunctions(M) # deltas for pressure and porosity
-(p, f) = Functions(M) # solution from last iteration
+(v, q) = TestFunctions(ME)
+du = TrialFunction(ME) # deltas for pressure and porosity
+u = Function(ME) # solution from last iteration
+
+# Split mixed functions
+dp, df = split(du)
+p, f = split(u)
#parameters and functions
@@ -32,9 +36,8 @@ hdt = 0.5*dt # half time step
# set domain size coefficient h=h/delta height in compaction lengths
hsquared = Constant(triangle) # (h/delta)^2
-
# facet normal for forcing term
-n = VectorConstant(triangle)
+n = P1.cell().n
# inverse shear viscosity
iEta = Function(P1)
@@ -47,33 +50,23 @@ PStar = Function(P1)
# regularization constants
-#eps = 5.e-3
eps = Constant(triangle)
zero = 0.
# permeability function
-def K(f, eps):
- return (f+eps)*(f+eps)
-
-# derivatives of permeability
-def dKdf(f,eps):
- return 2*(f+eps)
+K = (f+eps)*(f+eps) # regularized
+K0 = f*f # not regularized
# bulk viscosity function
-def Xi(f, eps, iEta):
- return iEta*(f + eps)
+Xi = hsquared*iEta*(f+eps)
-# derivative of bulk viscosity with respect to porosity
-def dXidf(iEta):
- return iEta
-
-a1 = (K(f,eps)*inner(grad(v),grad(dp)) + hsquared*v*Xi(f,eps,iEta)*dp)*dx + ((dKdf(f,eps)*inner(grad(v),grad(p)) + hsquared*dXidf(iEta)*v*p - dKdf(f,eps)*v.dx(1))*df)*dx + dKdf(f,eps)*df*v*n[1]*ds + dKdf(f,eps)*df*inner(grad(PStar),grad(v))*dx - dKdf(f,eps)*df*v*inner(grad(PStar),n)*ds
-a2 = u*df*dx - u*hdt*hsquared*(Xi(f,eps,iEta)*dp + p*dXidf(iEta)*df)*dx
-L1 = (K(f,eps)*inner(grad(v),grad(p)) + hsquared*v*Xi(f,eps,iEta)*p - K(f,zero)*v.dx(1))*dx + K(f,zero)*v*n[1]*ds + K(f,zero)*inner(grad(PStar),grad(v))*dx - K(f,zero)*v*inner(grad(PStar),n)*ds
-L2 = u*(f - hdt*(hsquared*Xi(f,eps,iEta)*p + Gamma))*dx
+# calculate linear form for residual F(u)
+L1 = (K*inner(grad(v),grad(p)) + v*Xi*p - K0*v.dx(1))*dx + K0*inner(grad(PStar),grad(v))*dx + K0*v*(n[1] - inner(grad(PStar),n))*ds
+L2 = q*(f - hdt*(Xi*p + Gamma))*dx
-#correct sign of L in Dolfin >=0.9.0
-a = a1 + a2
L = L1 + L2
+# bilinear form for weak form of the Jacobian J(u) (using functional differentiation of L)
+a = derivative(L, u, du)
+
diff -r 93fed2bfd487 -r 7feb2fbe6d2f MADDs-5/MADDs-5a/cpp/forms/MeltFlux_P1.ufl
--- a/MADDs-5/MADDs-5a/cpp/forms/MeltFlux_P1.ufl Fri Dec 11 13:33:52 2009 -0500
+++ b/MADDs-5/MADDs-5a/cpp/forms/MeltFlux_P1.ufl Fri Dec 11 14:31:48 2009 -0500
@@ -15,7 +15,7 @@
# The corresponding weak (variational problem) is a projection from FiniteElement spaces for V and p onto
# that of v assuming v can be described by BDM elements of degree q
#
-# Compile this form with FFC: ffc -l dolfin MeltFlux_P1.ufl
+# Compile this form with FFC: ffc -l dolfin MeltFlux_P2.ufl
#
#Modified from the PyDolfin MADDs-2.py demo
#
@@ -28,7 +28,9 @@
# set up pressure/velocity function spaces
VP2 = VectorElement("Lagrange",triangle, 2)
VP1 = VectorElement("Lagrange",triangle,1)
+P2 = FiniteElement("Lagrange",triangle, 2)
P1 = FiniteElement("Lagrange",triangle, 1)
+
@@ -42,8 +44,8 @@ V = Function(VP2)
# regularized permeability function
eps = Constant(triangle)
zero = 0.
-def K(f, eps):
- return (f+eps)*(f+eps)
+K = (f+eps)*(f+eps)
+K0 = f*f
# define new BDM function spaces
@@ -54,5 +56,5 @@ u = TrialFunction(VP1)
u = TrialFunction(VP1)
a = inner(w,u)*dx
-L = f*inner(w,V)*dx + K(f,eps)*(-inner(w,grad(p)))*dx + K(f,zero)*(-inner(w,grad(pstar))+w[1])*dx
+L = f*inner(w,V)*dx + K*(-inner(w,grad(p)))*dx + K0*(-inner(w,grad(pstar))+w[1])*dx
diff -r 93fed2bfd487 -r 7feb2fbe6d2f MADDs-5/MADDs-5a/cpp/forms/SemiLagrangianRHS_P1.ufl
--- a/MADDs-5/MADDs-5a/cpp/forms/SemiLagrangianRHS_P1.ufl Fri Dec 11 13:33:52 2009 -0500
+++ b/MADDs-5/MADDs-5a/cpp/forms/SemiLagrangianRHS_P1.ufl Fri Dec 11 14:31:48 2009 -0500
@@ -9,7 +9,7 @@ QE = FiniteElement("Quadrature", triangl
QE = FiniteElement("Quadrature", triangle,2) # quadrature element for Semi-Lagrangian source term
u = TestFunction(P1)
-gStar = Function(P1) # source term calculated from last time step by semi-Lagrangian update
+gStar = Function(QE) # source term calculated from last time step by semi-Lagrangian update
L = -u*gStar*dx
More information about the CIG-COMMITS
mailing list