Skip to content
Snippets Groups Projects
Commit 47a579b4 authored by David Seus's avatar David Seus
Browse files

add generate_exact_solution_expressions() method to helpersn

parent 070a9dd4
No related branches found
No related tags found
No related merge requests found
......@@ -8,6 +8,9 @@ import sys
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import typing as tp
import sympy as sym
# from termcolor import colored
# from colormap import Colormap
#import fenicstools as ft
......@@ -18,174 +21,105 @@ def print_once(*args,**kwargs):
print(*args,**kwargs)
sys.stdout.flush()
# def fit_IC(mesh,u):
# """
# Fit given data u on mesh to certain functions.
# Useful in the case of a single droplet.
# Returns the fitted data as class for IC.
# """
#
# # get mesh coordiantes
# xx = mesh.coordinates()[:,0]
# yy = mesh.coordinates()[:,1]
# nx = np.size(xx)
# ny = np.size(yy)
#
# # get function values
# rho_data = np.zeros(nx)
# phase_field_data = np.zeros(nx)
# mu_data = np.zeros(nx)
#
# for i in range(0,nx):
# rho_data[i] = u.sub(0)(xx[i],yy[i])
# mu_data[i] = u.sub(4)(xx[i],yy[i])
# phase_field_data[i] = u.sub(2)(xx[i],yy[i])
#
# # define functions to fit
# def rho(x,a,b,cx,cy,dx,dy):
# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
# return val
#
# def phase_field(x,a,b,cx,cy,dx,dy):
# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
# return val
#
# def mu(x,a,b,cx,cy,dx,dy):
# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
# return val
#
# xx1 = np.concatenate((xx[:,np.newaxis],yy[:,np.newaxis]),axis= 1)
# # fit data to function
# popt_rho = curve_fit(rho, xx1.T , rho_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
# [a,b,cx,cy,dx,dy] = popt_rho[0]
# # print(a,b,cx,cy,dx,dy)
#
# popt_phase_field = curve_fit(phase_field, xx1.T , phase_field_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
# [a1,b1,cx1,cy1,dx1,dy1] = popt_phase_field[0]
# # print(a1,b1,cx1,cy1,dx1,dy1)
#
# popt_mu = curve_fit(mu, xx1.T , mu_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
# [a2,b2,cx2,cy2,dx2,dy2] = popt_mu[0]
# # print(a2,b2,cx2,cy2,dx2,dy2)
#
# # plot result
# # res = 500;
# # [X,Y] = np.meshgrid(np.linspace(0.,1.,res),np.linspace(0.,1.,res))
# # Z = np.zeros((res,res))
# # for i in range(0,res):
# # for j in range(0,res):
# # Z[i,j] = phase_field([X[i,j],Y[i,j]],a1,b1,cx1,cy1,d1)
# # plt.imshow(Z,cmap = "viridis", origin ="lower")
# # plt.colorbar()
# # plt.show()
#
# class fitted_IC(Expression):
# def eval(self, values, x):
# # rho
# values[0] = rho(x,a,b,cx,cy,dx,dy)
#
# # velocity
# for i in range(2):
# values[i + 1] = 0.
# # c
# values[2 + 1] = phase_field(x,a1,b1,cx1,cy1,dx1,dy1)
# # tau
# values[2 + 2] = 0.0
# # mu
# values[2 + 3] = mu(x,a2,b2,cx2,cy2,dx2,dy2)
# # sigma
# for i in range(2):
# values[2 + 4 + i] = 0.
#
# def value_shape(self):
# return (4 + 2*2,)
#
# return fitted_IC
#
#
# def tensor_jump(v, n):
# """
# DG jump operator for vector valued functions
# """
# return df.outer(v("+"), n("+")) + df.outer(v("-"), n("-"))
#
# def print_constants(self,file=None):
# """
# Print the names and values of all df.Constant attributes.
# """
# for key in self.__dict__.keys():
# attribute = self.__dict__[key]
# if type(attribute) is type(df.Constant(0.)):
# print_once(attribute.name() + " = " + str(attribute.values()),file=file)
# for key in self.__dict__["_EOS"].__dict__:
# attribute = self.__dict__["_EOS"].__dict__[key]
# if type(attribute) is type(df.Constant(0.)):
# print_once(attribute.name() + " = " + str(attribute.values()),file=file)
#
# def get_interface_thickness(solution,ME,xmin=-.5,xmax=1.5,ymin=0.,ymax=1.,ny=1000):
# """
# Compute interface thickness from simulation result.
# The interface thickness is here defined as the region where 0.1 < c < 0.9.
# Parameter: solution, function space ME
# Meshrelated: xmin,xmax,ymin,ymax
# Number of generated samples: ny
# """
# xmid = .5*(xmin+xmax)
# yy = np.linspace(ymin,ymax,ny)
# xx = xmid*np.ones_like(yy)
# # evaluate phasefield along following points
# zz = np.column_stack((xx,yy))
# p = ft.Probes(zz.flatten(), ME.sub(2).collapse())
# p(solution.sub(2))
# # values collected on root process
# values= p.array(root=0)
# if df.MPI.rank(df.MPI.comm_world) == 0:
# try:
# # get interface_start
# interface_start = np.min(np.where(abs(values-0.1)<5.e-3))
# # get interface_end
# interface_end = np.min(np.where(abs(values-0.9)<5.e-3))
# # compute interface_thickness
# interface_thickness = abs(yy[interface_end]-yy[interface_start])
# return interface_thickness
# except:
# print_once(colored("Phasefield values out of range.","red"))
# return -1
#
#
# def plot_solution(u):
# """ Plot solution for poster or the like..."""
# # Define "CD"-Colormap
# red = [0.12549019607843137, 0.3843137254901961, 0.6431372549019608, 0.9019607843137255]
# blue = [0.6823529411764706, 0.788235294117647, 0.8941176470588236, 1.0]
# green = [0.33725490196078434, 0.5254901960784314, 0.7137254901960784, 0.9019607843137255]
# c = Colormap()
# mycmap = c.cmap( {'red':red[::-1], 'green':green[::-1], 'blue':blue[::-1]})
# matplotlib.cm.register_cmap(name='mycmap', cmap=mycmap)
# # plt.rcParams['image.cmap'] = 'mycmap'
# plt.rcParams['image.cmap'] = 'RdYlBu'
# matplotlib.rcParams.update({'font.size': 17})
# matplotlib.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
# matplotlib.rc('text', usetex=True)
#
# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
# # fixes to remove white contour lines
# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
#
# for c in cnt.collections:
# c.set_edgecolor("face")
#
# plt.xlim(0.45, 0.8)
# plt.ylim(0.2, 0.55)
# plt.axis("off")
# plt.annotate("$+$ phase",color='w',xy=(0.5, .45))
# plt.annotate("$-$ phase",xy=(0.7, .3))
# plt.annotate("$\Omega$",xy=(0.75, .5),fontsize=25)
# plt.annotate("$\gamma$",xy=(0.647, .358),fontsize=25)
# plt.annotate("",xy=(0.62,0.375),xytext=(0.66,0.33),arrowprops=dict(arrowstyle="<->",
# connectionstyle="arc3",facecolor='black',linewidth=1))
# #plt.show()
# plt.tight_layout()
# plt.savefig('diffuse_simulation.pdf',bbox_inches='tight')
# return 0
def generate_exact_solution_expressions(
symbols: tp.Dict[int, int],
isRichards: tp.Dict[int, bool],
symbolic_pressure: tp.Dict[int, tp.Dict[str, str]],
symbolic_capillary_pressure: tp.Dict[int, str],
viscosity: tp.Dict[int, tp.List[float]],#
porosity: tp.Dict[int, tp.Dict[str, float]],#
relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
relative_permeability_prime: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],
saturation_pressure_relationship: tp.Dict[int, tp.Callable[...,None]],#
saturation_pressure_relationship_prime: tp.Dict[int, tp.Callable[...,None]],#
densities: tp.Dict[int, tp.Dict[str, float]] = None,#
gravity_acceleration: float = 9.81,
include_gravity: bool = False,
) -> tp.Dict[str, tp.Dict[int, tp.Dict[str, str]] ]:
""" Generate exact solution, source and initial condition code from symbolic expressions"""
output = {
'source': dict(),
'initial_condition': dict(),
'exact_solution': dict()
}
x = symbols["x"]
y = symbols["y"]
t = symbols["t"]
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
dtS = dict()
div_flux = dict()
S_pc_sym = saturation_pressure_relationship
S_pc_sym_prime = saturation_pressure_relationship_prime
for subdomain, isR in isRichards.items():
dtS.update({subdomain: dict()})
div_flux.update({subdomain: dict()})
output['source'].update({subdomain: dict()})
output['exact_solution'].update({subdomain: dict()})
output['initial_condition'].update({subdomain: dict()})
if isR:
subdomain_has_phases = ["wetting"]
else:
subdomain_has_phases = ["wetting", "nonwetting"]
# conditional for S_pc_prime
pc = symbolic_capillary_pressure[subdomain]
dtpc = sym.diff(pc, t, 1)
dxpc = sym.diff(pc, x, 1)
dypc = sym.diff(pc, y, 1)
S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
for phase in subdomain_has_phases:
# Turn above symbolic expression for exact solution into c code
output['exact_solution'][subdomain].update(
{phase: sym.printing.ccode(symbolic_pressure[subdomain][phase])}
)
# save the c code for initial conditions
output['initial_condition'][subdomain].update(
{phase: sym.printing.ccode(symbolic_pressure[subdomain][phase].subs(t, 0))}
)
if phase == "nonwetting":
dtS[subdomain].update(
{phase: -porosity[subdomain]*dS*dtpc}
)
else:
dtS[subdomain].update(
{phase: porosity[subdomain]*dS*dtpc}
)
pa = symbolic_pressure[subdomain][phase]
dxpa = sym.diff(pa, x, 1)
dxdxpa = sym.diff(pa, x, 2)
dypa = sym.diff(pa, y, 1)
dydypa = sym.diff(pa, y, 2)
mu = viscosity[subdomain][phase]
ka = relative_permeability[subdomain][phase]
dka = relative_permeability_prime[subdomain][phase]
rho = densities[subdomain][phase]
g = gravity_acceleration
if phase == "nonwetting":
# x part of div(flux) for nonwetting
dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
# y part of div(flux) for nonwetting
if include_gravity:
dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(1-S)
else:
dydyflux = -1/mu*dka(1-S)*dS*dypc*dypa + 1/mu*dydypa*ka(1-S)
else:
# x part of div(flux) for wetting
dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
# y part of div(flux) for wetting
if include_gravity:
dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
else:
dydyflux = 1/mu*dka(S)*dS*dypc*(dypa) + 1/mu*dydypa*ka(S)
div_flux[subdomain].update({phase: dxdxflux + dydyflux})
contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
output['source'][subdomain].update(
{phase: sym.printing.ccode(contructed_rhs)}
)
return output
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment