<jupyter_start><jupyter_text>Solving for thermal advection and diffusion in 2D using the Fenics library<NAME>Abstract:This is an enxample note book to explore if we can make a submission process for a subsurface online reproducable journal. This notebook has been written to solve classic advection diffusion. The notebook uses for packages: fenics, mshr, and numpy.Introduction:This is an example of using the fenics library to solve for the thermal advection and diffusion in two dimensions. The equations are solved using the finite element technique within a cartesian box. The velocity of the material within the box is predefined ans a linear function of space. The aim of this notebook is simply to explore how to deal with such a notebook and how to test it for reproducability.<jupyter_code>%matplotlib inline
# I need to list all the funtions I use rather than a *
from fenics import Constant,Point,plot,FunctionSpace,VectorFunctionSpace,Expression,interpolate,SubDomain, \
                   DirichletBC,near,TrialFunction,TestFunction,Expression,dx,dot,grad,lhs,rhs,solve,Function
from mshr import Rectangle,generate_mesh
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>conda install -c conda-forge fenicsconda install -c conda-forge mshr I will set up the domain dimensions for the rectanular box and the model resolution.<jupyter_code># Domain dimensions in meters
lx = 1.5e5
lz = 1.0e5
# number of elements along z-edge
numz = 32<jupyter_output><empty_output><jupyter_text>Here I define the physical parameters that are used within the model and the time stepping for the numerical solution.<jupyter_code># Physical parameters
year      = 365*24*60*60
K         = 1e-6  # thermal diffusivity

# Time steps
num_steps = 20       # number of time steps
dtnom     = 1e5*year   # time step size (in years)
dtc       = Constant(dtnom)<jupyter_output><empty_output><jupyter_text>Now I am ready to define and solve the set of equations. The first step is to create the mesh using the library mshr and display the mesh.<jupyter_code># Create mesh and define function space
domain = Rectangle(Point(0,0), Point(lx,lz))
mesh = generate_mesh(domain,numz)
p = plot(mesh)
plt.title("Mesh")
plt.show()<jupyter_output><empty_output><jupyter_text>Then define the functions in which temperature and velocity will be stored.<jupyter_code># Define function space for temperature
V = FunctionSpace(mesh,'CG',1)
# Define function space for velocity
W = VectorFunctionSpace(mesh,'P',2)<jupyter_output><empty_output><jupyter_text>Define the intial conditions<jupyter_code>T_D = Expression('x[1]>=lz-10e3 ? 1450 : 0',lz=lz,degree=2)
T_n = interpolate(T_D, V)<jupyter_output><empty_output><jupyter_text>Create the boundary conditions<jupyter_code>class East(SubDomain):
    def inside(self, x , on_boundary):
       return near(x[0], lx)

class West(SubDomain):
    def inside(self, x , on_boundary):
       return near(x[0], 0.0)

class North(SubDomain):
    def inside(self, x , on_boundary):
       return near(x[1], lz)

class South(SubDomain):
    def inside(self, x , on_boundary):
       return near(x[1], 0.0)

Tbc = DirichletBC(V,T_D,North())<jupyter_output><empty_output><jupyter_text>We then define the test and trial functions.<jupyter_code># Define variational problem

#T Temperature
T = TrialFunction(V)
vt = TestFunction(V)<jupyter_output><empty_output><jupyter_text>The velocity of the material is defined as linarly increasing in the x direction and decreasing in the positive z direction.<jupyter_code># solid flow
u = Expression(('1e-2/year/lz*(lz-x[1])','-1e-2/year/lz*(lx-x[0])'),year=year,lz=lz,lx=lx,degree=2)<jupyter_output><empty_output><jupyter_text>Finally we form the PDE to be solved.<jupyter_code># PDE for temperature
Ft = ((T-T_n)/dtc)*vt*dx + dot(u,grad(T))*vt*dx + K*dot(grad(T),grad(vt))*dx
a1, L1 = lhs(Ft), rhs(Ft)<jupyter_output><empty_output><jupyter_text>We now use a for loop to march forward in time to solve for temperature.<jupyter_code># Time-stepping
T = Function(V)
#phi = Function(V)
t = 0
for n in range(num_steps):

    # Update current time
    t += dtc
      
    # Compute solution
    # solve(Ft == 0, T)
    solve(a1 == L1, T, Tbc)
             
    # Update previous solution
    T_n.assign(T)<jupyter_output><empty_output><jupyter_text>The final temperature field is then plotted below.<jupyter_code>p = plot(T)
# set colormap
p.set_cmap("viridis")
p.set_clim(0,1450)
# add a title to the plot
plt.title("Temperature")
# add a colorbar
plt.colorbar(p)
plt.ylim(100000,0)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09665508253692441, 'fraction_numerical': 0.015421372719374457, 'mean_word_length': 5.308219178082192, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '41918', 'license': 'Apache-2.0', 'n_tokens_mistral': 1583, 'n_tokens_neox': 1497, 'n_words': 476, 'path': 'submissions/armitage-05-2019/armitage-submission-05-2019.ipynb', 'repo_name': 'softwareunderground/subsurface-journal', 'size': '4614'}