Hi Abhishek,

I think your approach is fine and is supposed to work. The following code to some extent mimics what you are trying to do and it works fine.  I would gradually modify this small code and observe what specific piece of your code is actually causing the issue. I would do that in the following order:

Thanks,
Koki

```
from firedrake import *
 

mesh = UnitCubeMesh(2, 2, 2)
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
u = TrialFunction(V)

a = inner(grad(u), grad(v)) * dx
L = Constant(0) * v * dx

bc1 = DirichletBC(V, 0.0, 1)
bc2 = DirichletBC(V, 1.0, 2)

sol = Function(V)
problem1 = LinearVariationalProblem(a, L, sol, bcs=[bc1, bc2], constant_jacobian=True)
solver1 = LinearVariationalSolver(problem1, solver_parameters={"ksp_type": "cg", "pc_type": "jacobi"})
problem2 = LinearVariationalProblem(a, L, sol, bcs=[bc1], constant_jacobian=True)
solver2 = LinearVariationalSolver(problem2, solver_parameters={"ksp_type": "cg", "pc_type": "jacobi"})

solver1.solve()
print("after 1st solve:\n", sol.dat.data_ro)
solver2.solve()
print("after 2nd solve:\n", sol.dat.data_ro)
```

Outputs:
```
after 1st solve:
 [0.  0.  0.  0.5 0.  0.5 0.  0.  0.5 0.  0.  0.5 0.5 0.5 1.  0.5 0.5 0.
 1.  1.  1.  1.  1.  0.5 1.  1.  1. ]
after 2nd solve:
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 3.88578059e-16
 0.00000000e+00 5.55111512e-16 0.00000000e+00 0.00000000e+00
 5.55111512e-16 0.00000000e+00 0.00000000e+00 3.33066907e-16
 3.88578059e-16 3.88578059e-16 2.22044605e-16 3.33066907e-16
 3.33066907e-16 0.00000000e+00 4.44089210e-16 4.44089210e-16
 5.55111512e-16 5.55111512e-16 3.33066907e-16 3.33066907e-16
 4.44089210e-16 4.44089210e-16 3.33066907e-16]
```

From: Venketeswaran, Abhishek (FN) <Abhishek.Venketeswaran@netl.doe.gov>
Sent: Thursday, January 27, 2022 11:57 PM
To: firedrake <firedrake@imperial.ac.uk>; Sagiyama, Koki <k.sagiyama@imperial.ac.uk>
Subject: Help with implementing steps in simulation
 

This email from Abhishek.Venketeswaran@netl.doe.gov originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders list to disable email stamping for this address.

 

Hi,
I would like to prescribe a boundary condition for an initial time duration and then remove it for the remaining time duration. Commercial finite element solvers allow you to do this by specifying steps. I tried to implement it using the code below. This is giving nonsensical results at and after the time at which the BCs change. Can you please assist me in how I may implement this properly. I can share the full code if needed.

m = rho*inner((u - 2.0 * u_n + u_nm1),v) / Constant(dt * dt) * dxlump
k= inner(sigma(u_n),epsilon(v))*dx

def collar_load(t,T,freq,amp=1.0):
    return amp*hanning(t,T)*stepfn(t,T)*np.sin(freq*t)

# Define the loading as an expression depending on t
t=0
collar_disp= Constant([0.0, 0.0, 0.0])
F = m+k
a, r = lhs(F), rhs(F)
zero = Constant((0.0, 0.0, 0.0))
# Dirichlet BC that is valid for the entire time duration
bc1 = DirichletBC(V, zero, 11)
# Collar BC that remains on till T_pulse
bc2= DirichletBC(V, collar_disp, 21)
# My attempt of defining two problems to mimic two simulation steps
problem1 = LinearVariationalProblem(a, r, u_np1, bcs=[bc1,bc2], constant_jacobian=True)
solver1 = LinearVariationalSolver(problem1, solver_parameters=params)
problem2 = LinearVariationalProblem(a, r, u_np1, bcs=[bc1], constant_jacobian=True)
solver2 = LinearVariationalSolver(problem2, solver_parameters=params)

for i in range(len(time)):
    t = time[i]
    if (t<=T_pulse):
        solver1.solve()
    else:
        solver2.solve()
     u_nm1.assign(u_n)
     u_n.assign(u_np1)

--
Abhishek Venketeswaran
NETL Research Associate - ORISE
National Energy Technology Laboratory
Department of Energy
abhishek.venketeswaran@netl.doe.gov
Work: 412-386-4833
Mobile: 716-507-7890
www.netl.doe.gov