# Advanced finite difference

## Contents

# Advanced finite difference¶

This notebook assumes that you have completed the *finite difference operations* notebook.

```
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```

## 1 Spatially variable diffusivity¶

For most purposes we want \(\kappa\) to be spatially variable. The heat equation can be expressed as,

The corresponding finite difference approximations in the \(x\) and \(y\) coordinates are,

where \( \kappa_{i+1/2,j} \) can be averaged by,

EXERCISE 1Construct your own matrices using the finite difference approximation for non-constant diffusivity.

## 2. Modelling complex geometries¶

Now for the fun part. In this section we take a 2D geological cross section (from some random part of the world) and model the temperature variation across different rock types. To do this requires assignment of thermal properties, \( \kappa, H \), to each layer in the cross section.

Before you attempt this, here is a checklist of what you should have accomplished:

2D steady-state heat solver for non-constant \( \kappa \).

Neumann and Dirichlet boundary conditions.

Use of sparse matrices (optional, but recommended).

Our model domain is 800 x 260 nodes and contains integers from 1 to 28 that correspond to a unique rock layer. Your task will be to create a conductivity and heat production field, that vary with lithology, and pass them to your solver.

EXERCISE 2Assign thermal properties to the voxel data and solve steady-state diffusion

Create 3 plots: temperature, conductivity, and heat production. (You might want to give them distinct colourmaps to help distinguish them!)

```
voxel = np.load('voxel_data.npz')['data']
voxel.shape
```

```
fig = plt.figure(1, figsize=(20, 5))
ax1 = fig.add_subplot(111)
im1 = ax1.imshow(voxel, origin='lower', cmap='Paired', vmin=1, vmax=28)
fig.colorbar(im1, ax=ax1)
plt.show()
```

What we have is a 2D voxel model that was exported from GoCAD. Each colour represents a unique lithology. Pretty eh?

Now we will assign these thermal properties and solve temperature across them.