Hello.
I am struggling with a problem, where I want to compute distribution of some field J_s(x, y) when I know the distribution of field A(x, y) on the subdomain with integral constraint equal to I_t(t) at a given time step. The A and J_s can be either vector or scalar valued. I have chosen scalar for 2D example.
This can be expressed as: \int\int_{\Omega_1} (- \sigma \partial_t A + J_s) dV_1 = I_t(t) applied to subdomain \Omega_1 at a given time step.
The PDE that is solved is: (\nabla\times\frac{1}{\mu}\nabla\times A)\vert_{\Omega} = (J_s - \sigma \partial_t A)\vert_{\Omega_1}, the left side is solved in whole domain whereas right only in subdomain \Omega_1.
The first equation operates as an integral constraint and spits out J_s, both equations have to be solved simultaneously and I_t(t) is time-dependent.
I have approached this by defining P_1 and Real elements and combining those to from P_1 * R mixed function space. The test functions are v and r accordingly. The weak formulation I came up with is below:
(\nabla v, \frac{1}{\mu} \nabla A^{n})\vert_{\Omega} = (v, J_s)\vert_{\Omega_1} - (v, \frac{\sigma}{\Delta t} (A^{n} - A^{n-1})) - (r, J_t - J_s + \frac{\sigma}{\Delta t} (A^{n} - A^{n-1}))\vert_{\Omega_1}
The last term (r, J_t - J_s + \sigma (A^{n} - A^{n-1})/\Delta t)\vert_{\Omega_1} seems to be computed in whole domain \Omega instead of only in \Omega_1 with a code included below.
I don’t know if what I am doing is right or wrong and I don’t fully understand how Real elements work. To my understanding it substitutes the whole space with a single element and sums fields. That is why I substituted I_t which is constraint with J_t which would be a density of that field.
I have a feeling that this is similar type of problem as reaction-diffusion (there were couple questions about those) where the total amount of substances should be constant, but unfortunately it’s not my field of expertise. And over here, the amount of field should be constant in a certain subdomain only.
Could anybody help? How should I approach this?
The yellow circle is \Omega_1.
I’ve moved everything to one side and used nonlinear solver simply not to mess with moving A^{n-1} and A^n around.
import dolfin as df
import mshr
import matplotlib.pyplot as plt
import numpy as np
mm = 1e-3
radius = 3*mm
mu = 4*np.pi*1e-7
sigma = 58e6
periods = 1
f = 50e3
T = periods/f
nSteps = 100*periods
dt = T/nSteps
# create geometry
region = mshr.Rectangle(df.Point(-10*mm,-10*mm), df.Point(10*mm,10*mm))
wire = mshr.Circle(df.Point(0*mm,0*mm), radius)
region.set_subdomain(1, wire)
# mesh and subdomains
mesh = mshr.generate_mesh(region, 64)
markers = df.MeshFunction('size_t', mesh, 2, mesh.domains())
dx = df.Measure('dx', domain=mesh, subdomain_data=markers)
# function space
P1 = df.FiniteElement('P', mesh.ufl_cell(), 1)
R = df.FiniteElement('Real', mesh.ufl_cell(), 0)
ME = df.FunctionSpace(mesh, P1*R)
AnJn = df.TrialFunction(ME)
A_n, J_s = df.split(AnJn) # A_n = P1, J_s = R
UV = df.TestFunction(ME)
u, r = df.split(UV) # u = P1, r = R
An1Jn = df.Function(ME) # previous
A_n_1, _ = df.split(An1Jn)
# boundary conditions
bc = df.DirichletBC(ME.sub(0), df.Constant(0.0), 'on_boundary')
J_tot = df.Expression('J_max', J_max = 0.0, degree=0) # total current
# ----------------------------------------------------------------
F = df.inner(dt/mu*df.grad(A_n), df.grad(u))*dx \
+ sigma*u*(A_n-A_n_1)*dx(1) \
- (dt*u*J_s*dx(1) \
+ dt*r*(J_tot - J_s + sigma/dt*(A_n-A_n_1))*dx(1))
AnJn = df.Function(ME)
F = df.action(F, AnJn)
J = df.derivative(F, AnJn)
problem = df.NonlinearVariationalProblem(F, AnJn, bc, J)
solver = df.NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1.24E-7
prm['newton_solver']['relative_tolerance'] = 1E-8
prm['newton_solver']['maximum_iterations'] = 25
prm['newton_solver']['relaxation_parameter'] = 1.0
t = 0.0
pvdout = df.File('Out.pvd')
for i in range(nSteps):
solver.solve()
An1Jn.assign(AnJn)
pvdout << (AnJn.split()[1], t) # return Js
t += dt
J_tot.J_max = 2e6*np.sin(2*np.pi*f*t)