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: * modify the model for linear elasticity. * modify the model for elastodynamics and run only a few steps (with bcs=[bc1, bc2] and with bcs=[bc1]). * change BCs to something similar to what you really use and run a few steps. * adjust mesh size/timestep/parameters so that they would approximately match what you use in your actual problem. * change mesh. 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<https://spam.ic.ac.uk/SpamConsole/Senders.aspx> 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<http://www.netl.doe.gov> [cid:97ebd409-8c9c-4361-9753-e4c626433151]
