I have done it, but source function cant be recalculate each time step. How can i do it?
from dolfin import*
import numpy as np
Nt=1.0E+20
k=5.0E-22
p=1.E-7
dt = 1.8E+2
Tf = dt*20
t = float(dt)
arq_saida=“flux-Monophase”
o1 = open(arq_saida,‘w’)
Lcx=10.E+3
Lc=0.55E+3
mesh=RectangleMesh(Point(0., 0.), Point(Lc, Lcx), 50, 50, “left”)
boundaries = MeshFunction(“size_t”, mesh, mesh.topology().dim()-1)
V = FunctionSpace(mesh, “CG”, 1)
coordinates = V.tabulate_dof_coordinates()
dof_x = coordinates[:, 0]
X_max=np.amax(dof_x)
X_min=np.amin(dof_x)
Lc=X_max - X_min
class Rigth(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0] - X_max) < tol
class Left(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0]- X_min) < tol
saida=Rigth()
saida.mark(boundaries,1)
Entrada=Left()
Entrada.mark(boundaries,2)
dS = dS(subdomain_data=boundaries)
ds = ds(subdomain_data=boundaries)
D=3.9E-0
u0 = Constant(0.0)
bc1 = DirichletBC(V, u0, boundaries,1)
bcs = [bc1]
#Variational problem at each time
u0= Function(V)
u= Function(V)
f1 = Function(V)
f2 = Function(V)
v= TestFunction(V)
ut=Function(V)
n0=2.0
g2=-1.0E-9
du = TrialFunction(V)
f2=interpolate(Expression(“exp(g2 * pow(ut,n0))”,ut=ut, n0=n0, g2=g2, degree=2),V)
f0=exp(g2 * pow(ut,n0))
f1=Nt * k * exp(g2 * pow(ut,n0)) *u-p * ut
J0=3.0E+6
theta=0.5
dsN = Measure(“ds”, subdomain_id=2, subdomain_data=boundaries)
ds_out = Measure(“ds”, subdomain_id=1, subdomain_data=boundaries)
def Ftheta(u,v):
return D * inner(grad(u), grad(v)) * dx-J0 * v * dsN + f1 * v * dx\
F = (1./dt)(u-u0)vdx + thetaFtheta(u,v)
+(1.-theta)*Ftheta(u0,v)
J = derivative(F,u,du)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
fileu = File(“Monophase/u.pvd”)
fileut = File(“Monophase/ut.pvd”)
fileuf = File(“Monophase/f2.pvd”)
flu1=
tempo=
n1 = FacetNormal(mesh)
contador=0
while t <= Tf+ DOLFIN_EPS:
solver.solve()
f2=interpolate(Expression(“exp(g2pow(ut,n0))",ut=ut, n0=n0, g2=g2, degree=2),V)
ut.vector()[:] = dtNtkf2.vector()u.vector()-dtput.vector()+ut.vector()
u0.assign(u)
contador+=1
if contador%1==0:
fileu << (u,t)
fileut << (ut,t)
fileuf << (f2,t)
flux = -dot(grad(Du), n1)*ds_out
flux = assemble(flux)
if contador%1==0:
o1.write(”%s"%t)
o1.write(" %s \n" %flux)
t += float(dt)
print (“tempo = %s , contador = %s , fluxo = %s” %(t, contador, flux))
o1.close()