Time dependent problem with weak form depending on the solution of a secondary problem

variational formulation
1

Hi all,

I am solving a time dependent problem F1 where at each time-step I use its solution phi to solve a second problem F2. Then using the solution to this second problem x_sol I want to update the term dp in the weak form F1 and solve it at the next time step for these value. I am really having a hard time in order to somehow make this work. Can you please help me? For now I get the error at the weak form.

from fenics import *
import numpy as np

#  Function for solving the second problem:
def solve_second_problem(phi, t):
    global x_sol
    test_x= TestFunction(V2)
    x_trial = TrialFunction(V2)  

    F2 = + x_trial.dx(0)*test_x*dx\
         - stretch*(cos(psi))*test_x*dx 

    a = lhs(F2)
    L = rhs(F2)

    x = Function(V2)    
    solve(a == L, x, bcs2)

    x_sol = split(x)
    return [x_sol]
 

# Create mesh and set time step
L0=1
nx = 2**8 
mesh = IntervalMesh(2*nx,0,L0)
N  = 2**6     
dt  = 1/N   


# Function spaces for problem 1 and 2
P1      = FiniteElement('P',interval, 1)
element1 = MixedElement([P1, P1])
V1       = FunctionSpace(mesh,element1)
V2       = FunctionSpace(mesh, P1)

tol = 1E-14
# Boundary conditions  for problem 1
def boundary1(x, on_boundary):
    return on_boundary and near(x[0],0, tol) 
bc_psi = DirichletBC(V1.sub(0), Constant(np.pi/2), boundary1)
bcs1 = [bc_psi]

# Boundary conditions for problem 2
def boundary2(x, on_boundary):
    return on_boundary and near(x[0],0,tol)
bc_x = DirichletBC(V2,  Constant(0), boundary2)
bcs2 = [bc_x]


h  = Expression('h', h=dt, degree=1)
stretch =  Constant(2)


# Define functions for the problem 1
test_psi, test_u1S= TestFunctions(V1)
U_n = Function(V1) # Previous step unknowns
psi_n, u1S_n= split(U_n)
U = Function(V1) # current step unknowns
psi, u1S= split(U)
du = TrialFunction(V1)


t=0.0
i=0 

# Obtain an initial value for dp using phi=0 at t=0
solve_second_problem(0,0)

dp = x_sol # this is what I want to change at every time step, in F1 problem 

F1 = -test_psi.dx(0)*dx + stretch*(cos(psi)-sin(psi)*cos(psi))*test_psi*dx\
     +((u1S-u1S_n)+h*np.dot(dp,cos(psi)))*test_u1S*dx
J = derivative(F1,U,du)

# Initialization array of solutions
SolutionArray = []
SolutionArray.append(U) # append the most recent solution U
psi, u1S = U.split() # update variables

while t<50:

    t += dt    

    U_n.assign(SolutionArray[-1])    # update the previous step 

     
    problem = NonlinearVariationalProblem(F1, U, bcs1, J)
    solver = NonlinearVariationalSolver(problem)   
    solver.solve()   

    psi, u1S = U.split() # update variables
    SolutionArray = SolutionArray[1:]             # remove the first element
    SolutionArray.append(U)   # add the last solution
    
    solve_second_problem(psi,t)

    i += 1