# Lecture 11:  Scalar and Vector Fields, Accumulation 

## Reading and Reference
* Essential Mathematical Methods for Physicists, H. Weber and G. Arfken, Academic Press, 2003

## What to Learn?

* Define a scalar and vector field
* Define and use the gradient and divergence operators
* Develop a geometric intuition for the divergence and gradient
* Learn or review the definition of the chemical potential and ideal solution models
* The assumptions required to write the diffusion equation in a simplified form

## What to do?

* Compute the gradient of a scalar field and the divergence of a vector field
* Develop visualization tools
* Derive the accumulation formula from the flux into/out of a control volume
* Derive Fick's first law starting from the chemical potential

### Scalar Fields
----

Engineers and scientists are typically interested in how certain quantities vary within a domain.  A scalar field is a mapping of values to positions within some coordinate system.  In physical problems it is common to map values to three dimensions or fewer.  For example, scalar fields may be generated from temperature measurements within a chemical reactor, pressure data over the Earth's surface, or x-ray intensity from telescope observations.  Any mapping of a numerical quantity to a spatial coordinate is a scalar field.

Materials engineers interested in phase transformations typically deal with scalar fields corresponding to concentration as a function of position.  This could be indicated using the symbols $c(x,y)$ where the $c$ corresponds to the field of interest and the $(x,y)$ indicating the dimensions of and position within the field.

$$
f(x,y) = \sin(2x) \cdot \sin(y)
$$

Note that this function takes an $(x,y)$ position as input and returns a single scalar value.


```python
%matplotlib notebook

import sympy as sp
from sympy.plotting import plot3d
import numpy as np
import matplotlib.pyplot as plt
x, y = sp.symbols('x y')
```


```python
# Plot our scalar function over a specified range.
plot3d(sp.sin(2*x)*sp.sin(y), (x, -5, 5), (y, -5, 5));
```

We use `plot3d()` to visualize this field so that we use the height of the surface as a proxy for the magnitude of the scalar field.  


```python
# Here is another example using a more complex `sympy` expression:

plot3d(sp.series(sp.sin(2*x)*sp.sin(y), x, x0=0, n=7).removeO(), \
       (x, -5, 5), (y, -5, 5));
```

### Projections Onto Two Dimensions
----

* A field $f(x,y)$ contains three pieces of information:  f, x, y
* Colors and other glyphs can help access the additional information
* Contours and heatmaps are two such methods

An alternative to the three dimensional plot is to project the scalar values onto a two dimensional surface.  Two common options for this are _contour plots_ and _heat maps_.  Rather than using a third dimension to represent the scalar values a contour plot traces single-valued lines through the domain whereas a heat-map uses colors to represent the value of the scalar field.


### Visualizing the Scalar Field by Contours
----

As previously stated, _iso-lines_ are plotted within the domain and the color, annotations and position of the lines quantify the scalar field.  Topographic maps use this representation to help the reader understand the locations and incline of mountains and valleys.  `matplotlib` has a `.contour()` method that will generate contour plots.  We demonstrate this next.


```python
delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z = np.sin(2*X)*np.sin(Y)
```


```python
contours = 10
plt.figure()
CS = plt.contour(X, Y, Z, contours)
plt.clabel(CS, inline=1, fontsize=10)
plt.title(r'A Simple Contour Plot of Your Scalar Field')
plt.show()
```

### Visualizing the Scalar Field by Value
----

