get wrong output

[gfkdliucheng]

from devito import import matplotlib.pyplot as plt def solver(I, dt, T, eqs): dt = float(dt) Nt = int(round(T/dt)) t = TimeDimension('t', spacing=Constant('h_t')) u = TimeFunction(name='u', dimensions=(t,), shape=(Nt+1,), space_order=2, time_order=2) u.data[:] = I eqn = eval(eqs) stencil = Eq(u.forward, solve(eqn, u.forward)) op = Operator(stencil) op.apply(h_t=dt, t_M=Nt-1) return u.data, np.linspace(0, Ntdt, Nt+1)

u,t = solver(1, 0.01, 1,r'u.dt-t') print(u[-1]) # wrong output：50.499958

def solver1(I, dt, T):

dt = float(dt)
Nt = int(round(T/dt))
u = np.zeros(Nt+1)
t = np.linspace(0, Nt*dt, Nt+1)
u[0] = I
for n in range(Nt):
    u[n+1] = u[n] + dt*t[n]
return u, t

u,t = solver1(1, 0.01, 1) print(u[-1]) # correct output: 1.494999

[mloubout]

u[n+1] = u[n] + dt*t[n]

that's " u(t)' - t = 0" this has nothing to do with the other equation so not sure why comparing to this otr why you'd expect it to be the same.

"u'' + w2u = 0 " that's `u[n+1] = - u[n-1] + (2 + dt2 * w*2) u[n]`

[gfkdliucheng]

" u(t)' - t = 0" is my question, but how to implement this with devito, def solver(I, dt, T, eqs):
dt = float(dt) Nt = int(round(T/dt)) t = TimeDimension('t', spacing=Constant('h_t')) u = TimeFunction(name='u', dimensions=(t,), shape=(Nt+1,), space_order=2, time_order=2) u.data[:] = I eqn = eval(eqs) stencil = Eq(u.forward, solve(eqn, u.forward)) op = Operator(stencil) op.apply(h_t=dt, t_M=Nt-1) return u.data, np.linspace(0, Nt*dt, Nt+1)

              u,t = solver(1, 0.01, 1,r'u.dt-t')
              print(u[-1]) # wrong output：50.499958

[mloubout]

First order equation so you need

u = TimeFunction(name='u', dimensions=(t,), shape=(Nt+1,), space_order=2, time_order=1)

then r'u.dt-t' is the time index you need r'u.dt-t*t.spacing') to get the actual time

For future reference, please use triple quotes "```" to make code blocks so it's easier to read and test

In [8]: def solver(I, dt, T, eqs):
   ...:     dt = float(dt)
   ...:     Nt = int(round(T/dt))
   ...:     t = TimeDimension('t', spacing=Constant('h_t'))
   ...:     u = TimeFunction(name='u', dimensions=(t,),
   ...:     shape=(Nt+1,), space_order=2, time_order=1)
   ...:     u.data[:] = I
   ...:     eqn = eval(eqs)
   ...:     stencil = Eq(u.forward, solve(eqn, u.forward))
   ...:     print(stencil)
   ...:     op = Operator(stencil)
   ...:     op.apply(h_t=dt, t_M=Nt-1)
   ...:     return u.data, np.linspace(0, Nt*dt, Nt+1)
   ...: 
   ...: u,t = solver(1, 0.01, 1,r'u.dt-t*t.spacing')
   ...: print(u[-1]) # wrong output：50.499958
Eq(u(t + h_t), h_t*(t*h_t + u(t)/h_t))
Operator `Kernel` ran in 0.01 s
1.494997

[gfkdliucheng]

Thank you, bro! God bless you, but how（where） can i get the information “using t*t.spacing to get the actual time”