Trouble in extending a simple diffusion problem in time

General
1

I am solving a simple diffusion problem with Neumann boundary conditions in FEniCS. I am trying to extend the code in time, i.e., I ran it for time t =0 to t=0.5 and saved the solution in a file. Now I want to read the solution at t = 0.5 and pass it as an initial condition for the run from t=0.5 to t=1.

Here, the solver does a funny thing. When I print the value of the field upto machine precision for the extended run, it reads the initial condition correctly from the file but the field at the next timestep is off beyond the 12th decimal place. I am comparing the solution in the extended run with the one obtained by a direct run from t = 0 to t = 1.

I do not understand the mistake I’m making. Worryingly, I have to do this extension business for a non-linear problem where I’m interested in the long-time limit of the solution. I do not know how the errors would propagate for the non-linear problem, and how trustworthy the solution would be in the long-time limit.

Posting here the code where I have separated the run part from the solving of the differential equation part. The parameters are taken from the jsonfile.

SOLVE CODE:

import dolfin as df
import progressbar
import os
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
np.set_printoptions(precision=16)

class diffusion(object):
    def __init__(self, parameters):
        # read in parameters
        for key in parameters:
            setattr(self, key, parameters[key])
        
        self.mesh = df.IntervalMesh(self.resolution, 
                                    -self.L/2, 
                                     self.L/2)

        self.function_space = df.FunctionSpace(self.mesh,'P', 1)
        self.c = df.Function(self.function_space)
        self.tc = df.TestFunction(self.function_space)
        self.c0 = df.Function(self.function_space)
        
    def setup_initial_conditions(self, initc):
        if initc is None:
            self.c0.interpolate(df.Expression('1 + 0.2 * cos(2 * pi*x[0]/L)', pi=np.pi, L=self.L, degree=1))
        else:
            self.c0.interpolate(initc)
        print('c0=', self.c0.compute_vertex_values(self.mesh))

    def setup_weak_forms(self):
        self.form = (df.inner((self.c - self.c0)/self.timestep, self.tc) 
            + self.D * df.inner(df.nabla_grad(self.c), df.nabla_grad(self.tc))
            ) * df.dx

    def solve(self, initc=None, initTime=0.0, extend=False):

        # create file if it doesn't exist, open if it does
        self.cFile   = df.XDMFFile('%s_concentration.xdmf' % self.timestamp)
        
        # time-variables
        self.time = initTime
        print('time=', self.time)
        maxsteps = int(self.maxtime/self.timestep)

        self.setup_initial_conditions(initc)
        self.setup_weak_forms()

        # if fresh run, write initial condition
        if not extend:
            self.cFile.write_checkpoint(self.c0, 'concentration', self.time)

        for steps in progressbar.progressbar(range(1, maxsteps+1)):
            self.time += self.timestep
            print('time=', self.time)    
            df.solve(self.form == 0, self.c)
            self.cFile.write_checkpoint(self.c, 'concentration', 
                                            self.time, append=True)
            print('c=', self.c.compute_vertex_values(self.mesh))
            self.c0.assign(self.c)
            

        self.cFile.close()

RUN CODE:

import os
from extend_for_diffusion import diffusion
import  argparse
import dolfin as df
import h5py

parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='data file',
        default='parameters.json')
parser.add_argument('-t','--time', help='time to run', type=float)
args = parser.parse_args()
assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile

# load the parameters
with open(args.jsonfile) as jsonFile:                    
    parameters = json.load(jsonFile)

if args.jsonfile=='parameters.json':
    extend = False
    print('Fresh run...')
    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
    parameters['timestamp'] = timestamp
    parametersFile = timestamp + '_parameters.json'
    initc = None
    initTime = 0.0

else:
    extend = True
    print('Extending run...')
    parametersFile = args.jsonfile

    steps = int(parameters['maxtime']/parameters['timestep'])
    
    mesh = df.IntervalMesh(parameters['resolution'], 
                            - parameters['L']/2, 
                              parameters['L']/2)
    # Read data
    var = 'concentration'
    h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var)) 
    SFS = df.FunctionSpace(mesh, 'P', 1)
    initc = df.Function(SFS)
    cFile   = df.XDMFFile('%s_concentration.xdmf' % parameters['timestamp'])
    cFile.read_checkpoint(initc, 'concentration', steps)
    cFile.close()
    initTime = parameters['maxtime']
    parameters['maxtime'] = args.time

diff = diffusion(parameters)
diff.solve(initc, initTime, extend)
if extend:
    parameters['maxtime'] = initTime + parameters['maxtime']
with open(parametersFile, "w") as jsonFile:
    json.dump(parameters, jsonFile, indent=4, sort_keys=True)

PARAMETERS JSONFILE:

    "resolution": 10,
    "L": 1.0, 
    "D": 0.1,
    "timestep": 0.1,
    "maxtime": 0.5
}
2

Please format properly your code. It should be in a single cell with three backticks, i.e “```”, otherwise it will be hard for anyone to copy your code.

3

My apologies. Edited it now.