problem.py

"""
problem.py
"""
import numpy as np
import time
from scipy.interpolate import interp1d
from scipy.optimize import minimize, root

"""
 Sample1D class
"""
class Sample1D:
	def __init__(self,species=1,dim=(64),elementsize=1):
		self.species=species
		self.dim=dim
		self.elementsize=elementsize
		self.z=np.zeros(species)
		self.c1_bulk=np.zeros(species)
		self.c2_bulk=np.zeros(species)
		# The phenomena happens in 100nm circa and gets relevant in minutes.
		# Thus the scale we'll be using is going to be nm/s
		self.D1=np.zeros(species) # Last product accounts for cm^2/s to nm^2/s conversion
		self.D2=np.zeros(species) # Last product accounts for cm^2/s to nm^2/s conversion
		# Parameters
		self.L=0
		self.W0=0 # J/m
		self.c=np.zeros((dim,species))
		self.eta=np.zeros(dim)

	def prepare(self):
		# INTRODUCE ALL NECESSARY CONSISTENCY CHECKS HERE! THIS COMMAND WILL BE RUN BEFORE ANY SOLUTION ATTEMPT.
		# - vector dimension coherence
		# - warnings for zero values
		# - output of the considered parameters
		# Scaling factors
		self.L_s=self.L*(1/self.elementsize**3) # L is nm**3/(J*s)
		self.W0_s=self.W0*(self.elementsize) # W_0 is J/(nm)	
		# Prepare diffusion matrices
		self.D1mat=(1/self.elementsize**2)*(np.diag(self.D1[0:-1])-(self.z[0:-1]*self.c1_bulk[0:-1]*self.D1[0:-1])[:,np.newaxis]*\
			(self.z[0:-1]*(self.D1[0:-1]-self.D1[-1]))[np.newaxis,:]/np.sum(self.z**2*self.c1_bulk*self.D1))
		self.D2mat=(1/self.elementsize**2)*(np.diag(self.D2[0:-1])-(self.z[0:-1]*self.c2_bulk[0:-1]*self.D2[0:-1])[:,np.newaxis]*\
			(self.z[0:-1]*(self.D2[0:-1]-self.D2[-1]))[np.newaxis,:]/np.sum(self.z**2*self.c2_bulk*self.D2))
	
	def h(self):
	    return self.eta**3*(6*self.eta**2-15*self.eta+10)
	
	def D(self):
		return self.h()[...,np.newaxis,np.newaxis]*self.D1mat+(1-self.h())[...,np.newaxis,np.newaxis]*self.D2mat
	
	# Applying the boundary conditions by setting D=0 at the borders blocks
	# the entire diffusion process!
	# ~ def applyboundary(D):
		# ~ D[0,:,:,...]=0
		# ~ D[:,0,:,...]=0
		# ~ D[:,:,0,...]=0
		# ~ D[-1,:,:,...]=0
		# ~ D[:,-1,:,...]=0
		# ~ D[:,:,-1,...]=0
		# ~ return D
	
	# The free energy function is useful only during the study phase
	# ~ def f_f(c,c1_bulk,c2_bulk,eta):
		# ~ return h_f(eta)*np.sum((c[0:-1]-c1_bulk[0:-1])**2,axis=0) \
		# ~ +(1-h_f(eta))*np.sum((c[0:-1]-c2_bulk[0:-1])**2,axis=0) \
		# ~ +2*(eta**4-2*eta**3+eta**2)

	def dfdeta(self):
		# Expression obtained through the SymPy stack
		# See equations.py for explanation
		# Can probably be simplified a little bit
		# Input vectors must already be dimensionally consistent
		return (self.elementsize**3)*\
		(2*self.eta*(self.eta - 1)*\
		(15*np.sum((self.c[...,0:-1]-self.c1_bulk[0:-1])**2,axis=-1)*self.eta**2 \
		-15*np.sum((self.c[...,0:-1]-self.c1_bulk[0:-1])**2,axis=-1)*self.eta \
		-15*np.sum((self.c[...,0:-1]-self.c2_bulk[0:-1])**2,axis=-1)*self.eta**2 \
		+15*np.sum((self.c[...,0:-1]-self.c2_bulk[0:-1])**2,axis=-1)*self.eta \
		+4*self.eta - 2))
	
	def cahnhilliard(self):
		# D coefficients are applied BEFORE second derivative (APL104 method)
	    # NOTE: axis numbering goes from external to internal!
	    # Status: RC1
	    return np.gradient(np.sum(self.D()*np.gradient(self.c,axis=0)[...,np.newaxis,:],axis=-1),axis=0)

	def allencahn(self):
	    return self.L_c*self.L*(self.W0_c*self.W0*np.gradient(np.gradient(self.eta))-self.dfdeta())
	
	def fouriersolve(self,dt,steps):
		
		# Prepare static computational parameters
		self.prepare()
		
		# Prepare discretized space
		# STATUS: RC
		kx=np.concatenate([np.arange(0,np.pi,2*np.pi/self.dim),np.arange(-np.pi,0,2*np.pi/self.dim)])
		k=np.sqrt(kx**2)
		
		# Annotate starting time
		starttime=time.clock()
		
		# Prepare initial transformation
		c_fft=np.fft.fft(self.c[...,0:-1],axis=0)
		eta_fft=np.fft.fft(self.eta)
		# Semi-implicit solution algorithm
		for t in range(steps-1):
			print("Step {}".format(t), end='\r', flush=True)
			# Phase evolution equation in fourier space
			# STATUS: Beta
			dfdeta_fft=np.fft.fft(self.dfdeta()) # Can probably be expressed in fourier space
			eta_fft=(eta_fft-self.L_s*dt*dfdeta_fft)/(1+dt*self.L_s*self.W0_s*k**2)
			self.eta=np.real(np.fft.ifft(eta_fft)) # New eta
			# Concentration evolution equation in fourier space
			# STATUS: Beta
			den=np.eye(self.species-1)[np.newaxis,...]+dt*self.D()*k[...,np.newaxis,np.newaxis]**2
			c_fft=np.linalg.solve(den,c_fft) # Linear system of equations
			self.c[...,0:-1]=np.real(np.fft.ifft(c_fft,axis=0))
			self.c[...,-1]=1-np.sum(self.c[...,0:-1],axis=-1)
		endtime=time.clock()
		print("\nDone in {} seconds".format(endtime-starttime))

	def fouriersolvesilent(self,dt,steps):
		# Prepare static computational parameters
		self.prepare()
		
		# Prepare discretized space
		# STATUS: RC
		kx=np.concatenate([np.arange(0,np.pi,2*np.pi/self.dim),np.arange(-np.pi,0,2*np.pi/self.dim)])
		k=np.sqrt(kx**2)
				
		# Prepare initial transformation
		c_fft=np.fft.fft(self.c[...,0:-1],axis=0)
		eta_fft=np.fft.fft(self.eta)
		# Semi-implicit solution algorithm
		for t in range(steps-1):
			# Phase evolution equation in fourier space
			# STATUS: Beta
			dfdeta_fft=np.fft.fft(self.dfdeta()) # Can probably be expressed in fourier space
			eta_fft=(eta_fft-self.L_s*dt*dfdeta_fft)/(1+dt*self.L_s*self.W0_s*k**2)
			self.eta=np.real(np.fft.ifft(eta_fft)) # New eta
			# Concentration evolution equation in fourier space
			# STATUS: Beta
			den=np.eye(self.species-1)[np.newaxis,...]+dt*self.D()*k[...,np.newaxis,np.newaxis]**2
			c_fft=np.linalg.solve(den,c_fft) # Linear system of equations
			self.c[...,0:-1]=np.real(np.fft.ifft(c_fft,axis=0))
			self.c[...,-1]=1-np.sum(self.c[...,0:-1],axis=-1)
		
	