Color choice requires consideration of the medium in which you distribute your data as well as the capabilities of the reader.  An [essay](http://blogs.nature.com/onyourwavelength/2018/07/18/heuristics-for-better-figures/) has been written on the topic and provides some guidance for preparation of figures.  `matplotlib` has a module called `cm` that provides some standard color maps.  User specified maps are possible.

We use our same scalar field and the `.meshgrid()` method to create an array $Z$ that holds our field data.  These data are passed to `imshow` to display the data as an image on the plot axes.  The relationship between the field coordinates and the data is lost when plotted this way so we must specify that directly using the `extent` keyword argument.


```python
print(ax)
```


```python
from matplotlib import cm

delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z = np.sin(2*X)*np.sin(Y)

fig, ax = plt.subplots()
cax = ax.imshow(Z, interpolation='nearest', cmap=cm.coolwarm, extent=(-3,3,-2,2))
# Empirical hack to make colorbar not ugly.
fig.colorbar(cax, fraction=0.032, pad=0.04)
ax.set_title('Scalar Field')
plt.show()
```

### DIY:  Using Surface and Contour Plots
----

Plot the electric field around a point charge if the potential is given by:

$$
V(x,y) = \frac{1}{x^2+y^2}
$$

Use three dimensional surface plots and contour plots to help visualize this potential.


```python
# Your code here.
```

### Gradients
----

Scalar fields have associated vector fields.  One such vector field is known as the _gradient_ and is indicated with the symbol $\nabla$ and is called [nabla](https://en.wikipedia.org/wiki/Nabla_symbol).  The gradient is a type of directional derivative and is a vector quantity.  In a physical problem of three dimensions the gradient can be thought of as pointing in the direction of the maximum spatial rate of change.  Each basis vector magnitude is multiplied by the partial derivative of the field with respect to the basis vector's coordinate.  

You can visualize the gradient operator as a vector with components:

$$\overrightarrow{\nabla} = \frac{\partial}{\partial x} \hat{i}  + \frac{\partial}{\partial y} \hat{j} + \frac{\partial}{\partial z} \hat{k}  $$

and then proceed to perform multiplication as per vector algebra using the dot product.  

There are multiple ways to interact with Python and get at the gradient of a function.  In the first instance we can use the coordinate system capabilities of `sympy` so that we can access the built in method `.gradient()`.  We start by defining a coordinate system and then calling gradient on our scalar function.  Scalars and vectors are objects of the coordinate system.  [See this page](http://docs.sympy.org/latest/modules/vector/intro.html) for more information on the vector module.

For the purposes of numeric computing, `NumPy` has functions and operators that work as you would expect for vector operations.  What is presented below is for convenience, only.


```python
import sympy as sp
import sympy.vector as spv

C = spv.CoordSys3D('C')
spv.gradient(sp.sin(2*C.x)*sp.sin(C.y))
# The gradient function should return something that looks like u+v
# where u is a vector in the x direction and v is a vector in the
# y direction.
```

### Quiver Plots

One way to visualize a vector field in a physical problem is to use small oriented arrow to indicate direction and length to indicate magnitude.  These are known as quiver plots and they can be used in conjunction with heat maps to help visualize vector fields.  In the example below the scalar field is sampled with two different point densities.  A higher point density (100 points) for the filled contour plot and a lower point density (30 points) for the quivers.


```python
x0, x1 = (-3,3)
y0, y1 = (-2,2)

# Read the docstrings to understand why the numbers are given
# as complex quantities.  use: ?np.mgrid
Y, X =  np.mgrid[y0:y1:100j, x0:x1:100j]
Y1, X1 =  np.mgrid[y0:y1:30j, x0:x1:30j]

Z = np.sin(2*X)*np.sin(Y)

# u and v here are the results of applying the gradient operation
# to our scalar field.  Probably wise to check this in a seperate
# code block.
u = (2*np.sin(Y1)*np.cos(2*X1))
v = (np.sin(2*X1)*np.cos(Y1))
```


```python
fig, ax = plt.subplots()
filled_contour = plt.contourf(X,Y,Z,5)
fig.colorbar(filled_contour, ax=ax)
plt.quiver(X1,Y1,u,v, color='white')
plt.show()
```

### Divergence
----

If we interpret the gradient field as a measurement of the flow of some quantity through space, then the divergence of that vector field is a measurement of the accumulation of that quantity.  To compute the divergence we compute the dot product of _nabla_ and the vector field.  The resulting quantity is a scalar quantity:

$$
\overrightarrow{\nabla} \cdot \overrightarrow{F} = \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z}
$$


```python
gradientOfOurFunction = spv.gradient(sp.sin(2*C.x)*sp.sin(C.y))
gradientOfOurFunction
```


```python
spv.divergence(gradientOfOurFunction)
```

### Accumulation and Diffusion
----

One of the key problems in science and engineering is the diffusion of an extensive quantity such as energy, charge, or mass.  The gradient of a potential provides a driving force for this diffusion.  

Using the diffusion of mass as an example:

$$\overrightarrow{J} = -M \overrightarrow{\nabla} \cdot {\mu} $$

and then computing the accumulation based on the vector field:

$$\frac{\partial C(x,t)}{\partial t} = - \overrightarrow{\nabla} \cdot \overrightarrow{J} $$

where the minus sign indicates that accumulation occurs antiparallel to the gradient.


```python
# Define your example scalar field (a concentration like 
# variable C(x,y).
exampleField = sp.sin(sp.pi*C.x)*sp.cos(sp.pi*C.y)
exampleField
```

We then use the built in `sympy` function `gradient` to compute the gradient:


```python
# Compute the gradient.
gradientOfField = -spv.gradient(exampleField)
gradientOfField
```

We then compute the divergence of the gradient.  Note the absence of `C.i` and `C.j` in the answer indicating that these are not components of a vector.  (Compare this to the last slide.)


```python
# Compute the divergence.
accumulation = -spv.divergence(gradientOfField)
accumulation
```

### Plotting the Flux (or Accumulation)
----

In this section we drop the vector descriptions as provided by `sympy` and we keep track of the vector components by hand using the mathematical definitions.  Writing this out we have a scalar field:

$$
F(x,y)
$$

the flux written as minus the gradient of the field (with $M=1$):

$$
\overrightarrow{J} = -\nabla F(x,y) = -\left(\frac{\partial F}{\partial x} \hat{i} + \frac{\partial F}{\partial y}\hat{j} \right)
$$

and the accumulation as minus the divergence of the flux:

$$
A = -\nabla \cdot \overrightarrow{J} = \left(\frac{\partial^2 F}{\partial x^2} + \frac{\partial^2 F}{\partial y^2}\right)
$$


```python
x, y = sp.symbols('x y')

concentrationFunction = sp.sin(sp.pi*x)*sp.cos(sp.pi*y)
fluxX = -sp.diff(concentrationFunction,x)
fluxY = -sp.diff(concentrationFunction,y)
accumulationFunction = sp.diff(concentrationFunction,x,2) + sp.diff(concentrationFunction,y,2)

# We use lambdify to permit the functions to take arguments and vectorize
# the computations.
myConcentrationFunction = sp.lambdify((x,y), concentrationFunction, 'numpy')
myFluxX = sp.lambdify((x,y), fluxX, 'numpy')
myFluxY = sp.lambdify((x,y), fluxY, 'numpy')
myAccumulationFunction = sp.lambdify((x,y), accumulationFunction, 'numpy')
```

We use our quiver plotting capability to plot:

* the scalar value for the concentration (or accumulation) as `z`
* the $\hat{i}$ component of the flux as `u`
* the $\hat{j}$ component of the flux as `v`


```python
import numpy as np

x0, x1 = (-1,1)
y0, y1 = (-1,1)
plotResolution = 200

Y, X =  np.mgrid[y0:y1:200j, x0:x1:200j]
# Quivers are on a seperate grid since they clutter things up.
Y1, X1 =  np.mgrid[y0:y1:20j, x0:x1:20j]

Z = myConcentrationFunction(X,Y)
#Z = myAccumulationFunction(X,Y)
u = myFluxX(X1,Y1)
v = myFluxY(X1,Y1)
```


```python
fig, ax = plt.subplots()

plt.contourf(X,Y,Z,20, cmap=cm.coolwarm)
plt.colorbar()
plt.quiver(X1,Y1,u,v, color='white')

plt.show()
```

In this example we use the contour plot to show locations of high (red) and low (blue) concentrations.  From this concentration field we compute the flux and show that as a quiver plot overlaid on the contour plot.  Examination of this figure confirms our physical understanding that mass flows from high to low concentration areas when the chemical potential is proportional to the concentration of species.  

## DIY:  Using the ideal solution model for the Gibbs free energy, *derive* the accumulation of species in a diffusive problem.

The flux is the divergence of the chemical potential.  The chemical potential is the derivative of the Gibbs free energy with respect to composition.  Review the chain rule before you start this problem.  The ideal solution model in terms of mole fraction for a two component system is:

$$
G(X) = X_AG_A + X_BG_B + RT(X \ln X + (1-X) \ln (1-X))
$$

Compute the divergence of the gradient of this function.  Note that the concentration is a function of spatial coordinates.  In this problem assume that the mole fraction $X$ is only a function of $x$, i.e. $X(x)$.  $G_A$ and $G_B$ are constants that depend on the melting temperature and heat of fusion, $R$ is the gas constant and $T$ is the absolute temperature.

This is easiest if you use SymPy and ignore the coordinate system formalities above.  

## DIY:  Read [this](http://blogs.nature.com/onyourwavelength/2018/07/18/heuristics-for-better-figures/) article and write a 200 word essay on choosing colors for heat maps.

41457525b8ce5478d104d9787c24d131cbdec2cd
ipynb
Jupyter Notebook
Lecture-11-Scalars-Vectors.ipynb
juhimgupta/MTLE-4720
41797715111636067dd4e2b305a782835c05619f
MIT
Lecture-11-Scalars-Vectors.ipynb
juhimgupta/MTLE-4720
41797715111636067dd4e2b305a782835c05619f
MIT
Lecture-11-Scalars-Vectors.ipynb
juhimgupta/MTLE-4720
41797715111636067dd4e2b305a782835c05619f
MIT
Qwen/Qwen-72B
1. YES
2. YES