0

我目前正在开展一个项目，我们正在解决 FEniC 中的 PDE 系统。我创建了以下代码以解决系统问题，但出现无效语法错误

a = a0 + a1

我在 Python 方面不是那么好，而且我以前从未使用过 FEniCS。我还使用 Windows 子系统来运行它，这让我更难以理解我可能犯的任何错误。感谢您提出的任何建议，如果我提出明显的问题，我会提前道歉！

from fenics import *
 
# Create mesh and define function space
mesh = Mesh (" circle.xml ")

# Construct the finite element space
V = VectorFunctionSpace (mesh , 'P', 1)

# Define parameters :
T = 150
dt = 0.5
alpha = 0.4
beta = 2
gamma = 0.8
delta = 1

# Class representing the intial conditions
class InitialConditions ( UserExpression ):
    def eval (self , values , x):
        values [0] = Expression(("(4/25)-2*pow(10,-7)*(x[0]-0.1*x[1]-225)*(x[0]-0.1*x[1]-675)"),degree=2)
        values [1] = Expression(("(22/45)-3*pow(10,-5)*(x[0]-450)-1.2*pow(10,-4)*(x[1]-150)"),degree=2)

def value_shape ( self ):
    return (2 ,)

# Define initial condition
indata = InitialConditions(degree =2)
u0 = Function (V)
u0 = interpolate (indata,V)

# Test and trial functions
u,v = TrialFunction(V), TestFunction(V)

# Create bilinear and linear forms
a0 = (u[0]*v[0]*dx) + (0.5*delta*dt*inner(grad(u[0]),grad(v[0]))*dx)
a1 = (u[1]*v[1]*dx) + (0.5*delta*dt*inner(grad(u[1]),grad(v[1]))*dx)

L0 = (u0[0]*v[0]*dx) - (0.5*delta*dt*inner(grad(u0[0]),grad(v[0]))*dx) - (dt*u0[0]*v[0]*dx*(((u0[0]*u0[1])/(u0[0]+alpha))-u0[0]*(1-u0[0]))*dx) 
L1 = (u0[1]*v[1]*dx) - (0.5*delta*dt*inner(grad(u0[1]),grad(v[1]))*dx) - (dt*u0[1]*v[1]*dx*(-beta*((u0[0]*u0[1])/(u0[0]+alpha))-gamma*u0[1]*dx) 

a = a0 + a1
L = L0 + L1

#Set up boundary condition
g = Constant([0.0,0.0])
bc = DirichletBC(V,u_initial,DirichletBoundary())
bc = [] #NEUMANN

#Assemble matrix
A = assemble(a)

# Set an output file
out_file = File("Results.pvd","compressed")

# Set initial condition

u = Function(V)
u.assign(u0)

t = 0.0

out_file << (u,t)

u_initial = Function(V)
u_initial.assign(u0)

t_save = 0
num_samples = 20

# Time - stepping
while t < T:
                                                                          
    # assign u0
    u0.assign(u)
    
    #Assemble vector and apply boundary conditions
    A = assemble(a)
    b = asseble(L)

    t_save += dt

    if t_save > T/ num_samples or t >= T-dt:
        print("Saving!")
        
        #Save the solution to file
        out_file << (uv,t)
    
        t_save = 0
    
    #Move to next interval and adjust boundary condition
    t += dt
4

1 回答 1

0

您的行中有一个错字b = asseble(L)-->b=assemble(L)

也许这个小错误给你带来了问题？虽然我想错误信息会更具描述性。

于 2021-02-11T19:12:18.073 回答