I am trying to adapt a program that simulates a reaction-diffusion soliton made up of substrates X, Y and G. The program includes Navier-Stokes equation to give vortical motion to the substrates.

At present the code gives rotational/vortical motion to the X, Y and G in terms of velocity field u and v, representing velocity vector in x and y directions in this 2D simulation. The z axis being X, Y and G concentration potentials.

In 1D it looks like this:

[![enter image description here](https://i.stack.imgur.com/2iYiq.png)](https://i.stack.imgur.com/2iYiq.png)

I want to change the code so that substrates X, Y and G each have separate velocities (not just the same ones). I.e. X velocity field --> u\_X and v\_X etc. Please can someone help me adapt the code?

Here is what I have done so far:

Original [code](https://github.com/frostburn/tf2-model-g): and [theory](https://github.com/bjdarrer/tf2-model-g/blob/master/docs/overview.pdf) the code is based on.

The code being altered for separate velocity fields is shown [here](https://github.com/bjdarrer/tf2-model-g-velocity).

The code is run with e.g.: `python3 render_video.py ~/tf2-model-g/nucleation_and_motion_in_fluid_2D.mp4 --params params/nucleation_and_motion_in_fluid_2D.yaml`

Firstly, I changed [fluid\_model\_g\_\_\_\_\_\_\_\_1st\_attempt.py#L69](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/fluid_model_g________1st_attempt.py#L69) with:

```
    elif self.dims == 2:
        #self.u = tf.constant(u[0], 'float64') # original 2 lines
        #self.v = tf.constant(u[1], 'float64')

        #BJD adaptions - 6 lines below 9.6.2021 
        self.u_X = tf.constant(u[0], 'float64')
        self.v_X = tf.constant(u[1], 'float64')
        self.u_Y = tf.constant(u[2], 'float64')
        self.v_Y = tf.constant(u[3], 'float64')
        self.u_G = tf.constant(u[4], 'float64')
        self.v_G = tf.constant(u[5], 'float64')

```

and `def diffusion_advection_integrator` @ [fluid\_model\_g\_\_\_\_\_\_\_\_1st\_attempt.py#L159](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/fluid_model_g________1st_attempt.py#L159) :

```
        def diffusion_advection_integrator(G, X, Y, u_X, v_X, u_Y, v_Y, u_G, v_G, divergence1,divergence2, divergence3): # BJD original 10.6.2021

            f_G = self.fft(tf.cast(G, 'complex128'))
            f_X = self.fft(tf.cast(X, 'complex128'))
            f_Y = self.fft(tf.cast(Y, 'complex128'))

            f_G *= decay_G
            f_X *= decay_X
            f_Y *= decay_Y

            G = tf.cast(self.ifft(f_G), 'float64')
            X = tf.cast(self.ifft(f_X), 'float64')
            Y = tf.cast(self.ifft(f_Y), 'float64')

            G_dx = tf.cast(self.ifft(f_G * self.kernel_dx), 'float64')
            G_dy = tf.cast(self.ifft(f_G * self.kernel_dy), 'float64')
            X_dx = tf.cast(self.ifft(f_X * self.kernel_dx), 'float64')
            X_dy = tf.cast(self.ifft(f_X * self.kernel_dy), 'float64')
            Y_dx = tf.cast(self.ifft(f_Y * self.kernel_dx), 'float64')
            Y_dy = tf.cast(self.ifft(f_Y * self.kernel_dy), 'float64')

            #G -= (u*G_dx + v*G_dy + G*divergence) * self.dt # BJD original
            G -= (u_G*G_dx + v_G*G_dy + G*divergence1) * self.dt # BJD change 10.6.2021
            #X -= (u*X_dx + v*X_dy + X*divergence) * self.dt # BJD original
            #Y -= (u*Y_dx + v*Y_dy + Y*divergence) * self.dt # BJD original
            X -= (u_X*X_dx + v_X*X_dy + X*divergence2) * self.dt # BJD change 10.6.2021
            Y -= (u_Y*Y_dx + v_Y*Y_dy + Y*divergence3) * self.dt # BJD change 10.6.2021

            return G, X, Y

```

original code is shown here: [fluid\_model\_g.py#L57](https://github.com/frostburn/tf2-model-g/blob/master/fluid_model_g.py#L57) and here: [fluid\_model\_g.py#L119](https://github.com/frostburn/tf2-model-g/blob/master/fluid_model_g.py#L119)

3 separate densities were included in the new code: @ [fluid\_model\_g\_\_\_\_\_\_\_\_1st\_attempt.py#L297](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/fluid_model_g________1st_attempt.py#L297)

```
    # rho = tf.math.log(self.params['base-density'] + density_of_reactants) # BJD original here
    rho1 = tf.math.log(self.params['base-density1'] + density_of_reactants) # BJD added 10.6.2021
    rho2 = tf.math.log(self.params['base-density2'] + density_of_reactants) # BJD added 10.6.2021
    rho3 = tf.math.log(self.params['base-density3'] + density_of_reactants) # BJD added 10.6.2021

```

and `base-density1`, `base-density2` and `base-density3` were added to the [`yaml`](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/params/nucleation_and_motion_in_fluid_2D___1st_attempt.yaml#L25) file:

```
#base-density: 35.0 ~ BJD original here
base-density1: 35.0 # BJD added 10.6.2021
base-density2: 36.0 # BJD added 10.6.2021
base-density3: 37.0 # BJD added 10.6.2021

```

At the bottom of [`fluid_model_g.py`](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/fluid_model_g________1st_attempt.py#L301), these lines were changed:

```
    if self.dims == 2:
        #u, v = self.u, self.v  # Store unintegrated flow so that we're on the same timestep --- BJD original line 9.6.2021
        u_X, v_X = self.u_X, self.v_X  # Store unintegrated flow so that we're on the same timestep --- BJD added 9.6.2021
        u_Y, v_Y = self.u_Y, self.v_Y  # Store unintegrated flow so that we're on the same timestep --- BJD added 9.6.2021
        u_G, v_G = self.u_G, self.v_G  # Store unintegrated flow so that we're on the same timestep --- BJD added 9.6.2021

        #self.u, self.v, divergence = self.flow_integrator(rho, self.u, self.v) -BJD original line 9.6.2021
        self.u_X, self.v_X, divergence1 = self.flow_integrator(rho1, self.u_X, self.v_X) #--- BJD added 9.6.2021
        self.u_Y, self.v_Y, divergence2 = self.flow_integrator(rho2, self.u_Y, self.v_Y) #--- BJD added 9.6.2021
        self.u_G, self.v_G, divergence3 = self.flow_integrator(rho3, self.u_G, self.v_G) #--- BJD added 9.6.2021

        #self.G, self.X, self.Y = self.diffusion_advection_integrator(self.G, self.X, self.Y, u, v, divergence) # BJD original line 10.6.2021
        self.G, self.X, self.Y = self.diffusion_advection_integrator(self.G, self.X, self.Y, u_X, v_X, u_Y, v_Y, u_G, v_G, divergence1, divergence2, divergence3) #--- BJD added 10.6.2021

```

Looking at line: [fluid\_model\_g\_\_\_\_\_\_\_\_1st\_attempt.py#L47](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/fluid_model_g________1st_attempt.py#L47)

```
    #super().__init__(dx, dt, concentration_G.shape) # BJD original 13.6.2021
    super().__init__(dx, dt, concentration_G.shape, concentration_X.shape, concentration_Y.shape) # BJD change 13.6.2021

```

I started to make changes in file `pde_solver.py` @ [pde\_solver\_\_\_\_\_\_\_\_1st\_attempt.py#L10](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/pde_solver________1st_attempt.py#L10) and [pde\_solver\_\_\_\_\_\_\_\_1st\_attempt.py#L16](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/pde_solver________1st_attempt.py#L16) and [pde\_solver\_\_\_\_\_\_\_\_1st\_attempt.py#L55](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/pde_solver________1st_attempt.py#L55); [pde\_solver\_\_\_\_\_\_\_\_1st\_attempt.py#L81](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/pde_solver________1st_attempt.py#L81) and [pde\_solver\_\_\_\_\_\_\_\_1st\_attempt.py#L99](https://github.com/bjdarrer/tf2-model-g-velocity/blob/master/pde_solver________1st_attempt.py#L99)

```
def __init__(self, dx, dt, shape_G, shape_X, shape_Y): # shape_G, shape_X, shape_Y ::: put in by David J. 13.6.2021
    self.dx = dx
    self.dt = dt
    self.t = 0

    omega = []
    for s in shape_G:  # shape_G, shape_X, shape_Y, i.e. 3 for loops made; original 1 for loop ::: put in by David J. 13.6.2021
        wave_numbers = np.arange(s)
        wave_numbers -= s * (2*wave_numbers > s)  # Deal with TensorFlow's uncentered FFT
        expected_span = 2*np.pi
        actual_span = s*dx
        omega.append(wave_numbers * expected_span / actual_span)
    for s in shape_X:
        wave_numbers = np.arange(s)
        wave_numbers -= s * (2*wave_numbers > s)  # Deal with TensorFlow's uncentered FFT
        expected_span = 2*np.pi
        actual_span = s*dx
        omega.append(wave_numbers * expected_span / actual_span)
    for s in shape_Y:
        wave_numbers = np.arange(s)
        wave_numbers -= s * (2*wave_numbers > s)  # Deal with TensorFlow's uncentered FFT
        expected_span = 2*np.pi
        actual_span = s*dx
        omega.append(wave_numbers * expected_span / actual_span)
    self.omega = np.meshgrid(*omega, indexing='ij')
    self.dims = len(shape_G) + len(shape_X) + len(shape_Y)
    # The naming is a bit off. These are not actual 'kernels'.
    # They are discrete fourier transforms of the periodic versions of the kernels
    if self.dims == 1:
        self.fft = tf.signal.fft
        self.ifft = tf.signal.ifft
        self.omega_x = self.omega[0]
    elif self.dims == 2:
        self.fft = tf.signal.fft2d
        self.ifft = tf.signal.ifft2d

        self.omega_x = self.omega[0]
        self.omega_y = self.omega[1]
    elif self.dims == 3:
        self.fft = tf.signal.fft3d
        self.ifft = tf.signal.ifft3d

        self.omega_x = self.omega[0]
        self.omega_y = self.omega[1]
        self.omega_z = self.omega[2]
    elif self.dims == 6: # With 3 fluids 2 dimensions # :::: put in by David J. 13.6.2021
        self.fft = tf.signal.fft2d
        self.ifft = tf.signal.ifft2d

        self.omega_x = self.omega[0]
        self.omega_y = self.omega[1]
        self.omega_X_x = self.omega[2]
        self.omega_X_y = self.omega[3]
        self.omega_Y_x = self.omega[4]
        self.omega_Y_y = self.omega[5]
    else:
        raise ValueError('{} dimensions not supported'.format(self.dims))

class PDESolverDx(PDESolver):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        if self.dims == 1:
            self.kernel_dx = tf.constant(1j * self.omega_x, 'complex128')
        elif self.dims == 2:
            self.kernel_dx = tf.constant(1j * self.omega_x, 'complex128')
            self.kernel_dy = tf.constant(1j * self.omega_y, 'complex128')
        elif self.dims == 3:
            self.kernel_dx = tf.constant(1j * self.omega_x, 'complex128')
            self.kernel_dy = tf.constant(1j * self.omega_y, 'complex128')
            self.kernel_dz = tf.constant(1j * self.omega_z, 'complex128')
        elif self.dims == 6:   # :::: put in by David J. 13.6.2021
            self.kernel_dx = tf.constant(1j * self.omega_x, 'complex128')
            self.kernel_dy = tf.constant(1j * self.omega_y, 'complex128')
            self.kernel_X_dx = tf.constant(1j * self.omega_X_x, 'complex128')
            self.kernel_X_dy = tf.constant(1j * self.omega_X_y, 'complex128')
            self.kernel_Y_dx = tf.constant(1j * self.omega_Y_x, 'complex128')
            self.kernel_Y_dy = tf.constant(1j * self.omega_Y_y, 'complex128')

class PDESolverDx2(PDESolverDx):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        if self.dims == 1:
            self.kernel_laplacian = tf.constant(-self.omega_x**2, 'complex128')
        elif self.dims == 2:
            self.kernel_laplacian = tf.constant(-(self.omega_x**2 + self.omega_y**2), 'complex128')
        elif self.dims == 3:
            self.kernel_laplacian = tf.constant(-(self.omega_x**2 + self.omega_y**2 + self.omega_z**2), 'complex128')
        elif self.dims == 6:  # :::: put in by David J. 13.6.2021
            self.kernel_laplacian = tf.constant(-(self.omega_x**2 + self.omega_y**2), 'complex128')
            self.kernel_laplacian_X = tf.constant(-(self.omega_X_x**2 + self.omega_X_y**2), 'complex128')
            self.kernel_laplacian_Y = tf.constant(-(self.omega_Y_x**2 + self.omega_Y_y**2), 'complex128')

```

The original is shown here: [pde\_solver.py#L10](https://github.com/frostburn/tf2-model-g/blob/master/pde_solver.py#L10)

This is as far as I got! Can anyone help me solve this program change for including separate velocity fields for X, Y and G substrates?