I am trying to solve a system of pde. I am having problems updating my solutions. Bellow is my code and the corresponding error. Thank you in advance for your answers:
I think the problem is in the expressions:
p_k.assign(p)
h_k.assign(h)
in the while loop at the end of the program
mesh2d = UnitSquareMesh(40,40)
mesh2d = UnitIntervalMesh(120)
mesh with geothermal well in the center
mesh2d = Mesh(domain, 45)
plot(mesh2d, "2D mesh")
setting up mixt function space
P = FunctionSpace(mesh2d, "Lagrange",2)
H = FunctionSpace(mesh2d, "Lagrange",2)
M = MixedFunctionSpace([H,P])
Test = TestFunction(M)
(v,q) = split(Test)
constant and parameters for the problem
phi = 0.1 #porosity
k = 1e-15 #permeability, m^2
k_th = 2 #thermal conductivity, W/m/K
rho_s = 2900. #rock density, kg/m^3
c_s = 840.
h0 = Expression("1790.47") #boundary enthalpy
p0 = Expression("19290.69") # boundary pressure
boundary condition for enthalpy h and pressure p
pressure p and enthalpy are set constant at the boundary
def boundary_h(x,on_boundary):
return on_boundary
bc1 = DirichletBC(M.sub(0),h0,boundary_h)
def boundary_p(x,on_boundary):
return on_boundary
bc2 = DirichletBC(M.sub(1),p0,boundary_p)
bcs =[bc1,bc2]
h_k and p_k are the solutionat the preview time level
h_init and p_init are the initial enthalpy and pressure in the domain
h_init = Expression("1790754.47") # corresponding enthalpy to hydrostatic pressure
p_init = Expression("19290097.69 ") # hydrostatic pressure
h_k = project(h_init,H)
p_k = project(p_init,P)
h and p are the unknow at the new time level
solution = Function(M)
h,p = split(solution)
source term - correspond to production
source term + correspond to injection
a = 0.3
f = Expression("100")
g = Expression("0.0")
dt = 0.3
unsteady state variational formulation
#
variational formulation for conservation of mass
a1 =(phidrdpp)vdx + dt( inner( ((rhok)/mu)grad(p),grad(v) ))dx
L1 = (phidrdpp_k+phidrdhh_k)vdx +dtfvdx -phidrdhh_kv*dx
variational formulation of energy equation
a2 = (phih_kdrdp+(1-phi)rho_sc_sdTdp)pqdx # h was linearised to h_k
+(phirho+phihdrdh)qdx +((1-phi)rho_sc_sdTdh)hq*dx
+dtinner( ((rhokh)/mu)grad(p),grad(q))*dx
+dtinner( k_thdTdhgrad(h),grad(q) )dx
dtinner( k_thdTdpgrad(p),grad(q) )dx
L2 = (phihdrdp+(1-phi)rho_sc_sdTdp)p_kqdx
+( phirho+phihdrdh+(1-phi)rho_sc_sdTdh)h_kqdx
+dtgqdx
function F for Newton method
a = a1+a2
L = L1+L2
F = a-L
h_k = Function(H)
p_k = Function(P)
T=1
t=0
while(t<=T):
h0.t = t
p0.t = t
solve(F==0,solution,bcs)
h,p = split(solution)
p_k.assign(p)
h_k.assign(h)
t +=dt
plot(p)
plot(h)
interactive()
Traceback (most recent call last):
File "unsteady.py", line 114, in
p_k.assign(p)
File "/usr/lib/python2.7/dist-packages/dolfin/functions/function.py", line 410, in assign
"Expects a Function or linear combinations of "\
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/common.py", line 2294, in dolfin_error
return _common.dolfin_error(*args)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics@fenicsproject.org
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to function assignment.
*** Reason: Expects a Function or linear combinations of Functions in the same FunctionSpaces.
*** Where: This error was encountered inside function.py.
*** Process: 0
*** DOLFIN version: 1.3.0
*** Git changeset: unknown
