Building mass and stiffness matrix for time dependent PDE with Gridap

Specific Domains Modelling & Simulations
,
1

Hi,

I have implemented an algorithm for the solution of system of ODEs that I want to test against a method of line discretization of a PDE.

My routine needs to take as input the mass matrix (assumed to be time-independent) and the right-hand side. Something that should be applied to:

M y’ = K y + f(t)
y[0] = y0

Looking at the Poisson tutorial I have tried doing something of the form:

u0f(x) = 2.0
ftrue(x,t) = 1.0 + 0.5*sin(2*pi*t)
f(x) = ftrue(x,0.0)
h(x) = 3.0
m(u,v) = ∫( v*u )*dΩ
a(u,v) = ∫( ∇(v)⋅∇(u) )*dΩ
b(v) = ∫( v*f )*dΩ + ∫( v*h )*dΓ
u0(v) = ∫( v*u0f )*dΩ
massop = AffineFEOperator(m,u0,Ug,V0)
op = AffineFEOperator(a,b,Ug,V0)

Mass = get_matrix(massop)
y0 = get_vector(massop)
L = -get_matrix(op)

function fx(t,ftrue,a,Ug,V0,h,dΩ,dΓ)
    f(x) = ftrue(x,t)
    b(v) = ∫( v*f )*dΩ + ∫( v*h )*dΓ
    op = AffineFEOperator(a,b,Ug,V0)
    return get_vector(op)
end

ff(t) = fx(t,ftrue,a,Ug,V0,h,dΩ,dΓ)

so that I can pass to the ODE solver the different pieces:

tspan = [0,1]
n = 100
t,y = odesolve(Mass,L,ff,tspan,n,y0)

The solution seems to numerically diverge. I have tried the same code against a finite difference discretization of an analogous problem, and things seem to work fine. I suppose I’m doing something wrong with the mass matrix/boundary assembly.

I have also tried looking at the

TransientLinearFEOperator( )

but those object do not seem to have matrix assembly routines.

Any pointers would be very helpful.

Thank you.