how to add an exact, time-dependent right-hand-side vector for the method of manufactured solution (MMS)?

[aph4499]

Hello,

In the example below, I want to add an exact right-hand-side vector to my equations (exact_rhs). The exact right-hand-side does not depend on the variables solved for. It is given by say:
exact_rhs = u_exact(z)*u_exact(t)

The exact_rhs vector is time-dependent such that it needs to be updated at each time iteration

How would you go about adding that term to the equation below?

Thanks,

Anthony

# create grid
m = Grid1D(nx=len(z), Lx=zmax-zmin)

# create variables (system of equations)
c = CellVariable(name="c", mesh=m, hasOld=True, value=c0_vec)
...

# Exact right-hand-side vector for MMS: needs to be recomputed at each time-step
exact_rhs = u_exact(z)*u_exact(t)

# create equations
eq_c = TransientTerm(var=c) == DiffusionTerm(...) + ... + exact_rhs
...

full_eq = eq_c & ...

# Solve
for t in range(nnt-1):

    update exact_rhs
    c.updateOld()
    ...

    iter = 0
    res = 1.
    while ( res > 1e-7 and iter < 50 ):
        res = full_eq.sweep(dt=dt_soln)
        iter = iter + 1

    soln.c[t+1] = c
    ...

[guyer]

Make these changes

# offset the mesh so that it goes from zmin to zmax
# instead of 0 to (zmax - zmin)
m = Grid1D(nx=len(z), Lx=zmax-zmin) + [[zmin]]

# Declare a time variable
t = Variable(name="t", value=0.)

# Declare exact_rhs as a function of mesh coordinates and time
# It will automatically update when t changes
exact_rhs = u_exact(z=mesh.x, t=t)

# Solve
for step in range(nnt-1):

    # no need to update exact_rhs; it happens automatically
    c.updateOld()
    ...

    iter = 0
    res = 1.
    while ( res > 1e-7 and iter < 50 ):
        res = full_eq.sweep(dt=dt_soln)
        iter = iter + 1

    t.value = t.value + dt_soln
    soln.c[step+1] = c
    ...

[aph4499]

Hello,
Since my rhs are updated by calling a function, I integrated it in the following way (see below please) using CellVariables for the rhs. It does lead to the proper result. Is that acceptable?
Thanks,

# Fipy grid
mesh_num = Grid1D(nx=len(zgrid), Lx=zmax-zmin) + [[zmin]]

# Fipy dependent variables
Rxx  = CellVariable(name="Rxx", mesh=mesh_num, hasOld=True, value=Rxx0)
Rzz  = CellVariable(name="Rzz", mesh=mesh_num, hasOld=True, value=Rzz0)
az   = CellVariable(name="az", mesh=mesh_num, hasOld=True, value=az0)
b    = CellVariable(name="b", mesh=mesh_num, hasOld=True, value=cc0)
cbar = CellVariable(name="cbar", mesh=mesh_num, hasOld=True, value=c0)
ee   = CellVariable(name="ee", mesh=mesh_num, hasOld=True, value=eps0)

tvar = Variable(name="tvar", value=0.)

s_rhs_mms_rxx_  = CellVariable(mesh=mesh_num)
s_rhs_mms_rzz_  = CellVariable(mesh=mesh_num)
s_rhs_mms_az_   = CellVariable(mesh=mesh_num)
s_rhs_mms_cc_   = CellVariable(mesh=mesh_num)
s_rhs_mms_ee_   = CellVariable(mesh=mesh_num)
s_rhs_mms_cbar_ = CellVariable(mesh=mesh_num)

# System of equations (later will add coupling terms)
eq_cbar = TransientTerm(var=cbar) == (
    + DiffusionTerm(coeff=1/Re, var=cbar)
    + s_rhs_mms_cbar_
)

eq_Rxx = TransientTerm(var=Rxx) == (
    + DiffusionTerm(coeff=1/Re, var=Rxx)
    + s_rhs_mms_rxx_
)

