Commit 9550322b authored by sfritschi's avatar sfritschi
Browse files

Also extended periodic solver

parent e938c604
......@@ -4,7 +4,7 @@ include ${PETSC_DIR}/lib/petsc/conf/rules
# Linking flags needed in setup.py
export PETSC_LIB
.PHONY: all, run, cleanup
.PHONY: all, run, test, cleanup
all: setup.py
python3 $< build_ext -i
......
......@@ -8,8 +8,14 @@ from typing import List, Set, Dict, Tuple # for type hints in argument lists
from .netgen import *
from enum import Enum
class Solver(Enum):
AMG = 0,
CG = 1,
ILU = 2,
PETSC = 3
def solve_flow_periodic(network: Network, mu: float, c: int, \
def solve_flow_periodic(network: Network, solver: Solver, mu: float, c: int, \
P2L: float, tol: float = 1e-20, maxiter: int = 1e3, resinfo: bool = False) \
-> Tuple[Dict[Pore,float],Dict[Throat,float]]:
"""Solve flow problem for periodic network.
......@@ -23,44 +29,55 @@ def solve_flow_periodic(network: Network, mu: float, c: int, \
The first return value is a pore/pressure dictionary and the second a
throat/flux dictionary. Fluxes going out of pore1 of a throat are positive.
"""
L = [ub - lb for ub, lb in zip(network.ub,network.lb)]
# initialize
iterator = iter(network.pores)
seed = next(iterator)
tree = __grow({seed})
# find sufficiently large tree of connected pores
while len(tree) < len(network.pores)/2:
# find seed outside of tree being too small
for seed in iterator:
if (seed not in tree): break
else: raise AssertionError('Fragmented network!')
# grow tree based on new seed
tree = __grow({seed})
# assemble sparse pressure matrix
pore_dict = {pore:k for k, pore in enumerate(tree)}
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
root = 0
A = {}; b = {} # pressure matrix and right hand side vector
for p_i, i in pore_dict.items(): # loop over pores i
if i == 0: # boundary condition
A[(0,0)] = 1.0; b[0] = 0.0
continue
for t_ij in p_i.throats:
# pressure system matrix A
T,p_j_pos = __add_coefficient_to(A, t_ij, pore_dict, p_i, i, mu, L)
# right hand side b
__add_to(b, i, T * (p_i.pos[c-1] - p_j_pos[c-1]) * P2L)
p = __amg_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
pore_dict = {}
Q = {}
if rank == root:
L = [ub - lb for ub, lb in zip(network.ub,network.lb)]
# initialize
iterator = iter(network.pores)
seed = next(iterator)
tree = __grow({seed})
# find sufficiently large tree of connected pores
while len(tree) < len(network.pores)/2:
# find seed outside of tree being too small
for seed in iterator:
if (seed not in tree): break
else: raise AssertionError('Fragmented network!')
# grow tree based on new seed
tree = __grow({seed})
# assemble sparse pressure matrix
pore_dict = {pore:k for k, pore in enumerate(tree)}
for p_i, i in pore_dict.items(): # loop over pores i
if i == 0: # boundary condition
A[(0,0)] = 1.0; b[0] = 0.0
continue
for t_ij in p_i.throats:
# pressure system matrix A
T,p_j_pos = __add_coefficient_to(A, t_ij, pore_dict, p_i, i, mu, L)
# right hand side b
__add_to(b, i, T * (p_i.pos[c-1] - p_j_pos[c-1]) * P2L)
if solver == Solver.PETSC:
p = __petsc_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
elif solver == Solver.AMG:
p = __amg_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
elif solver == Solver.CG:
p = __cg_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
else:
p = __ilu_gmres_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
# calculate resulting fluxes
pm = lambda x: P2L*L[c-1] - (x[c-1] - network.lb[c-1]) * P2L # press. drop
Q = __calculate_fluxes(network, pore_dict, p, mu, L, pm)
if rank == root:
pm = lambda x: P2L*L[c-1] - (x[c-1] - network.lb[c-1]) * P2L # press. drop
Q = __calculate_fluxes(network, pore_dict, p, mu, L, pm)
return p, Q
from enum import Enum
class Solver(Enum):
AMG = 0,
CG = 1,
ILU = 2,
PETSC = 3
def solve_flow_inout(network: Network, pin: float, pout: float, \
inpores: Set[Pore], outpores: Set[Pore], solver: Solver, \
mu: float, tol: float = 1e-20, maxiter: int = 1e4, resinfo: bool = False) \
......@@ -74,7 +91,6 @@ def solve_flow_inout(network: Network, pin: float, pout: float, \
The first return value is a pore/pressure dictionary and the second a
throat/flux dictionary. Fluxes going out of pore1 of a throat are positive.
"""
from time import perf_counter
from mpi4py import MPI
comm = MPI.COMM_WORLD
......@@ -113,7 +129,7 @@ def solve_flow_inout(network: Network, pin: float, pout: float, \
else:
p = __ilu_gmres_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
except:
print("Error in AMG solver: switching to ILU-GMRES solver")
print("Error in solver: switching to ILU-GMRES solver")
p = __ilu_gmres_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
# calculate resulting fluxes (only on root)
if rank == root:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment