Heat MHD Cavity flow simulation error

Errors
1

Thank you for watching. This is the program used to simulate the cavity flow, but it still makes mistakes after a long time of modification. Please publish it again. I hope you can find the mistakes or make comments.

from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
import sympy as sym
from sympy import cos, sin
import numpy as np
import math
import time
from mshr import *

def boundary(x):
    tol = 1E-12
    return abs(x[0]) < tol or abs(x[1]) < tol or abs(x[0] - 1) < tol or abs(x[1] - 1) < tol


# Define two-dimensional versions of cross and curl operators
def scross(x, y):
    return x[0] * y[1] - x[1] * y[0]  # 向量叉乘


def vcross(x, y):
    return as_vector([x[1] * y, -x[0] * y])


def scurl(x):
    return x[1].dx(0) - x[0].dx(1)


def vcurl(x):
    return as_vector([x.dx(1), -x.dx(0)])


def acurl(x):
    return as_vector([x[2].dx(1),-x[2].dx(0),x[1].dx(0) - x[0].dx(1)])


def E_L2_H1(n):

    channel = Rectangle(Point(0, 0), Point(1.0, 1.0))
    mesh = generate_mesh(channel, n)

    P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
    B = FiniteElement("Bubble", mesh.ufl_cell(), 3)
    V = VectorElement(P1 + B)  # velocity field space
    # V = VectorElement("P", mesh.ufl_cell(), 2)
    Q = FiniteElement("P", mesh.ufl_cell(), 1)  # pressure field space
    M = FiniteElement("N1curl", mesh.ufl_cell(), 1)  # magnetic filed space
    G = FiniteElement("P", mesh.ufl_cell(), 1)
    H = FiniteElement("P", mesh.ufl_cell(), 1)  # temperture

    We = MixedElement([V, Q, M, G, H])  # u p B r T
    W = FunctionSpace(mesh, We)

    T = 11  # final time
    Re = 1000
    Sc = 1
    Rm = 1
    dt = 1 / 100
    t = dt

    f = Constant((0.0, 0.0))
    g = Constant((0.0, 0.0))
    beta = Constant((1.0, 0.0))
    k = Constant(1.0)
    ps=Constant((0.0, 0.0))

    u1_d1 = '''0'''
    u2_d1 = '''0'''
    B_d1 = '''0'''
    B_d2 = '''0'''
    p_d1 = '''0'''
    r_x = '''0'''
    T_D = '''0'''


    W_ex = Expression((u1_d1, u2_d1, B_d1, B_d2, p_d1, r_x, T_D), degree=4)
    U = interpolate(W_ex, W)
    Uk = Function(W)  # 储存上一步值的函数空间 定义
    start = time.time()
    while t <= T:
        inflow = 'near(x[0],0)'
        outflow = 'near(x[0],1)'
        walls = 'near(x[1],0) || near (x[1],1)'
        inflow_profile = ('6.0', '0')
        bcu_inflow = DirichletBC(W.sub(0), Expression(inflow_profile, degree=2), inflow)
        bcu_outflow = DirichletBC(W.sub(0), Expression(inflow_profile, degree=2), outflow)
        bcu_walls = DirichletBC(W.sub(0), Constant((0.0, 0.0)), walls)
        bcB_inflow = DirichletBC(W.sub(2), Constant((0.0, 0.0)), inflow)
        bcB_walls = DirichletBC(W.sub(2), Constant((1.0, 0.0)), walls)
        bcB_outflow = DirichletBC(W.sub(2), Constant((0.0, 0.0)), outflow)
        bcr = DirichletBC(W.sub(3), Constant(0.0), boundary)
        bcT_inflow = DirichletBC(W.sub(4), Constant(1.0), inflow)
        #bcp_outflow = DirichletBC(W.sub(1), Constant(0.0), outflow)
        bcs = [bcu_inflow, bcu_outflow, bcu_walls, bcB_inflow, bcB_walls, bcB_outflow, bcr , bcT_inflow]

        u, p, B, r, TTT = TrialFunctions(W)
        v, q, s, w, ttt = TestFunctions(W)
        (uk, pk, Bk, rk, TTTk) = Uk.split(deepcopy=True)

        # temperture-inner(beta*TTT,v)
        F = (dot((u - uk) / dt, v) + inner(grad(u), grad(v)) * pow(Re, -1) + 0.5 * inner(dot(grad(u), uk),v)- 0.5 * inner(dot(grad(v), uk), u) - p * div(v) - dot(scurl(B), scross(Bk, v)) * Sc + div(u) * q + 1.E-10 * inner(p, q)- inner(beta * TTT, v)) * dx \
            + (dot((B - Bk) / dt, s) + dot(scurl(B), scurl(s)) * Sc * pow(Rm, -1) - dot(scross(u, Bk),scurl(s)) * Sc - inner(grad(r),s)+ inner(grad(w), B)) * dx \
            + (dot((TTT - TTTk) / dt, ttt) + inner(grad(TTT), grad(ttt)) * k + 0.5 * inner(dot(grad(TTT), uk),ttt) - 0.5 * inner(dot(grad(ttt), uk), TTT)) * dx - (inner(f, v) + inner(g, s) + inner(ps, ttt)) * dx

        A0 = lhs(F)
        L0 = rhs(F)
        U = Function(W)
        solve(A0 == L0, U, bcs)
        (u, p, B, r, TTT) = U.split(deepcopy=True)
        plot(u, title='t=%.d Velocity' % t)
        plt.show()

    end = time.time()
    print("time：%。2f秒" % (end - start))

n = 16
E_L2_H1(n)

···

The error 
Shapes do not match: <Coefficient id=139856827226432> and <Indexed id=139856696760336>.
Traceback (most recent call last):
  File "/mnt/c/Users/ALSACE/PycharmProjects/pythonProject3/main.py", line 118, in <module>
    E_L2_H1(n)
  File "/mnt/c/Users/ALSACE/PycharmProjects/pythonProject3/main.py", line 104, in E_L2_H1
    + (dot((TTT - TTTk) / dt, ttt) + inner(grad(TTT), grad(ttt)) * k + 0.5 * inner(dot(grad(TTT), uk),ttt) - 0.5 * inner(dot(grad(ttt), uk), TTT)) * dx - (inner(f, v) + inner(g, s) + inner(ps, ttt)) * dx
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
    return Inner(a, b)
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
    error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <Coefficient id=139856827226432> and <Indexed id=139856696760336>.

Thanks again!
2

This is your porblem, ps is a 2D vector

TTT is from

a scalar function space.
Thus the inner product of these functions does not make sense.

1 Like
3

Thanks Mr.Dokken
Problem solved, ready to run.

4

Now I met a new problem, change the coefficient but the image does not change, I suspect that there is something wrong with the cycle, but I can’t find it, anyone have any good opinions?

from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
import sympy as sym
from sympy import cos, sin
import numpy as np
import math
import time
from mshr import *


def boundary(x):
    tol = 1E-12
    return abs(x[0]) < tol or abs(x[1]) < tol or abs(x[0] - 1) < tol or abs(x[1] - 1) < tol


# Define two-dimensional versions of cross and curl operators
def scross(x, y):
    return x[0] * y[1] - x[1] * y[0]  # 向量叉乘

def vcross(x, y):
    return as_vector([x[1] * y, -x[0] * y])

def scurl(x):
    return x[1].dx(0) - x[0].dx(1)

def vcurl(x):
    return as_vector([x.dx(1), -x.dx(0)])

def acurl(x):
    return as_vector([x[2].dx(1),-x[2].dx(0),x[1].dx(0) - x[0].dx(1)])


def E_L2_H1(n):

    #channel = Rectangle(Point(0, 0), Point(1.0, 1.0))
    #mesh = generate_mesh(channel, n)
    mesh = UnitSquareMesh(n,n)
    #plot(mesh)

    P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
    B = FiniteElement("Bubble", mesh.ufl_cell(), 3)
    V = VectorElement(P1 + B)  # velocity field space
    # V = VectorElement("P", mesh.ufl_cell(), 2)
    Q = FiniteElement("P", mesh.ufl_cell(), 1)  # pressure field space
    M = FiniteElement("N1curl", mesh.ufl_cell(), 1)  # magnetic filed space
    G = FiniteElement("P", mesh.ufl_cell(), 1)
    H = FiniteElement("P", mesh.ufl_cell(), 1)  # temperture

    We = MixedElement([V, Q, M, G, H])  # u p B r T
    W = FunctionSpace(mesh, We)

    T = 2  # final time
    #Re = 1000
    Pru = 1
    Ra = 100
    Sc  = 1
    #Rm = 1
    Prb = 1
    Prt = 0.71
    dt  = 1 / 100
    t = dt

    f = Constant((0.0, 0.0))
    g = Constant((0.0, 0.0))

    #beta = Expression((0.0, 1.0), degree=10, t=t)

    beta = Constant((0.0, 1.0))
    k = Constant(1.0)
    ps = Constant(0.0)

    u1_d1 = '''0'''
    u2_d1 = '''0'''
    B_d1 = '''0'''
    B_d2 = '''0'''
    p_d1 = '''0'''
    r_x = '''0'''
    T_D = '''0'''

    W_ex = Expression((u1_d1, u2_d1, B_d1, B_d2, p_d1, r_x, T_D), degree=4)
    U = interpolate(W_ex, W)
    Uk = Function(W)  # 储存上一步值的函数空间
    start = time.time()

    while t <= T:
        inflow = 'near(x[0],0)'
        outflow = 'near(x[0],1)'
        downwalls = 'near(x[1],0)'
        upwalls = 'near(x[1],1)'
        #inflow_profile = ('6.0*sin(pi*t/8)', '0')

        #bcu_inflow = DirichletBC(W.sub(0), Expression(inflow_profile, degree=2,t=t), inflow)
        #bcu_outflow = DirichletBC(W.sub(0), Expression(inflow_profile, degree=2,t=t), outflow)
        bcu_inflow = DirichletBC(W.sub(0), Constant((0.0, 0.0)), inflow)
        bcu_outflow = DirichletBC(W.sub(0), Constant((0.0, 0.0)), outflow)
        bcu_downwalls = DirichletBC(W.sub(0), Constant((1.0, 0.0)), downwalls) #下底
        bcu_upwalls = DirichletBC(W.sub(0), Constant((0.0, 0.0)), upwalls)

        bcB_inflow = DirichletBC(W.sub(2), Constant((0.0, 0.0)), inflow)
        bcB_downwalls = DirichletBC(W.sub(2), Constant((1.0, 0.0)), downwalls)
        bcB_upwalls = DirichletBC(W.sub(2), Constant((1.0, 0.0)), upwalls)
        bcB_outflow = DirichletBC(W.sub(2), Constant((0.0, 0.0)), outflow)

        bcr = DirichletBC(W.sub(3), Constant(0.0), boundary)

        bcT_inflow = DirichletBC(W.sub(4), Constant(1.0), inflow)
        bcT_outflow  = DirichletBC(W.sub(4), Constant(0.0), outflow )
        #bcp_outflow = DirichletBC(W.sub(1), Constant(0.0), outflow)

        bcs = [bcu_inflow,bcu_outflow, bcu_downwalls,bcu_upwalls, bcB_inflow, bcB_downwalls,bcB_upwalls,bcB_outflow,bcr,bcT_inflow,bcT_outflow]

        u, p, B, r, TTT = TrialFunctions(W)
        v, q, s, w, ttt = TestFunctions(W)
        (uk, pk, Bk, rk, TTTk) = Uk.split(deepcopy=True)

        # temperture-inner(beta*TTT,v)
        F = (dot((u - uk) / dt, v) + inner(grad(u), grad(v)) * Pru + 0.5 * inner(dot(grad(u), uk),v)- 0.5 * inner(dot(grad(v), uk), u) - p * div(v) - dot(scurl(B), scross(Bk, v)) * Sc + div(u) * q + 1.E-10 * inner(p, q)- inner(Pru*Ra*beta*TTT, v)) * dx\
            + (dot((B - Bk) / dt, s) + dot(scurl(B), scurl(s)) * Sc * Prb - dot(scross(u, Bk),scurl(s)) * Sc - inner(grad(r),s)+ inner(grad(w), B)) * dx\
            + (dot((TTT - TTTk) / dt, ttt) + inner(grad(TTT), grad(ttt)) * k + 0.5 * inner(dot(grad(TTT), uk),ttt) - 0.5 * inner(dot(grad(ttt), uk), TTT))* Prt * dx - (inner(f, v) + inner(g, s) + inner(ps, ttt)) * dx

        A0 = lhs(F)
        L0 = rhs(F)
        U = Function(W)

        solve(A0 == L0, U, bcs)
        (u, p, B, r, TTT) = U.split(deepcopy=True)
        plot(u, title='t=%.d Velocity' % t)
        plt.show()
        #plot(B, title='t=%.d B' % t)
        #plt.show()
        t += dt

    end = time.time()
    print("time：%.2f秒" % (end - start))

n = 30
E_L2_H1(n)

Thanks again!