eq_Rzz = TransientTerm(var=Rzz) == (
    + DiffusionTerm(coeff=1/Re, var=Rzz)
    + s_rhs_mms_rzz_
)

eq_ee = TransientTerm(var=ee) == (
    + DiffusionTerm(coeff=1/Re, var=ee)
    + s_rhs_mms_ee_
)

eq_az = TransientTerm(var=az) == (
    
    + DiffusionTerm(coeff=1/Re, var=az)
    + s_rhs_mms_az_
)

eq_b = TransientTerm(var=b) == (
    + DiffusionTerm(coeff=1/Re, var=b)
    + s_rhs_mms_cc_
)

# Fipy system of equations
transition = eq_cbar & eq_Rxx & eq_Rzz & eq_ee & eq_az & eq_b

# Viewer
vi = Viewer((Rxx, Rzz, az, b, ee))

for t in range(1, nnt-1):

    print("")
    print("t = %15.5f, (it = %i/%i)" % (ttimes[t], t, nnt) )
    print("=========================================")

    # Compute exact solution for comparison with numerical
    exact_sol.Rxx[t]  = Rxx0/ftime[0]*ftime[t]
    exact_sol.Rzz[t]  = Rzz0/ftime[0]*ftime[t]
    exact_sol.az[t]   = az0/ftime[0]*ftime[t]
    exact_sol.b[t]    = cc0/ftime[0]*ftime[t]
    exact_sol.ee[t]   = eps0/ftime[0]*ftime[t]
    exact_sol.cbar[t] = scipy.special.erf( zgrid/( 1 + coefdd*ttimes[t] ) )

    cbar.updateOld()
    Rxx.updateOld()
    Rzz.updateOld()
    ee.updateOld()
    az.updateOld()
    b.updateOld()

    # Update MMS rhs
    rhs_mms_rxx_, rhs_mms_rzz_, rhs_mms_az_, rhs_mms_cc_, rhs_mms_ee_, rhs_mms_cbar_ = \
        compute_rhs_terms_mms(zgrid,Re,Rxx0,Rzz0,az0,cc0,eps0,c0,t,ttimes,ftime,dftimedt,dRxx0dz,d2Rxx0dz2,dRzz0dz,d2Rzz0dz2,daz0dz,d2az0dz2,dcc0dz,d2cc0dz2,deps0dz,d2eps0dz2,dcbar0dz,d2cbar0dz2)

    s_rhs_mms_rxx_[0:]  = rhs_mms_rxx_
    s_rhs_mms_rzz_[0:]  = rhs_mms_rzz_
    s_rhs_mms_az_[0:]   = rhs_mms_az_
    s_rhs_mms_cc_[0:]   = rhs_mms_cc_ 
    s_rhs_mms_ee_[0:]   = rhs_mms_ee_
    s_rhs_mms_cbar_[0:] = rhs_mms_cbar_
    
    iter = 0
    res = 1.e50

    while ( res > 1e-4 and iter < 20 ):
    
        # Solve
        res = transition.sweep(dt=dt, solver=LinearLUSolver(tolerance = 1.e-10)) ### Direct solver does not work in parallel
        iter = iter + 1
        print("res = ", res)
        
    vi.plot()

    transoln.cbar[t] = cbar
    transoln.Rxx[t] = Rxx
    transoln.Rzz[t] = Rzz
    transoln.ee[t] = ee
    transoln.az[t] = az
    transoln.b[t] = b

[guyer]

If you're getting the answer you want, I guess it's fine, but you're doing too much work

    # Update MMS rhs
    rhs_mms_rxx_, rhs_mms_rzz_, rhs_mms_az_, rhs_mms_cc_, rhs_mms_ee_, rhs_mms_cbar_ = \
        compute_rhs_terms_mms(zgrid,Re,Rxx0,Rzz0,az0,cc0,eps0,c0,t,ttimes,ftime,dftimedt,dRxx0dz,d2Rxx0dz2,dRzz0dz,d2Rzz0dz2,daz0dz,d2az0dz2,dcc0dz,d2cc0dz2,deps0dz,d2eps0dz2,dcbar0dz,d2cbar0dz2)

    s_rhs_mms_rxx_[0:]  = rhs_mms_rxx_
    s_rhs_mms_rzz_[0:]  = rhs_mms_rzz_
    s_rhs_mms_az_[0:]   = rhs_mms_az_
    s_rhs_mms_cc_[0:]   = rhs_mms_cc_ 
    s_rhs_mms_ee_[0:]   = rhs_mms_ee_
    s_rhs_mms_cbar_[0:] = rhs_mms_cbar_

I can't tell you how to fix this without knowing what compute_rhs_terms_mms() looks like, but there's almost certainly a simpler way to do this.

[aph4499]

Hello,
This is what is inside compute_rhs_terms_mms (see below please), with for instance (zcoord is my grid in z):

dRzz0dz = -2*zcoord*exp_rzz/sig_rzz**2.
dRzz0dz = dRzz0dz*Arzz*ftime[0]

d2Rzz0dz2 = 4*zcoord**2.*exp_rzz/sig_rzz**4. - 2*exp_rzz/sig_rzz**2.
d2Rzz0dz2 = d2Rzz0dz2*Arzz*ftime[0]

with

exp_rzz = exp(-( zcoord/sig_rzz )**2.), sig_rzz = constant, Arzz = constant

and

Rzz0 = Arzz * numpy.exp(-(zcoord/sig_rzz)**2.) * ftime[0]

and the time function is defined as:

coeft = 1
ftime = numpy.exp(coeft*ttimes)
dftimedt = coeft*numpy.exp(coeft*ttimes)

Thanks

=================================================

def compute_rhs_terms_mms(zcoord,Re,Rxx0,Rzz0,az0,cc0,eps0,c0,i,ttimes,ftime,dftimedt,dRxx0dz,d2Rxx0dz2,dRzz0dz,d2Rzz0dz2,daz0dz,d2az0dz2,dcc0dz,d2cc0dz2,deps0dz,d2eps0dz2,dcbar0dz,d2cbar0dz2):

    # Compute rhs terms                                                                                                                                                                                                                                    
    dRxxdz_exact_    = dRxx0dz/ftime[0]*ftime[i]
    d2Rxxdz2_exact_  = d2Rxx0dz2/ftime[0]*ftime[i]
    d2Rzzdz2_exact_  = d2Rzz0dz2/ftime[0]*ftime[i]
    d2azdz2_exact_   = d2az0dz2/ftime[0]*ftime[i]
    d2ccdz2_exact_   = d2cc0dz2/ftime[0]*ftime[i]
    d2epsdz2_exact_  = d2eps0dz2/ftime[0]*ftime[i]

    dcbardz_exact_   = 2*numpy.exp(-zcoord**2./( 1 + coefdd*ttimes[i] )**2.)/( numpy.sqrt(numpy.pi)*( 1 + coefdd*ttimes[i] ) )
    d2cbardz2_exact_ = -4*zcoord*numpy.exp(-zcoord**2./( 1 + coefdd*ttimes[i] )**2.)/( numpy.sqrt(numpy.pi)*( 1 + coefdd*ttimes[i] )**3. )

    dcbardt       = -2*coefdd*zcoord*numpy.exp(-zcoord**2./( 1 + coefdd*ttimes[i] )**2.)/( numpy.sqrt(numpy.pi)*( 1 + coefdd*ttimes[i] )**2. )

    rhs_mms_rxx_  = Rxx0/ftime[0]*dftimedt[i] - 1/Re*d2Rxxdz2_exact_
    rhs_mms_rzz_  = Rzz0/ftime[0]*dftimedt[i] - 1/Re*d2Rzzdz2_exact_
    rhs_mms_az_   = az0/ftime[0]*dftimedt[i] - 1/Re*d2azdz2_exact_
    rhs_mms_cc_   = cc0/ftime[0]*dftimedt[i] - 1/Re*d2ccdz2_exact_
    rhs_mms_ee_   = eps0/ftime[0]*dftimedt[i] - 1/Re*d2epsdz2_exact_
    rhs_mms_cbar_ = dcbardt - 1/Re*d2cbardz2_exact_

