Multi Steady State Pressure difference driven

python

from opentidalfarm import *
set_log_level(INFO)

# Some domain information extracted from the geo file
basin_x = 2000.
basin_y = 500.

headland_x = 500.
headland_y = 200.

site_x = 300.
site_y = basin_y-headland_y
site_x_start = basin_x/2-site_x/2
site_x_end = basin_x/2+site_x/2

site_y_start = headland_y
site_y_end = headland_y + site_y


def load_turbine_pos_from_csv(csv_file):
    print "Reading turbine positions from %s." % csv_file
    f = open(csv_file, "r")
    f.readline()  # Ingore header file
    turbine_positions = []
    turbine_friction = []
    for l in f:
        pos_x, pos_y, friction = [float(l.split(',')[i]) for i in (0, 1, -1)]
        turbine_positions.append((pos_x, pos_y))
        turbine_friction.append(friction)
    return turbine_positions, turbine_friction



config = UnsteadyConfiguration("mesh.xml", inflow_direction = [1, 0])
config.set_site_dimensions(site_x_start, site_x_end, site_y_start, site_y_end)

# Change the parameters such that in fact two steady state problems are solved consecutively
config.params['initial_condition'] = ConstantFlowInitialCondition(config, val=[1, 1, 1])
config.params['theta'] = 1
config.params['start_time'] = 0 
config.params['dt'] = 1 
config.params['finish_time'] = 2 
config.params['include_time_term'] = False
config.params['functional_quadrature_degree'] = 0

config.params["save_checkpoints"] = True
config.params['controls'] = ['turbine_pos', 'turbine_friction']

config.params['depth'] = 20
config.params['diffusion_coef'] = Constant(5)
config.params["print_individual_turbine_power"] = True
config.params['cache_forward_state'] = True

# Work out the expected delta eta for a free-stream of 2.5 m/s (without turbines) 
# by assuming balance between the pressure and friction terms
u_free_stream = 3.25
print0("Target free-stream velocity (without turbines): %f " % u_free_stream)
delta_eta = config.params["friction"](())/config.params["depth"]/config.params["g"]
if config.params["quadratic_friction"]: 
	delta_eta *= u_free_stream**2
else:
	delta_eta *= u_free_stream
delta_eta *= basin_x
print0("Derived head-loss difference to achieve target free-stream: %f." % delta_eta)

# Set Boundary conditions
bc = DirichletBCSet(config)
expl = Expression("delta_eta/2*cos(pi*(t-1))", delta_eta=delta_eta, t=0)
expr = Expression("-delta_eta/2*cos(pi*(t-1))", delta_eta=delta_eta, t=0)
bc.add_analytic_eta(1, expl)
bc.add_analytic_eta(2, expr)
bc.bcs.append(DirichletBC(bc.function_space.sub(0).sub(1), Constant(0), config.domain.boundaries, 1))
bc.bcs.append(DirichletBC(bc.function_space.sub(0).sub(1), Constant(0), config.domain.boundaries, 2))
bc.add_noslip_u(3)
config.params['strong_bc'] = bc

# Place some turbines 
# Temporarly change the site dimensions to put the turbines closer to towards the central symmetry line
#config.set_site_dimensions(basin_x/2-50, basin_x/2+50, site_y_start, site_y_end)
#deploy_turbines(config, nx=2, ny=6)
# Reset site dimensions
#config.set_site_dimensions(site_x_start, site_x_end, site_y_start, site_y_end)

turbine_positions, turbine_friction = load_turbine_pos_from_csv("turbine_info.csv")

config.params["turbine_pos"] = turbine_positions
config.params["turbine_friction"] = turbine_friction
config.info()

rf = ReducedFunctional(config)
rf.load_checkpoint() 

lb_f, ub_f = friction_constraints(config, lb=0., ub=100.)

# Redefine the site dimension to stretch over the whole y direction
config.set_site_dimensions(site_x_start, site_x_end, 0, basin_y)
lb, ub = position_constraints(config) 

bounds = [lb_f + lb, ub_f + ub]

class TurbineSite(SubDomain):
    def inside(self, x, on_boundary):
        return True if site_x_start <= x[0] <= site_x_end else False

domains = MeshFunction('size_t', config.domain.mesh, 2)
domains.set_all(0)
turbine_site = TurbineSite()
turbine_site.mark(domains, 1)
config.site_dx = Measure("dx")[domains]

feasible_area = get_distance_function(config, domains)
ineq_domain = get_domain_constraints(config, feasible_area, attraction_center=(site_x_start+site_x/2, site_y_start+site_y/2))
ineq_min_dist = get_minimum_distance_constraint_func(config, min_distance=20)
ineq = merge_constraints(ineq_min_dist, ineq_domain)

maximize(rf, bounds=bounds, constraints=ineq, method="SLSQP", options={"maxiter": 400})

The mesh.geo for this example is:

"""
basin_x = 2000;
basin_y = 500;

headland_x = 500;
headland_y = 200;

site_x = 300;
site_x_start = basin_x/2-site_x/2;
site_x_end = basin_x/2+site_x/2;

element_size = 3;
element_size_coarse = 30;

Point(1) = {0, basin_y, 0, element_size_coarse};
Point(2) = {basin_x, basin_y, 0, element_size_coarse};

// Generate nodes for the headland
res = 100;
For k In {0:res:1}
x = basin_x/res_k;
b = 100;
y = headland_y_Exp(-0.5*((x-basin_x/2)/b)^2);
Point(10+k) = {x, y, 0, element_size_coarse};
EndFor

// Generate lines for the headland

BSpline(100) = { 10 : res+10 };

Line(101) = {10, 1};
Line(102) = {1, 2};
Line(103) = {2, res+10};
Line Loop(104) = {100, -103, -102, -101};

// Set mesh element size

// Define the attractor
Field[1] = Attractor;
Field[1].EdgesList = {101, 103};

Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = element_size_coarse;
Field[2].LcMax = element_size;
Field[2].DistMin = basin_x/2-site_x/2 - 100;
Field[2].DistMax = basin_x/2-site_x/2;

Background Field = 2;

// Generate site nodes
Plane Surface(111) = {104};
Physical Surface(112) = {111, 109};
Physical Line(1) = {101};
Physical Line(2) = {103};
Physical Line(3) = {100, 102};
"""