Toddler in Fenics

variational formulation
1

Hi Folks,

I am trying to adapt two-phase flow in a poroelastic medium. I am just a beginner and I am actually struggling to understand the bits and bobs of dolfin and ufl. Anyway, I really appreciate it if I get some advice for this code (it gives an error on the assemble point for A matrix), the error is UFLException: Component and shape length don’t match.:

def WeakForm(U,V,X_n,t):

p, q, us, sw = split(U)
Pt, Qt, Ust, Swt = split(V)
p_n, q_n, us_n, sw_n = split(X_n)
p_m, q_m, us_m, sw_m = split(X_m)


### Weak Forms ###
# Conservation of Mass (Pressure Equation)
TwoPhase_P_1 = (dt/(rho_w(p_m,p_0,p_h,rho0_w,beta_w)*phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)))*\
               div(rho_w(p_m,p_0,p_h,rho0_w,beta_w)*k_mod(k0, exponent, alpha,us,p_m,beta_s,phi_0,phi_min,t)*\
               mobility_w(vis_w, sw_m, lambda_kr, Sw_ir)*div(p - Pc(sw_m, m , a, Sw_ir)))*Pt*dx 
TwoPhase_P_2 = -((1-sw_n)/rho_g(p_m,p_0,p_h,rho0_g,beta_g))*drho_gP(p_m,p_0,p_h,rho0_g,beta_g)*p*Pt*dx
TwoPhase_P_3 = (dt/(rho_g(p_m,p_0,p_h,rho0_g,beta_g)*phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)))*\
               div(rho_g(p_m,p_0,p_h,rho0_g,beta_g)*k_mod(k0, exponent, alpha,us,p_m,beta_s,phi_0,phi_min,t)*\
               mobility_g(vis_g, sw_m, lambda_kr, Sw_ir)*div(p))*Pt*dx

TwoPhase_Pr1 = -((1-sw_n)/rho_g(p_m,p_0,p_h,rho0_g,beta_g))*drho_gP(p_m,p_0,p_h,rho0_g,beta_g)*p_n*Pt*dx
TwoPhase_Pr2 = - (dt/(rho_w(p_m,p_0,p_h,rho0_w,beta_w) * phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)))*Qww *Pt*dx
TwoPhase_Pr3 = - (dt/(rho_g(p_m,p_0,p_h,rho0_g,beta_g) * phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)))*Qgg *Pt*dx

TwoPhase_P_lhs = TwoPhase_P_1 + TwoPhase_P_2 + TwoPhase_P_3
TwoPhase_P_rhs = TwoPhase_Pr1 + TwoPhase_Pr2 + TwoPhase_Pr3
    
# Conservation of Mass (Saturation Equation)
TwoPhase_S_1 = sw*Swt*dx
TwoPhase_S_2 = -(dt/(rho_w(p_m,p_0,p_h,rho0_w,beta_w) * phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)))*\
               div(rho_w(p_m,p_0,p_h,rho0_w,beta_w) * k_mod(k0, exponent, alpha,us,p_m,beta_s,phi_0,phi_min,t)*\
               mobility_w(vis_w, sw_n, lambda_kr, Sw_ir)* div(p - Pc(sw_n, m , a, Sw_ir))) * Swt*dx

TwoPhase_Sr1 = sw_n*Swt*dx
TwoPhase_Sr2 = (dt/(rho_w(p_m,p_0,p_h,rho0_w,beta_w) * phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)))*Qww *Swt*dx

TwoPhase_S_lhs = TwoPhase_S_1 + TwoPhase_S_2
TwoPhase_S_rhs = TwoPhase_Sr1 + TwoPhase_Sr2


CoMass_l_1 = rho_f(p_m,p_0,p_h,rho_0,beta_f)*Constant(alpha)*epsilon_v(us)*Pt*dx
CoMass_l_2 = rho_f(p_m,p_0,p_h,rho_0,beta_f)*((Constant(alpha) - phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t))*Constant(beta_s)\
    + phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)*Constant(beta_f))*p*Pt*dx
CoMass_l_3 = dt*rho_f(p_m,p_0,p_h,rho_0,beta_f)*div(q)*Pt*dx
CoMass_l_4 = dt*Constant(weight)*inner(q,grad(rho_f(p_m,p_0,p_h,rho_0,beta_f)*ones_func))*Pt*dx
CoMass_l = CoMass_l_1 + CoMass_l_2 + CoMass_l_3 + CoMass_l_4

CoMass_r_1 = dt*(Qww+Qgg)*Pt*dx
CoMass_r_2 = rho_f(p_m,p_0,p_h,rho_0,beta_f)*Constant(alpha)*epsilon_v(us_n)*Pt*dx
CoMass_r_3 = rho_f(p_m,p_0,p_h,rho_0,beta_f)*((Constant(alpha)\
    -phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t))*Constant(beta_s) \
    + phi(alpha,us_m,p_m,beta_s,phi_0,phi_min,t)*Constant(beta_f))*p_n\
    *Pt*dx
CoMass_r = CoMass_r_1 + CoMass_r_2 + CoMass_r_3

# Darcy's Law
DL_l = mu/k_mod(k0, exponent, alpha,us,p,beta_s,phi_0,phi_min,t)*inner(q,Qt)*dx - p*div(Qt)*dx
DL_r = -Constant(100)*inner(Qt,norm)*ds(6)

# Conservation of Momentum
CoMom_l = inner(sigma(us,p,alpha,K,G),grad(Ust))*dx
# CoMom_r = inner(Constant((0.0,***loading***)),Ust)*ds(-)

A =  TwoPhase_P_lhs + TwoPhase_S_lhs + CoMass_l + CoMom_l + DL_l    
B =  TwoPhase_P_rhs + TwoPhase_S_rhs + CoMass_r + DL_r #+ CoMom_r 

return A,B

def LinearSolver(U,V,X_n,t,bcs):
a,L = WeakForm(U,V,X_n,t)
b = assemble(L)
A = assemble(a)
A, b = assemble_system(a,L,bcs)
solve(A, X.vector(), b, ‘mumps’)
return X

2

Please make sure that you follow Read before posting: How do I get my question answered?
and look at the minimal working examples there.

It is very hard to help you with your problem ( since the code is not executable as the function spaces, meshes etc are not defined).