[guyer]

My recommendation is

# Fipy grid
mesh_num = Grid1D(nx=len(zgrid), Lx=zmax-zmin) + [[zmin]]

zcoord = mesh_num.x

tvar = Variable(name="tvar", value=0.)


coeft = 1
ftime0 = numerix.exp(coeft*ttimes[0])
ftime = numerix.exp(coeft*tvar)
dftimedt = coeft*numerix.exp(coeft*tvar)

Rzz0 = Arzz * numerix.exp(-(zcoord/sig_rzz)**2.) * ftime0

# Fipy dependent variables

...

exp_rzz = numerix.exp(-( zcoord/sig_rzz )**2.)
dRzz0dz = -2*zcoord*exp_rzz/sig_rzz**2.
dRzz0dz = dRzz0dz*Arzz*ftime0

d2Rzz0dz2 = 4*zcoord**2.*exp_rzz/sig_rzz**4. - 2*exp_rzz/sig_rzz**2.
d2Rzz0dz2 = d2Rzz0dz2*Arzz*ftime0

dRxxdz_exact_    = dRxx0dz/ftime0*ftime
d2Rxxdz2_exact_  = d2Rxx0dz2/ftime0*ftime
d2Rzzdz2_exact_  = d2Rzz0dz2/ftime0*ftime
d2azdz2_exact_   = d2az0dz2/ftime0*ftime
d2ccdz2_exact_   = d2cc0dz2/ftime0*ftime
d2epsdz2_exact_  = d2eps0dz2/ftime0*ftime

dcbardz_exact_   = 2*numerix.exp(-zcoord**2./( 1 + coefdd*tvar )**2.)/( numerix.sqrt(numerix.pi)*( 1 + coefdd*tvar ) )
d2cbardz2_exact_ = -4*zcoord*numerix.exp(-zcoord**2./( 1 + coefdd*tvar )**2.)/( numerix.sqrt(numerix.pi)*( 1 + coefdd*tvar )**3. )

s_rhs_mms_rxx_  = Rxx0/ftime0*dftimedt - 1/Re*d2Rxxdz2_exact_
s_rhs_mms_rzz_  = Rzz0/ftime0*dftimedt - 1/Re*d2Rzzdz2_exact_
s_rhs_mms_az_   = az0/ftime0*dftimedt - 1/Re*d2azdz2_exact_
s_rhs_mms_cc_   = cc0/ftime0*dftimedt - 1/Re*d2ccdz2_exact_
s_rhs_mms_ee_   = eps0/ftime0*dftimedt - 1/Re*d2epsdz2_exact_
s_rhs_mms_cbar_ = Cdcbardt - 1/Re*d2cbardz2_exact_

# System of equations (later will add coupling terms)

...

for t in range(1, nnt-1):

    tvar.value = ttimes[t]

    print("")
    print("t = %15.5f, (it = %i/%i)" % (tvar, t, nnt) )
    print("=========================================")

    # Compute exact solution for comparison with numerical
   ...
 
    cbar.updateOld()
    Rxx.updateOld()
    Rzz.updateOld()
    ee.updateOld()
    az.updateOld()
    b.updateOld()

    iter = 0
    res = 1.e50

    while ( res > 1e-4 and iter < 20 ):
    
        # Solve
        res = transition.sweep(dt=dt, solver=LinearLUSolver(tolerance = 1.e-10)) ### Direct solver does not work in parallel
        iter = iter + 1
        print("res = ", res)
        
    vi.plot()

    ...

You can have arrays ttime's and ftime's if you want, but it's important that these terms get defined in terms of a tvar and an ftime that are FiPy Variable objects. That way, everything automatically updates whenever you call tvar.value = ttimes[t].