"""
 Problem function definition
"""
def problem(guess,target_c):
	sample=Sample1D(6,128,3)
	sample.z=np.array([4,3,3,3,2,-2])
	sample.c1_bulk=np.array([0.3,0.05,0,0,0,0.65])
	sample.c2_bulk=np.array([0,0,0.2,0.14,0.06,0.6])
	sample.L=1
	sample.W0=4 # J/m
	# Set guesses in sample
	sample.D1=guess[0:6]
	sample.D2=guess[6:12]
	# Set initial condition
	sample.c[...]=sample.c2_bulk
	sample.c[int(128*1/4):int(128*3/4)]=sample.c1_bulk
	sample.eta[...]=0
	sample.eta[int(128*1/4):int(128*3/4)]=1
	# Solve sample
	sample.fouriersolve(0.5,3600*2)
	# Evaluate error
	# CONSIDER OTHER METRICS TO EVALUATE THE ERROR IN THE CURVES
	# RMSE all over the curves
	return np.sqrt(np.sum((sample.c-target_c)**2)/(sample.dim*sample.species))
	
"""
 Script
"""
# Prepare a sample target curve
target=Sample1D(6,128,3)
target.z=np.array([4,3,3,3,2,-2])
target.c1_bulk=np.array([0.3,0.05,0,0,0,0.65])
target.c2_bulk=np.array([0,0,0.2,0.14,0.06,0.6])
target.L=1
target.W0=4 # J/m
target.D1=np.array([1.3,11.3,10,1.1,2.07,1e8])*1e-16*1e14 # Last product accounts for cm^2/s to nm^2/s conversion
target.D2=np.array([1.3,11.3,26,2.41,1.16,1e7])*1e-16*1e14 # Last product accounts for cm^2/s to nm^2/s conversion
### 1D Problem
# Symmetric, for fourier solver
target.c[...]=target.c2_bulk
target.c[int(128*1/4):int(128*3/4)]=target.c1_bulk
target.eta[...]=0
target.eta[int(128*1/4):int(128*3/4)]=1
target.fouriersolve(0.5,3600*2)

# Prepare first guess using exact guess with some noise
guess=np.concatenate([target.D1,target.D2])

# Optimize
res=root(problem,guess,target.c,options={'maxfev':1})