diff --git a/data_assimilation/ARTSS.py b/data_assimilation/ARTSS.py
index 77e44922..28665b95 100644
--- a/data_assimilation/ARTSS.py
+++ b/data_assimilation/ARTSS.py
@@ -1,11 +1,61 @@
#!/usr/bin/env python3
import os
import xml.etree.ElementTree as ET
-from typing import List, Dict, Union, Set
+from typing import List, Dict, Union, Set, TextIO
import numpy as np
from obstacle import Obstacle, PATCHES, STATE
+from utility import log
+
+
+def start_new_instance(output_file: str, directory: str, artss_exe_path: str,
+ artss_exe: str = 'artss_data_assimilation_serial'):
+ cwd = os.getcwd()
+ os.chdir(directory)
+ exe_command = f'mpirun --np 2 {os.path.join(artss_exe_path, artss_exe)} {output_file}'
+ print(os.getcwd(), exe_command)
+ os.system(exe_command)
+ os.chdir(cwd)
+
+
+def change_xml(change: Dict[str, float], input_file: str, output_file: str, artss_root_path: str, artss_data_path: str, file_debug: TextIO):
+ out_file = os.path.join(artss_data_path, output_file)
+ log(f'out_file: {out_file}', file_debug)
+
+ command = ''
+ for key in change:
+ command += f' --{key} {change[key]}'
+ command = f'{os.path.join(artss_root_path, "tools", "change_xml.sh")} -i {input_file} -o {out_file} --loglevel off' + command
+ log(f'{command}', file_debug)
+ os.system(command)
+
+
+def create_start_xml(change: dict, input_file: str, output_file: str, t_revert: float, file: str, t_end: float,
+ directory: str, artss_path: str, file_debug: TextIO, dir_name: str = '') -> [str, str]:
+ f = os.path.join('..', file)
+ f = f.replace("/", "\/")
+
+ command = ''
+ name = ''
+ for key in change:
+ command += f' --{key} {change[key]}'
+ name += f'{key}_'
+ name = name[:-1]
+
+ if not dir_name == '':
+ name = dir_name
+ output_dir = os.path.join(directory, '.vis', name)
+ if not os.path.exists(output_dir):
+ os.makedirs(output_dir)
+
+ out_file = os.path.join(output_dir, output_file)
+ log(f'out_file: {out_file}', file_debug)
+ log(f'out_dir: {output_dir}', file_debug)
+ command = f'{os.path.join(artss_path, "tools", "change_xml.sh")} -i {input_file} -o {out_file} --da {t_revert} "{f}" 7779 --tend {t_end} --loglevel off' + command
+ log(f'{command}', file_debug)
+ os.system(command)
+ return out_file, output_dir
class XML:
@@ -52,7 +102,7 @@ def get_output_time_step(self) -> float:
return 1
def get_temperature_source(self) -> dict:
- if self.temperature_source is None:
+ if not self.temperature_source:
# check if temperature source is even there ? or crash
root = self.xml_tree.getroot()
source_tree = root.find('solver').find('temperature').find('source')
@@ -89,8 +139,7 @@ def __init__(self, domain_param: Dict[str, float], enable_computational_domain:
for o in obstacles:
self.obstacles[o.name] = o
self.domain_boundary_cells: Dict[str, Set[int]] = {}
- self.domain_inner_cells: Dict[str, Set[int]] = {}
- self.obstacle_cells: Dict[str, Dict[str, Set[int]]] = {}
+ self.domain_inner_cells: Set[int] = set()
self.computational_domain(domain_param, enable_computational_domain)
self.calculate_domain()
self.calculate_indices()
@@ -166,12 +215,12 @@ def calculate_obstacle_cells(self, o: Obstacle):
index_y2: int = self.get_matching_obstacle_index(o.geometry['oy2'], 'y')
index_z1: int = self.get_matching_obstacle_index(o.geometry['oz1'], 'z') + 1
index_z2: int = self.get_matching_obstacle_index(o.geometry['oz2'], 'z')
+ o.index = {'x1': index_x1, 'x2': index_x2, 'y1': index_y1, 'y2': index_y2, 'z1': index_z1, 'z2': index_z2}
for i in range(index_x1 + 1, index_x2):
for j in range(index_y1 + 1, index_y2):
for k in range(index_z1 + 1, index_z2):
inner_cells.append(self.calculate_index(i, j, k))
- self.obstacle_cells[o.name]: Dict[str, Set[int]] = {}
- self.obstacle_cells[o.name]['inner cells'] = set(inner_cells)
+ self.obstacles[o.name].cells['inner cells'] = set(inner_cells)
# indices of obstacle boundary cells
front = []
@@ -180,8 +229,8 @@ def calculate_obstacle_cells(self, o: Obstacle):
for j in range(index_y1, index_y2 + 1):
front.append(self.calculate_index(i, j, index_z1))
back.append(self.calculate_index(i, j, index_z2))
- self.obstacle_cells[o.name]['front'] = set(front)
- self.obstacle_cells[o.name]['back'] = set(back)
+ self.obstacles[o.name].cells['front'] = set(front)
+ self.obstacles[o.name].cells['back'] = set(back)
bottom = []
top = []
@@ -189,8 +238,8 @@ def calculate_obstacle_cells(self, o: Obstacle):
for k in range(index_z1, index_z2 + 1):
bottom.append(self.calculate_index(i, index_y1, k))
top.append(self.calculate_index(i, index_y2, k))
- self.obstacle_cells[o.name]['bottom'] = set(bottom)
- self.obstacle_cells[o.name]['top'] = set(top)
+ self.obstacles[o.name].cells['bottom'] = set(bottom)
+ self.obstacles[o.name].cells['top'] = set(top)
left = []
right = []
@@ -198,8 +247,8 @@ def calculate_obstacle_cells(self, o: Obstacle):
for k in range(index_z1, index_z2 + 1):
left.append(self.calculate_index(index_x1, j, k))
right.append(self.calculate_index(index_x2, j, k))
- self.obstacle_cells[o.name]['left'] = set(left)
- self.obstacle_cells[o.name]['right'] = set(right)
+ self.obstacles[o.name].cells['left'] = set(left)
+ self.obstacles[o.name].cells['right'] = set(right)
def calculate_domain_boundaries(self):
front = []
@@ -236,9 +285,9 @@ def calculate_domain_inner(self):
for k in range(self.domain_param['start z index CD'], self.domain_param['end z index CD']):
inner.append(self.calculate_index(i, j, k))
inner = set(inner)
- for oc in self.obstacle_cells:
- for type_of_cell in self.obstacle_cells[oc]:
- inner = inner - self.obstacle_cells[oc][type_of_cell]
+ for obstacle in self.obstacles.values():
+ for type_of_cell in obstacle.cells:
+ inner = inner - obstacle.cells[type_of_cell]
self.domain_inner_cells = inner
def calculate_index(self, i, j, k):
@@ -250,14 +299,14 @@ def print_info(self):
f"Domain size inner cells: {self.domain_param['nx']} {self.domain_param['ny']} {self.domain_param['nz']}\n"
f"step size (x|y|z): ({self.domain_param['dx']}|{self.domain_param['dy']}|{self.domain_param['dz']})")
for o in self.obstacles.values():
- print(f"-- Obstacle {o.name}\n"
+ print(f"-- Obstacle {o.name} ({o.state})\n"
f" size of slices (Front|Back Bottom|Top Left|Right): "
- f"{len(self.obstacle_cells[o.name]['front'])}|"
- f"{len(self.obstacle_cells[o.name]['back'])} "
- f"{len(self.obstacle_cells[o.name]['bottom'])}|"
- f"{len(self.obstacle_cells[o.name]['top'])} "
- f"{len(self.obstacle_cells[o.name]['left'])}|"
- f"{len(self.obstacle_cells[o.name]['right'])}\n"
+ f"{len(self.obstacles[o.name].cells['front'])}|"
+ f"{len(self.obstacles[o.name].cells['back'])} "
+ f"{len(self.obstacles[o.name].cells['bottom'])}|"
+ f"{len(self.obstacles[o.name].cells['top'])} "
+ f"{len(self.obstacles[o.name].cells['left'])}|"
+ f"{len(self.obstacles[o.name].cells['right'])}\n"
f" coords (x y z): "
f"({o.geometry['ox1']}|{o.geometry['ox2']}) "
f"({o.geometry['oy1']}|{o.geometry['oy2']}) "
@@ -279,23 +328,27 @@ def print_debug(self):
print(f"index X PD: ({self.domain_param['start x index PD']}|{self.domain_param['end x index PD']}) "
f"index Y PD: ({self.domain_param['start y index PD']}|{self.domain_param['end y index PD']}) "
f"index Z PD: ({self.domain_param['start z index PD']}|{self.domain_param['end z index PD']})")
+ print(f"inner cells index: ({min(self.domain_inner_cells)}|{max(self.domain_inner_cells)})")
+ for patch in PATCHES:
+ print(f"boundary cells index {patch}: ({min(self.domain_boundary_cells[patch])}|{max(self.domain_boundary_cells[patch])})")
+ print()
for o in self.obstacles.values():
print(f"-- Obstacle {o.name}\n"
f" size of slices (Front|Back Bottom|Top Left|Right): "
- f"{len(self.obstacle_cells[o.name]['front'])}|"
- f"{len(self.obstacle_cells[o.name]['back'])} "
- f"{len(self.obstacle_cells[o.name]['bottom'])}|"
- f"{len(self.obstacle_cells[o.name]['top'])} "
- f"{len(self.obstacle_cells[o.name]['left'])}|"
- f"{len(self.obstacle_cells[o.name]['right'])}\n"
+ f"{len(self.obstacles[o.name].cells['front'])}|"
+ f"{len(self.obstacles[o.name].cells['back'])} "
+ f"{len(self.obstacles[o.name].cells['bottom'])}|"
+ f"{len(self.obstacles[o.name].cells['top'])} "
+ f"{len(self.obstacles[o.name].cells['left'])}|"
+ f"{len(self.obstacles[o.name].cells['right'])}\n"
f" coords (x y z): "
f"({o.geometry['ox1']}|{o.geometry['ox2']}) "
f"({o.geometry['oy1']}|{o.geometry['oy2']}) "
f"({o.geometry['oz1']}|{o.geometry['oz2']})")
str_obstacle_cells = ''
for patch in PATCHES:
- start = min(self.obstacle_cells[o.name][patch])
- end = max(self.obstacle_cells[o.name][patch])
+ start = min(self.obstacles[o.name].cells[patch])
+ end = max(self.obstacles[o.name].cells[patch])
str_obstacle_cells += f" ({patch}: {start}|{end})\n"
k = self.get_k(start)
j = self.get_j(start, k=k)
@@ -305,7 +358,7 @@ def print_debug(self):
j = self.get_j(end, k=k)
i = self.get_i(end, k=k, j=j)
str_obstacle_cells += f" {patch} end: {i}|{j}|{k}\n"
- print(str_obstacle_cells)
+ print(str_obstacle_cells[:-2])
def get_ijk_from_xyz(self, coord_x, coord_y, coord_z):
return self.get_i_from_x(coord_x), self.get_j_from_y(coord_y), self.get_k_from_z(coord_z)
@@ -336,11 +389,11 @@ def get_type(self, index):
for p in PATCHES.keys():
if index in self.domain_boundary_cells[p]:
matches.append(f'domain boundary cell ({p})')
- for oc in self.obstacle_cells:
+ for oc in self.obstacles:
for p in PATCHES.keys():
- if index in self.obstacle_cells[oc][p]:
+ if index in self.obstacles[oc].cells[p]:
matches.append(f'obstacle boundary cell ({oc}/{p})')
- if index in self.obstacle_cells['inner cells']:
+ if index in self.obstacles[oc].cells['inner cells']:
matches.append(f'obstacle inner cell ({oc})')
if len(matches) == 0:
return "cell type unknown"
@@ -356,17 +409,22 @@ def change_obstacle(self, obstacle: Obstacle):
self.obstacles[obstacle.name] = obstacle
self.calculate_obstacle_cells(obstacle)
- def remove_obstacle(self, name: str):
- self.obstacles[name].state = 'deleted'
+ def remove_obstacle(self, name: str) -> Obstacle:
+ obstacle = self.obstacles.pop(name)
+ obstacle.state = 'deleted'
+ return obstacle
def set_value_of_obstacle_cells(self, value: float, field: np.ndarray, obstacle_name: str):
cells: Set[float] = set()
- for val in self.obstacle_cells[obstacle_name].values():
+ for val in self.obstacles[obstacle_name].cells.values():
cells.update(val)
for cell in cells:
field[cell] = value
def set_value_of_obstacle_patch(self, value: float, field: np.ndarray, obstacle_name: str, patch: str):
- for cell in self.obstacle_cells[obstacle_name][patch]:
+ for cell in self.obstacles[obstacle_name].cells[patch]:
field[cell] = value
+
+ def get_obstacles(self):
+ return list(self.obstacles.values())
diff --git a/data_assimilation/TCP_client.py b/data_assimilation/TCP_client.py
index a205afb5..0e643936 100644
--- a/data_assimilation/TCP_client.py
+++ b/data_assimilation/TCP_client.py
@@ -1,11 +1,14 @@
+import datetime
import socket
+import struct
class TCPClient:
def __init__(self):
# Create a TCP/IP socket
self.socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
-
+ timeval = struct.pack('ll', 20, 0)
+ self.socket.setsockopt(socket.SOL_SOCKET, socket.SO_SNDTIMEO, timeval)
# Connect the socket to the port where the server is listening
self.server_address = ('localhost', 7777)
@@ -23,7 +26,8 @@ def send_message(self, message: bin):
# Look for the response
expected_response = "message was received: Rollback done"
response = ''
-
+ self.socket.settimeout(30)
+ print('time', datetime.datetime.now())
while len(response) < len(expected_response):
try:
reply = self.socket.recv(16)
diff --git a/data_assimilation/data_assimilation.py b/data_assimilation/data_assimilation.py
index 9ab9e511..9903d1d3 100644
--- a/data_assimilation/data_assimilation.py
+++ b/data_assimilation/data_assimilation.py
@@ -2,6 +2,7 @@
import os
import re
import struct
+import time
from datetime import datetime
import locale
from typing import List, Dict
@@ -16,7 +17,7 @@
def create_message(t_cur: float, config_file_name: str) -> bin:
- print(f'send message with time step "{t_cur}" and xml "{config_file_name}"')
+ print(f'create message with time step "{t_cur}" and xml "{config_file_name}"')
package = DAPackage(t_cur, config_file_name)
return package.pack()
@@ -69,11 +70,20 @@ def pack(self) -> bin:
class DAFile:
def __init__(self):
self.xml_root = ET.Element('ARTSS')
-
- def write_obstacle_changes(self, list_obstacles: List[Obstacle], obstacle_enabled: bool):
+ self.data_assimilation: ET.SubElement = ET.SubElement(self.xml_root, 'data_assimilation')
+
+ def write_obstacle_changes(self, list_obstacles: List[Obstacle]):
+ """
+ write part of the xml for the obstacle changes. State UNMODIFIED not implemented in ARTSS, use XML and provide the
+ necessary details
+ :param list_obstacles: obstacles to be written out
+ """
# create obstacle part
- #
- #
+ #
+ #
+ #
+ #
+ #
#
#
#
@@ -84,13 +94,13 @@ def write_obstacle_changes(self, list_obstacles: List[Obstacle], obstacle_enable
#
#
#
- #
-
- obstacle_root = ET.SubElement(self.xml_root, 'obstacles',
- enabled='Yes' if obstacle_enabled else 'No')
+ #
+ tag = 'obstacle_changer'
+ ET.SubElement(self.data_assimilation, 'class', {'name': 'ObstacleChanger', 'tag': tag})
+ obstacle_root = ET.SubElement(self.xml_root, tag)
for obstacle in list_obstacles:
single_obstacle = ET.SubElement(obstacle_root, 'obstacle', name=obstacle.name, state=obstacle.state)
- if obstacle.state != 'unmodified':
+ if obstacle.state != 'deleted': # && obstacle.state != 'unmodified'
# convert from Dict[str, float] to Dict[str, str]
geometry = dict(zip(obstacle.geometry.keys(), [str(a) for a in obstacle.geometry.values()]))
ET.SubElement(single_obstacle, 'geometry', geometry)
@@ -104,15 +114,19 @@ def write_obstacle_changes(self, list_obstacles: List[Obstacle], obstacle_enable
def create_obstacle_changes(self, list_obstacles: List[Obstacle], obstacle_enabled: bool):
# create obstacle part
- #
+ #
+ #
+ #
+ #
#
#
#
#
#
- #
- obstacle_root = ET.SubElement(self.xml_root, 'obstacles',
- enabled='Yes' if obstacle_enabled else 'No')
+ #
+ tag = 'obstacle_changer'
+ ET.SubElement(self.data_assimilation, 'class', {'name': 'ObstacleChanger', 'tag': tag})
+ obstacle_root = ET.SubElement(self.xml_root, tag)
for obstacle in list_obstacles:
single_obstacle = ET.SubElement(obstacle_root, 'obstacle', name=obstacle.name)
# convert from Dict[str, float] to Dict[str, str]
@@ -130,23 +144,31 @@ def create_temperature_source_changes(self, source_type: dict,
temperature_source: dict,
random: dict):
# create temperature source part
- #
- # 25000.
- # 1.023415823
- # 40.
- # -3.
- # 0.
- # 1.0
- # 1.5
- # 1.0
- # 5.
- #
- # 0
- # 0.1
- # 1
- #
- #
- source = ET.SubElement(self.xml_root, 'source', type=source_type['type'], dir=source_type['dir'],
+ #
+ #
+ #
+ #
+ #
+ # 25000.
+ # 1.023415823
+ # 40.
+ # -3.
+ # 0.
+ # 1.0
+ # 1.5
+ # 1.0
+ # 5.
+ #
+ # 0
+ # 0.1
+ # 1
+ #
+ #
+ #
+ tag = 'temperature_source'
+ ET.SubElement(self.data_assimilation, 'class', {'name': 'TemperatureSourceChanger', 'tag': tag})
+ source_root = ET.SubElement(self.xml_root, tag)
+ source = ET.SubElement(source_root, 'source', type=source_type['type'], dir=source_type['dir'],
temp_fct=source_type['temp_fct'],
dissipation='Yes' if source_type['dissipation'] == 'Yes' else 'No',
random='Yes' if source_type['random'] == 'Yes' else 'No')
@@ -163,26 +185,30 @@ def create_temperature_source_changes(self, source_type: dict,
if random['custom_seed'] == 'Yes':
ET.SubElement(random_tree, 'seed').text = str(random['seed'])
- def create_config(self, fields: dict, field_file_name=''):
- # create config file. format:
- #
+ def create_field_changes(self, fields: Dict[str, bool], field_file_name=''):
+ # create field changes part
+ #
+ #
+ #
+ #
#
- #
- changed = False
+ #
+ tag: str = 'field_changes'
+ keys = ['u', 'v', 'w', 'p', 'T', 'C']
+ ET.SubElement(self.data_assimilation, 'class', {'name': 'FieldChanger', 'tag': tag})
field_values = {}
- for key in fields.keys():
- if fields[key]:
- field_values[key] = 'Yes'
- changed = True
+ for key in keys:
+ if key in fields.keys():
+ field_values[key] = 'Yes' if fields[key] else 'No'
else:
field_values[key] = 'No'
-
- if changed:
- ET.SubElement(self.xml_root, 'fields_changed', u=field_values['u'], v=field_values['v'],
+ subsection = ET.SubElement(self.xml_root, tag)
+ if 'Yes' in field_values.values():
+ ET.SubElement(subsection, 'fields_changed', u=field_values['u'], v=field_values['v'],
w=field_values['w'], p=field_values['p'],
T=field_values['T'], concentration=field_values['C'], filename=field_file_name)
else:
- ET.SubElement(self.xml_root, 'fields_changed', u=field_values['u'], v=field_values['v'],
+ ET.SubElement(subsection, 'fields_changed', u=field_values['u'], v=field_values['v'],
w=field_values['w'], p=field_values['p'],
T=field_values['T'], concentration=field_values['C'])
@@ -235,16 +261,16 @@ def print_header(self):
@staticmethod
@retry(delay=1, tries=6)
- def get_t_current(path: str = '.') -> float:
- with open(os.path.join(path, '.vis/meta'), 'r') as inp:
+ def get_t_current(path: str = '.', dir_name='.vis') -> float:
+ with open(os.path.join(path, dir_name, 'meta'), 'r') as inp:
t = float([x for x in inp.readlines() if x.startswith('t:')][0][2:])
return t
@staticmethod
@retry(delay=1, tries=6)
- def get_all_time_steps(path: str = '.') -> List[float]:
+ def get_all_time_steps(path: str = '.', dir_name='.vis') -> List[float]:
pattern = re.compile('[0-9]+\.[0-9]{5}e[+|-][0-9]+')
- f = os.listdir(os.path.join(path, '.vis'))
+ f = os.listdir(os.path.join(path, dir_name))
files = []
for p in f:
if pattern.match(p):
@@ -254,8 +280,19 @@ def get_all_time_steps(path: str = '.') -> List[float]:
return files
@staticmethod
- def get_xml_file_name(path: str = '.') -> str:
- fpath = os.path.join(path, '.vis/meta')
+ def get_time_step_artss(t_sensor: float, artss_data_path: str, dt: float, time_back: float = 10) -> [float, float]:
+ files = FieldReader.get_all_time_steps(path=artss_data_path)
+ times = list(filter(lambda x: x >= t_sensor, files))
+ while len(times) == 0:
+ time.sleep(1)
+ files = FieldReader.get_all_time_steps(path=artss_data_path)
+ times = list(filter(lambda x: x >= t_sensor, files))
+ t_revert = ([dt] + list(filter(lambda x: x <= max(0., t_sensor - time_back), files)))[-1]
+ return times[0], t_revert
+
+ @staticmethod
+ def get_xml_file_name(path: str = '.', dir_name='.vis') -> str:
+ fpath = os.path.join(path, dir_name, 'meta')
with open(fpath, 'r') as inp:
xml_file_name = [x for x in inp.readlines() if x.startswith('xml_name:')][0][len('xml_name:'):]
return xml_file_name.strip()
@@ -283,8 +320,10 @@ def read_field_data(self) -> Dict[str, np.ndarray]:
return fields
@staticmethod
- def write_field_data_keys(file_name: str, data: Dict[str, np.ndarray], field_keys: List[str]):
- with h5py.File(file_name, 'w') as out:
+ def write_field_data_keys(file_name: str, data: Dict[str, np.ndarray], field_keys: List[str], path: str = './'):
+ file = os.path.join(path, file_name)
+ print(f"write file {file}")
+ with h5py.File(file, 'w') as out:
for key in field_keys:
if key not in data.keys():
continue
@@ -293,5 +332,5 @@ def write_field_data_keys(file_name: str, data: Dict[str, np.ndarray], field_key
dset = out.create_dataset(key, (len(field),), dtype='d')
dset[:] = field
- def write_field_data(self, file_name: str, data: Dict[str, np.ndarray]):
- FieldReader.write_field_data_keys(file_name=file_name, data=data, field_keys=self.fields)
+ def write_field_data(self, file_name: str, data: Dict[str, np.ndarray], path: str = './'):
+ FieldReader.write_field_data_keys(file_name=file_name, data=data, field_keys=self.fields, path=path)
diff --git a/data_assimilation/downhill_simplex.py b/data_assimilation/downhill_simplex.py
new file mode 100644
index 00000000..d3abffda
--- /dev/null
+++ b/data_assimilation/downhill_simplex.py
@@ -0,0 +1,439 @@
+import math
+import os
+import time
+from pprint import pprint
+from typing import TextIO, Tuple, Dict, List
+
+import numpy as np
+import pandas
+import pandas as pd
+import scipy.optimize as op
+
+import TCP_client
+import data_assimilation
+import fds_utility
+from ARTSS import XML, Domain
+from data_assimilation import FieldReader, DAFile
+from utility import wait_artss, log, write_da_data, kelvin_to_celsius
+
+
+def nelder_mead_special(client: TCP_client,
+ file_da: TextIO, file_debug: TextIO,
+ sensor_times: pandas.Index, devc_info: dict,
+ cur: dict, delta: dict, keys: list,
+ artss_data_path: str, artss: XML, domain: Domain,
+ fds_data: pandas.DataFrame,
+ heat_source: dict,
+ n_iterations: int,
+ parallel=True,
+ precondition=False):
+ t_artss = 0.0
+ t_revert = 0.0
+ t_sensor = 0.0
+
+ def f(x):
+ print(x)
+ print(keys)
+ new_para = {}
+ for i in range(len(keys)):
+ new_para[keys[i]] = x[i]
+ diff_cur, _ = do_rollback(client=client,
+ t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
+ new_para=new_para, heat_source=heat_source.values(),
+ sub_file_name='scipy',
+ artss_data_path=artss_data_path,
+ fds_data=fds_data,
+ devc_info=devc_info,
+ file_da=file_da, file_debug=file_debug)
+ log(f'diff: {diff_cur["T"]}', file_debug)
+ return diff_cur['T']
+
+ def call_back(xk) -> bool:
+ print(f'xk: {xk}')
+ file_debug.write(f'xk: {xk}')
+ return True
+
+ for count, t_sensor in enumerate(sensor_times):
+ pprint(cur)
+
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor,
+ artss_data_path,
+ dt=artss.get_dt(),
+ time_back=6)
+ wait_artss(t_artss, artss_data_path)
+
+ start_time = time.time()
+ log(f't_sensor: {t_sensor} t_artss: {t_artss} t_revert: {t_revert}', file_debug)
+ field_reader = FieldReader(t_artss, path=artss_data_path)
+ file_da.write(f'original data;time_sensor:{t_sensor};time_artss:{t_artss}\n')
+ write_da_data(file_da=file_da, parameters=cur)
+ diff_orig, minima_x = comparison_sensor_simulation_data(devc_info, fds_data, field_reader, t_sensor, file_da)
+
+ file_da.write(f'original: t_artss:{t_artss};t_sensor:{t_sensor};differenceT:{diff_orig["T"]};HRR:{cur["HRR"]};x0:{cur["x0"]}\n')
+
+ log(f'org: {diff_orig["T"]}', file_debug)
+ log(f't_revert: {t_revert}', file_debug)
+
+ if diff_orig['T'] < 1e-5:
+ log(f'skip, difference: {diff_orig["T"]}', file_debug)
+ continue
+
+ if precondition:
+ log('precondition for x0', file_debug)
+ precondition = False
+ safe_keys = keys
+ key_precon = ['x0']
+ keys = key_precon
+ x0 = [cur[x] for x in key_precon]
+ initial_simplex = np.array([x0] * (len(key_precon) + 1))
+ for i in range(len(key_precon)):
+ initial_simplex[i, i] += delta[key_precon[i]]
+ interim_start = time.time()
+ res = op.minimize(f,
+ x0=np.array(x0),
+ callback=call_back,
+ tol=1e-5,
+ method='Nelder-Mead',
+ options={
+ 'maxiter': 5,
+ 'disp': True,
+ 'initial_simplex': initial_simplex
+ })
+ interim_end = time.time()
+ log(f't_artss: {t_artss}; interim time: {interim_end - interim_start}', file_debug)
+ log(f'org: {diff_orig}', file_debug)
+ log(f'res: {res}', file_debug)
+ for i in range(len(keys)):
+ cur[keys[i]] = res.x[i]
+ log(f'res_x (new parameter): {list(res.x)}\n', file_debug)
+ log(f'res_sim (final simplex): {res.final_simplex}\n', file_debug)
+ log(f'res_fun: {res.fun}\n', file_debug)
+ log(f'res_n (number of iterations): {res.nit}\n', file_debug)
+ log(f'res_nfev: {res.nfev}\n', file_debug)
+ log(f'res_suc (exit successfully): {res.success}\n', file_debug)
+ log(f'res_msg (message if exit was unsuccessful): {res.message}\n', file_debug)
+ log(f'final rollback with {cur}', file_debug)
+ diff_cur, _ = do_rollback(client=client,
+ t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
+ new_para=cur, heat_source=heat_source.values(),
+ sub_file_name='final',
+ artss_data_path=artss_data_path,
+ fds_data=fds_data,
+ devc_info=devc_info,
+ file_da=file_da, file_debug=file_debug)
+ file_da.write(f'final: t_artss:{t_artss};t_sensor:{t_sensor};differenceT:{diff_cur["T"]};HRR:{cur["HRR"]};x0:{cur["x0"]}\n')
+ file_da.flush()
+ file_debug.flush()
+ end_time = time.time()
+ log(f'time: {end_time - start_time}', file_debug)
+ keys = safe_keys
+
+ x0 = [cur[x] for x in keys]
+ delta['HRR'] = cur['HRR'] * 0.05,
+ initial_simplex = np.array([x0] * (len(keys) + 1))
+ for i in range(len(keys)):
+ initial_simplex[i, i] += delta[keys[i]]
+ interim_start = time.time()
+ res = op.minimize(f,
+ x0=np.array(x0),
+ callback=call_back,
+ tol=1e-5,
+ method='Nelder-Mead',
+ options={
+ 'maxiter': n_iterations,
+ 'disp': True,
+ 'initial_simplex': initial_simplex
+ })
+ interim_end = time.time()
+ log(f't_artss: {t_artss}; interim time: {interim_end - interim_start}', file_debug)
+ log(f'org: {diff_orig}', file_debug)
+ log(f'res: {res}', file_debug)
+ for i in range(len(keys)):
+ cur[keys[i]] = res.x[i]
+ log(f'res_x (new parameter): {list(res.x)}\n', file_debug)
+ log(f'res_sim (final simplex): {res.final_simplex}\n', file_debug)
+ log(f'res_fun: {res.fun}\n', file_debug)
+ log(f'res_n (number of iterations): {res.nit}\n', file_debug)
+ log(f'res_nfev: {res.nfev}\n', file_debug)
+ log(f'res_suc (exit successfully): {res.success}\n', file_debug)
+ log(f'res_msg (message if exit was unsuccessful): {res.message}\n', file_debug)
+ log(f'final rollback with {cur}', file_debug)
+ diff_cur, _ = do_rollback(client=client,
+ t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
+ new_para=cur, heat_source=heat_source.values(),
+ sub_file_name='final',
+ artss_data_path=artss_data_path,
+ fds_data=fds_data,
+ devc_info=devc_info,
+ file_da=file_da, file_debug=file_debug)
+ file_da.write(f'final: t_artss:{t_artss};t_sensor:{t_sensor};differenceT:{diff_cur["T"]};HRR:{cur["HRR"]};x0:{cur["x0"]}\n')
+ file_da.flush()
+ file_debug.flush()
+ end_time = time.time()
+ log(f'time: {end_time - start_time}', file_debug)
+
+
+def opt_scipy(client: TCP_client,
+ file_da: TextIO, file_debug: TextIO,
+ sensor_times: List[float], devc_info: dict,
+ cur: dict, delta: dict, keys: list,
+ artss_data_path: str, artss: XML, domain: Domain,
+ fds_data: pandas.DataFrame,
+ heat_source: dict,
+ n_iterations: int,
+ offset: float = 0,
+ parallel=True):
+ t_artss = 0.0
+ t_revert = 0.0
+ t_sensor = 0.0
+
+ bounds = op.Bounds(
+ lb=np.asarray([0, domain.domain_param['X1']]),
+ ub=np.asarray([np.inf, domain.domain_param['X2']])
+ )
+
+ def f(x):
+ print(x)
+ print(keys)
+ new_para = {}
+ for i in range(len(keys)):
+ new_para[keys[i]] = x[i]
+ diff_cur, _ = do_rollback(client=client,
+ t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
+ new_para=new_para, heat_source=heat_source.values(),
+ sub_file_name='scipy',
+ artss_data_path=artss_data_path,
+ fds_data=fds_data,
+ devc_info=devc_info,
+ file_da=file_da, file_debug=file_debug)
+ log(f'diff: {diff_cur["T"]}', file_debug)
+ return diff_cur['T']
+
+ def call_back(xk) -> bool:
+ print(f'xk: {xk}')
+ file_debug.write(f'xk: {xk}')
+ return True
+
+ iterations = np.ones(len(sensor_times)) * n_iterations
+ # iterations[0] = 15
+ for count, t_sensor in enumerate(sensor_times):
+ t_projected = t_sensor - offset + artss.get_dt()
+ pprint(cur)
+
+ wait_artss(t_projected, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_projected,
+ artss_data_path,
+ dt=artss.get_dt(),
+ time_back=6)
+ wait_artss(t_artss, artss_data_path)
+ if count == 0:
+ t_revert = 0.2
+ start_time = time.time()
+ log(f't_sensor: {t_sensor} sensor_time: {t_projected} t_artss: {t_artss} t_revert: {t_revert}', file_debug)
+ field_reader = FieldReader(t_artss, path=artss_data_path)
+ file_da.write(f'original data;time_sensor:{t_sensor};time_projected:{t_projected};time_artss:{t_artss}\n')
+ write_da_data(file_da=file_da, parameters=cur)
+ diff_orig, minima_x = comparison_sensor_simulation_data(devc_info, fds_data, field_reader, t_sensor, file_da)
+
+ file_da.write(f'original: t_artss:{t_artss};t_sensor:{t_sensor};time_projected:{t_projected};differenceT:{diff_orig["T"]};HRR:{cur["HRR"]};x0:{cur["x0"]};z0:{cur["z0"]}\n')
+
+ log(f'org: {diff_orig["T"]}', file_debug)
+ log(f't_revert: {t_revert}', file_debug)
+
+ if diff_orig['T'] < 1e-5:
+ log(f'skip, difference: {diff_orig["T"]}', file_debug)
+ continue
+ x0 = [cur[x] for x in keys]
+ delta['HRR'] = cur['HRR'] * 0.05,
+ initial_simplex = np.array([x0] * (len(keys) + 1))
+ for i in range(len(keys)):
+ initial_simplex[i, i] += delta[keys[i]]
+ interim_start = time.time()
+ res = op.minimize(f,
+ x0=np.array(x0),
+ callback=call_back,
+ tol=1e-5,
+ method='Nelder-Mead',
+ bounds=bounds,
+ options={
+ 'maxiter': iterations[count],
+ 'disp': True,
+ 'initial_simplex': initial_simplex
+ })
+ interim_end = time.time()
+ log(f't_artss: {t_artss}; interim time: {interim_end - interim_start}', file_debug)
+ log(f'org: {diff_orig}', file_debug)
+ log(f'res: {res}', file_debug)
+ for i in range(len(keys)):
+ cur[keys[i]] = res.x[i]
+ log(f'res_x (new parameter): {list(res.x)}\n', file_debug)
+ log(f'res_sim (final simplex): {res.final_simplex}\n', file_debug)
+ log(f'res_fun: {res.fun}\n', file_debug)
+ log(f'res_n (number of iterations): {res.nit}\n', file_debug)
+ log(f'res_nfev: {res.nfev}\n', file_debug)
+ log(f'res_suc (exit successfully): {res.success}\n', file_debug)
+ log(f'res_msg (message if exit was unsuccessful): {res.message}\n', file_debug)
+ log(f'final rollback with {cur}', file_debug)
+ diff_cur, _ = do_rollback(client=client,
+ t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
+ new_para=cur, heat_source=heat_source.values(),
+ sub_file_name='final',
+ artss_data_path=artss_data_path,
+ fds_data=fds_data,
+ devc_info=devc_info,
+ file_da=file_da, file_debug=file_debug)
+ file_da.write(f'final: t_artss:{t_artss};t_sensor:{t_sensor};time_projected:{t_projected};differenceT:{diff_cur["T"]};HRR:{cur["HRR"]};x0:{cur["x0"]};z0:{cur["z0"]}\n')
+ file_da.flush()
+ file_debug.flush()
+ end_time = time.time()
+ log(f'time: {end_time - start_time}', file_debug)
+
+
+def do_rollback(client: TCP_client,
+ t_sensor, t_artss, t_revert,
+ new_para: dict, heat_source: [dict, dict, dict],
+ sub_file_name: str,
+ artss_data_path: str,
+ fds_data,
+ devc_info: dict,
+ file_da: TextIO, file_debug: TextIO) -> [dict, float]:
+ config_file_name, _ = write_changes_xml(
+ new_para,
+ heat_source,
+ f'{t_artss}_{sub_file_name}',
+ file_da=file_da,
+ path=artss_data_path
+ )
+ print(f'make adjustment with {new_para}')
+ file_debug.write(f'make adjustment with: {new_para}\n')
+ write_da_data(file_da=file_da, parameters=new_para)
+ client.send_message(data_assimilation.create_message(t_revert, config_file_name))
+ wait_artss(t_artss, artss_data_path)
+
+ field_reader = FieldReader(t_artss, path=artss_data_path)
+ diff_cur, min_pos_x = comparison_sensor_simulation_data(
+ devc_info,
+ fds_data,
+ field_reader,
+ t_sensor,
+ file_da
+ )
+ return diff_cur, min_pos_x
+
+
+def comparison_sensor_simulation_data(devc_info: dict, sensor_data: pd.DataFrame, field_reader: FieldReader,
+ t_sensor: float, file_da: TextIO) -> Tuple[Dict[str, float], float]:
+ diff: dict = {'T': [], 'C': []}
+ fields_sim = field_reader.read_field_data()
+ min_pos_x: float = 0
+ min_sensor_val = math.inf
+ for key in devc_info:
+ # sensor data
+ type_sensor: str = devc_info[key]['type']
+ index_sensor: int = devc_info[key]['index']
+ value_sensor: float = sensor_data[key][t_sensor]
+
+ value_sim = fields_sim[type_sensor][index_sensor]
+ difference = kelvin_to_celsius(value_sim) - value_sensor
+ if abs(difference) > 2:
+ diff[type_sensor].append(difference)
+ file_da.write(
+ f'sensor:{key};time_sensor:{t_sensor};time_artss:{field_reader.t};sensor_val:{value_sensor};artss_val:{value_sim};diff:{difference};considered:{abs(difference) > 2};position:{devc_info[key]["XYZ"]}\n')
+ file_da.flush()
+
+ if difference < min_sensor_val:
+ min_sensor_val = difference
+ min_pos_x = devc_info[key]['XYZ'][0]
+ print(f'diff: {diff["T"]}')
+ result = {}
+ for key in diff:
+ if len(diff[key]) == 0:
+ result[key] = 0
+ file_da.write(f'differences:{key}:{result[key]};no_of_sensors:{len(diff[key])}\n')
+ continue
+ result[key] = np.sqrt(sum(np.array(diff[key]) ** 2))
+ file_da.write(f'differences:{key}:{result[key]};no_of_sensors:{len(diff[key])}\n')
+ return result, min_pos_x
+
+
+def write_changes_xml(change: dict, source: list, file_name: str, file_da: TextIO, path='.') -> [str, str]:
+ source_type, temperature_source, random = data_assimilation.change_heat_source(*source,
+ changes={'source_type': {},
+ 'temperature_source': change,
+ 'random': {}})
+ write_da_data(file_da=file_da, parameters=temperature_source)
+ da = DAFile()
+ da.create_temperature_source_changes(
+ source_type=source_type,
+ temperature_source=temperature_source,
+ random=random)
+
+ config_file_name = f'config_{file_name}.xml'
+ config_file_path = os.path.join(path, config_file_name)
+ da.write_xml(config_file_path)
+ return config_file_name, config_file_path
+
+
+def start(fds_data_path: str, fds_input_file_name: str, artss_data_path: str, parallel=True, port=7777):
+ artss_data_path = os.path.abspath(artss_data_path)
+ client = TCP_client.TCPClient()
+ client.set_server_address('localhost', port)
+ # client.set_server_address('172.17.0.2', 7777)
+ client.connect()
+
+ xml = XML(FieldReader.get_xml_file_name(artss_data_path), path=artss_data_path)
+ xml.read_xml()
+ domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles,
+ enable_computational_domain=xml.computational_domain)
+ domain.print_info()
+ domain.print_debug()
+
+ devices, fds_data = fds_utility.read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devc_info_temperature = devices['temperature']
+ print('devices:')
+ pprint(devc_info_temperature)
+
+ sensor_times = fds_data.index[3:]
+ print('sensor times:', fds_data.index)
+
+ heat_source = xml.get_temperature_source()
+
+ file_da = open(os.path.join(artss_data_path, 'da_details.csv'), 'w')
+ file_debug = open(os.path.join(artss_data_path, 'da_debug_details.dat'), 'w')
+ pprint(devc_info_temperature, stream=file_debug)
+
+ delta = {
+ 'HRR': float(heat_source['temperature_source']['HRR']) * 0.05,
+ 'x0': 1.5,
+ 'y0': domain.domain_param['dy'] * 1,
+ 'z0': 0.5
+ }
+
+ cur = {
+ 'HRR': float(heat_source['temperature_source']['HRR']),
+ 'x0': float(heat_source['temperature_source']['x0']),
+ 'y0': float(heat_source['temperature_source']['y0']),
+ 'z0': float(heat_source['temperature_source']['z0'])
+ }
+
+ keys = ['HRR', 'x0', 'z0']
+
+ log(f'cur: {cur}', file_debug)
+ log(f'delta: {delta}', file_debug)
+ log(f'keys: {keys}', file_debug)
+ # nelder_mead_special(client=client, file_da=file_da, file_debug=file_debug,
+ # sensor_times=sensor_times,
+ # devc_info=devc_info_temperature, fds_data=fds_data,
+ # artss_data_path=artss_data_path,
+ # domain=domain, heat_source=heat_source,
+ # cur=cur, delta=delta, keys=keys, n_iterations=5, artss=xml,
+ # parallel=parallel,
+ # precondition=True)
+ opt_scipy(client=client, file_da=file_da, file_debug=file_debug,
+ sensor_times=sensor_times,
+ devc_info=devc_info_temperature, fds_data=fds_data,
+ artss_data_path=artss_data_path,
+ domain=domain, heat_source=heat_source,
+ cur=cur, delta=delta, keys=keys, n_iterations=5, artss=xml,
+ parallel=parallel)
diff --git a/data_assimilation/example/test_door.xml b/data_assimilation/example/test_door.xml
index 5b7e9fe1..db3b0819 100644
--- a/data_assimilation/example/test_door.xml
+++ b/data_assimilation/example/test_door.xml
@@ -77,7 +77,7 @@
-
+
0.05
diff --git a/data_assimilation/example/tunnel_30.xml b/data_assimilation/example/tunnel_30.xml
index ed6916e6..320e8473 100644
--- a/data_assimilation/example/tunnel_30.xml
+++ b/data_assimilation/example/tunnel_30.xml
@@ -53,9 +53,9 @@
30
0.
3.
- 0.5
- 0.1
- 0.5
+ 1.5
+ 3
+ 2.5
5.
diff --git a/data_assimilation/fds_utility.py b/data_assimilation/fds_utility.py
new file mode 100644
index 00000000..34ee1292
--- /dev/null
+++ b/data_assimilation/fds_utility.py
@@ -0,0 +1,154 @@
+import os
+from typing import Dict, List
+
+import fdsreader
+import pandas
+import pandas as pd
+
+from ARTSS import Domain
+
+
+def read_fds_data(input_path: str, input_file_name: str, artss: Domain) -> [Dict[str, Dict[str, any]], pd.DataFrame]:
+ full_name = os.path.join(input_path, input_file_name)
+ str_devc = read_devc_from_input_file(full_name + '.fds')
+ dict_devc = parse_devc(str_devc)
+ chid_file_name = get_chid_from_input_file(full_name + '.fds')
+
+ fds_data = read_devc_from_csv_file(os.path.join(input_path, chid_file_name + '_devc.csv'))
+
+ devices: Dict[str, Dict[str, any]] = {}
+ devc_temperature: Dict[str, any] = {}
+ devc_velocity: Dict[str, any] = {}
+ devc_pressure: Dict[str, any] = {}
+ devc_visibility: Dict[str, any] = {}
+ for d in dict_devc:
+ dict_devc[d]['type']: str = 'T'
+
+ ijk = artss.get_ijk_from_xyz(dict_devc[d]['XYZ'][0], dict_devc[d]['XYZ'][2], dict_devc[d]['XYZ'][1])
+ dict_devc[d]['index']: int = artss.get_index(*ijk)
+ d_lower = d.lower()
+ if d_lower.startswith('temperatur'):
+ devc_temperature[d] = dict_devc[d]
+ elif d_lower.startswith('velo'):
+ devc_velocity[d] = dict_devc[d]
+ elif d_lower.startswith('vis'):
+ devc_visibility[d] = dict_devc[d]
+ elif d_lower.startswith('pressure'):
+ devc_pressure[d] = dict_devc[d]
+
+ devices['temperature'] = devc_temperature
+ devices['velocity'] = devc_velocity
+ devices['pressure'] = devc_pressure
+ devices['visibility'] = devc_visibility
+ return devices, fds_data
+
+
+def read_devc_from_csv_file(file_name: str) -> pd.DataFrame:
+ df = pd.read_csv(file_name, skiprows=1)
+ df.index = df.iloc[:, 0]
+ df.drop('Time', axis=1, inplace=True)
+ return df
+
+
+def get_chid_from_input_file(input_file: str) -> str:
+ with open(input_file, 'r') as inp:
+ for line in inp:
+ if line.startswith('&HEAD'):
+ line = line[len('&HEAD'):].strip()
+ tmp_arr = line.split('=')
+ index = tmp_arr.index('CHID')
+ chid_str = tmp_arr[index + 1].strip()
+ index_end = chid_str[1:].index("'")
+ chid = chid_str[1:index_end + 1]
+ break
+ return chid
+
+
+def parse_devc(str_devc: list, start='&DEVC') -> dict:
+ devc = {}
+ for line in str_devc:
+ line = line[len(start):-1]
+ line.strip()
+ tmp_arr = line.split(',')
+ parameters = []
+ for i in tmp_arr:
+ if type(i) == str and '=' in i:
+ parameters.append(i)
+ else:
+ parameters[-1] += f',{i}'
+
+ param_dict = {}
+ for p in parameters:
+ tmp_arr = p.split('=')
+ param_dict[tmp_arr[0].strip()] = tmp_arr[1].strip().strip("'").strip('"')
+ id = param_dict.pop('ID')
+ if 'XYZ' in param_dict.keys():
+ param_dict['XYZ'] = [float(i) for i in param_dict['XYZ'].split(',')]
+ devc[id] = param_dict
+ return devc
+
+
+def read_devc_from_input_file(fds_file_name: str, start="&DEVC", end="/") -> list:
+ """
+ Reads FDS input out of a list of string and collects namelist definitions
+ of a given type, e.g. "&DEVC ".
+ Reads all definitions line by line, without processing.
+
+ :param fds_file_name: List of string, of a single namelist definition.
+ :param start: String, indicate the namelist, default "&DEVC ".
+ :param end: String, indicate the end of the namelist, default "/".
+ :return: list with lines
+ inspired by Tristan Hehnen
+ """
+
+ namelist_lines = []
+
+ collect_line = False
+ fds_input_content = open(fds_file_name, 'r')
+ for line in fds_input_content:
+ line = line.strip()
+ if collect_line is True:
+ namelist_lines[-1] += f' {line}'
+ elif line.startswith(start):
+ namelist_lines.append(line)
+ collect_line = True
+ collect_line = collect_line and (end not in line)
+ fds_input_content.close()
+ return namelist_lines
+
+
+def read_fds_file(data_path: str, artss: Domain) -> pd.DataFrame:
+ data = fdsreader.Simulation(data_path)
+ devc = data.devices
+ return_val: pd.DataFrame = pd.DataFrame([], index=devc['Time'].data)
+ return_val.index.names = ['Time']
+ return_val.columns.names = ['Position of Thermocouple']
+ for key in devc:
+ if key.startswith('Thermocouple'):
+ if devc[key] is not None:
+ print(key, devc[key].xyz, *devc[key].xyz, devc[key].xyz[1])
+ # be careful, in FDS the third parameter indicates the height, but in ARTSS it is the second
+ i, j, k = artss.get_ijk_from_xyz(devc[key].xyz[0], devc[key].xyz[2], devc[key].xyz[1])
+ index = artss.get_index(i, j, k)
+ print(f'index: {index} of {i}|{j}|{k}')
+ if index in return_val:
+ print(f'index {index} already exists, take average')
+ tmp = (devc[key].data + return_val[index]) / 2.
+ return_val[index] = tmp
+ else:
+ return_val[index] = devc[key].data
+ else:
+ print(f'{key} is None, skipped')
+ return return_val
+
+
+def get_starting_time(fds_data : pd.DataFrame, threshold: float, keys: List[str]):
+ columns = []
+ for i in fds_data:
+ for key in keys:
+ if key in i.lower():
+ columns.append(i)
+ data_extract = fds_data[columns]
+ for index, row in data_extract.iterrows():
+ if any(row > threshold):
+ return index
diff --git a/data_assimilation/full_corridor/full_corridor.xml b/data_assimilation/full_corridor/full_corridor.xml
new file mode 100644
index 00000000..54025130
--- /dev/null
+++ b/data_assimilation/full_corridor/full_corridor.xml
@@ -0,0 +1,224 @@
+
+
+
+ 1000.
+ 0.05
+ 3.1e-5
+ 3.34e-3
+ -9.81
+ 1
+ 4.2e-5
+
+
+
+
+
+
+ 100
+ 1e-07
+ 1
+
+
+ 0.2
+
+
+ 293.15
+
+
+ 1
+ 2
+ 100
+ 1e-07
+ 4
+
+ 100
+ 1e-07
+ 0.6666666667
+
+
+
+
+
+
+ 100
+ 1e-07
+ 1
+
+
+ 0.5
+
+
+ 400
+ 1.
+ 1.5
+ 0.1
+ 2
+ 2
+ 0.1
+ 1.5
+ 5.
+
+
+
+
+
+
+
+ 0
+ 20
+ 0
+ 2.4
+ 0
+ 12.4
+ 160
+ 48
+ 248
+
+
+
+
+
+
+ 0.05
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 293.15
+
+ 0
+ 0.1
+ 1
+
+
+
+
+ 1
+
+
+
+
+
diff --git a/data_assimilation/gradient_based_optimisation.py b/data_assimilation/gradient_based_optimisation.py
index 401ea06c..e067ab2a 100644
--- a/data_assimilation/gradient_based_optimisation.py
+++ b/data_assimilation/gradient_based_optimisation.py
@@ -1,9 +1,8 @@
import math
import os
import time
-from typing import Dict, TextIO
+from typing import Dict, TextIO, Tuple
-import fdsreader
import multiprocessing
import numpy as np
import pandas
@@ -13,105 +12,59 @@
from copy import copy
from pprint import pprint
-from tqdm import tqdm
-
import ARTSS
import TCP_client
import data_assimilation
+import fds_utility
from ARTSS import Domain, XML
from data_assimilation import FieldReader, create_message, DAFile
+from utility import wait_artss, log, write_da_data, kelvin_to_celsius
-def write_da_data(file_da: TextIO, parameters: dict):
- for key in parameters:
- file_da.write(f';{key}:{parameters[key]}')
- file_da.write('\n')
-
-
-def map_minima(client, artss_data_path, cur, delta, sensor_times, devc_info_thermocouple, devc_info_temperature,
- fds_data, source_type, temperature_source, random, file_da, cwd, artss, xml):
+def map_minima(client, artss_data_path, cur, sensor_times, devc_info_thermocouple, devc_info_temperature,
+ fds_data: pd.DataFrame, heat_source, file_da: TextIO, cwd: str, xml: ARTSS.XML):
def f(t_end, param):
- t_artss, t_revert = get_time_step_artss(t_end, artss_data_path, xml.get_dt())
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_end, artss_data_path, xml.get_dt())
config_file_name, _ = write_changes_xml(
param,
- [source_type, temperature_source, random],
+ [heat_source['source_type'], heat_source['temperature_source'], heat_source['random']],
f'{t_artss}_f',
file_da=file_da,
path=cwd
)
client.send_message(create_message(t_revert, config_file_name))
- wait_artss(t_end, artss_data_path, artss)
+ wait_artss(t_artss, artss_data_path)
field_reader = FieldReader(t_artss, path=artss_data_path)
- ret, _ = comparison_sensor_simulation_data(
- devc_info_thermocouple,
- fds_data,
- field_reader,
- t_artss,
- t_end
- )
+ ret, _ = comparison_sensor_simulation_data(devc_info=devc_info_temperature, sensor_data=fds_data,
+ field_reader=field_reader, t_sensor=t, file_da=file_da)
return ret['T']
print('org: {f(t_artss, t_revert, dict())}')
- t = sensor_times[2]
+ t = sensor_times[3]
pprint(cur)
- wait_artss(t, artss_data_path, artss)
+ wait_artss(t, artss_data_path)
- t_artss, t_revert = get_time_step_artss(t, artss_data_path, xml.get_dt())
+ t_artss, t_revert = FieldReader.get_time_step_artss(t, artss_data_path, xml.get_dt())
pprint((t, t_artss, t_revert))
+ wait_artss(t_artss, artss_data_path)
+
field_reader = FieldReader(t_artss, path=artss_data_path)
- diff_orig, _ = comparison_sensor_simulation_data(devc_info_temperature, fds_data, field_reader, t_artss, t)
+ diff_orig, _ = comparison_sensor_simulation_data(devc_info=devc_info_temperature, sensor_data=fds_data,
+ field_reader=field_reader, t_sensor=t, file_da=file_da)
print(f'org: {diff_orig["T"]}')
print(t_revert)
- with open('map_2.csv', 'w') as out:
- for x in range(20, 85, 10):
- for y in [delta['y0'] * t for t in range(8)]:
- ret = f(t, {'x0': x, 'y0': y})
- print(f'map: {x},{y},{ret}')
- out.write(f'{x},{y},{ret}\n')
+ with open(os.path.join(cwd, 'map_2.csv'), 'w') as out:
+ for hrr in range(4000, 7000, 200):
+ ret = f(t, {'HRR': hrr})
+ print(f'map: {hrr},{ret}')
+ out.write(f'{hrr},{ret}\n')
return
-def start_new_instance(output_file: str, directory: str, artss_exe_path: str,
- artss_exe: str = 'artss_data_assimilation_serial'):
- cwd = os.getcwd()
- os.chdir(directory)
- exe_command = f'mpirun --np 2 {os.path.join(artss_exe_path, artss_exe)} {output_file}'
- print(exe_command)
- os.system(exe_command)
- os.chdir(cwd)
-
-
-def create_start_xml(change: dict, input_file: str, output_file: str, t_revert: float, file: str, t_end: float,
- directory: str, artss_path: str, file_debug: TextIO, dir_name: str = '') -> [str, str]:
- f = os.path.join('..', file)
- f = f.replace("/", "\/")
-
- command = ''
- name = ''
- for key in change:
- command += f' --{key} {change[key]}'
- name += f'{key}_'
- name = name[:-1]
-
- if not dir_name == '':
- name = dir_name
- output_dir = os.path.join(directory, '.vis', name)
- if not os.path.exists(output_dir):
- os.makedirs(output_dir)
-
- out_file = os.path.join(output_dir, output_file)
- log(f'out_file: {out_file}', file_debug)
- log(f'out_dir: {output_dir}', file_debug)
- command = f'{os.path.join(artss_path, "tools", "change_xml.sh")} -i {input_file} -o {out_file} --da {t_revert} "{f}" 7779 --tend {t_end} --loglevel off' + command
- log(f'{command}', file_debug)
- os.system(command)
- return out_file, output_dir
-
-
def calc_diff(t_artss: float, t_sensor: float,
data_path: str, fds_data: pd.DataFrame, devc_info: dict,
file_da: TextIO) -> [Dict[str, float], float]:
@@ -131,7 +84,7 @@ def do_rollback(client: TCP_client,
t_sensor, t_artss, t_revert,
new_para: dict, heat_source: [dict, dict, dict],
sub_file_name: str,
- artss_data_path: str, artss: XML,
+ artss_data_path: str,
fds_data,
devc_info: dict,
file_da: TextIO, file_debug: TextIO) -> [dict, float]:
@@ -146,7 +99,7 @@ def do_rollback(client: TCP_client,
file_debug.write(f'make adjustment with: {new_para}\n')
write_da_data(file_da=file_da, parameters=new_para)
client.send_message(create_message(t_revert, config_file_name))
- wait_artss(t_sensor, artss_data_path, artss)
+ wait_artss(t_artss, artss_data_path)
field_reader = FieldReader(t_artss, path=artss_data_path)
diff_cur, min_pos_x = comparison_sensor_simulation_data(
@@ -159,11 +112,6 @@ def do_rollback(client: TCP_client,
return diff_cur, min_pos_x
-def log(message: str, file_debug: TextIO):
- print(message)
- file_debug.write(f'{message}\n')
-
-
def continuous_gradient_parallel(client: TCP_client,
file_da: TextIO, file_debug: TextIO,
sensor_times: pandas.Index, devc_info: dict,
@@ -178,11 +126,11 @@ def continuous_gradient_parallel(client: TCP_client,
artss_exe_path = os.path.join(artss_path, 'build', 'bin')
xml_input_file = os.path.join(artss_data_path, FieldReader.get_xml_file_name(path=artss_data_path))
n = max(1, n_iterations)
- for t_sensor in sensor_times[3:]:
+ for t_sensor in sensor_times:
pprint(cur)
- wait_artss(t_sensor, artss_data_path, artss)
-
- t_artss, t_revert = get_time_step_artss(t_sensor, artss_data_path, dt=artss.get_dt(), time_back=6)
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=artss.get_dt(), time_back=6)
+ wait_artss(t_artss, artss_data_path)
log(f't_sensor: {t_sensor} t_artss: {t_artss} t_revert: {t_revert}', file_debug)
field_reader = FieldReader(t_artss, path=artss_data_path)
@@ -205,7 +153,7 @@ def continuous_gradient_parallel(client: TCP_client,
t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
new_para={'x0': cur['x0']}, heat_source=heat_source.values(),
sub_file_name='precondition',
- artss_data_path=artss_data_path, artss=artss,
+ artss_data_path=artss_data_path,
fds_data=fds_data,
devc_info=devc_info,
file_da=file_da, file_debug=file_debug)
@@ -216,15 +164,15 @@ def continuous_gradient_parallel(client: TCP_client,
directories = {}
counter = 1
for p in keys:
- output_file, output_dir = create_start_xml(change={p: cur[p] + delta[p]},
- input_file=xml_input_file, output_file=f'config_{p}.xml',
- t_revert=t_revert, t_end=t_artss,
- directory=artss_data_path,
- file=f'{t_revert:.5e}',
- artss_path=artss_path,
- file_debug=file_debug)
+ output_file, output_dir = ARTSS.create_start_xml(change={p: cur[p] + delta[p]},
+ input_file=xml_input_file, output_file=f'config_{p}.xml',
+ t_revert=t_revert, t_end=t_artss,
+ directory=artss_data_path,
+ file=f'{t_revert:.5e}',
+ artss_path=artss_path,
+ file_debug=file_debug)
directories[p] = output_dir
- job = multiprocessing.Process(target=start_new_instance, args=(
+ job = multiprocessing.Process(target=ARTSS.start_new_instance, args=(
output_file, output_dir, artss_exe_path, 'artss_data_assimilation_serial'))
job.start()
counter += 1
@@ -305,7 +253,7 @@ def continuous_gradient_parallel(client: TCP_client,
t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
new_para=new_para, heat_source=heat_source.values(),
sub_file_name=n,
- artss_data_path=artss_data_path, artss=artss,
+ artss_data_path=artss_data_path,
fds_data=fds_data,
devc_info=devc_info,
file_da=file_da, file_debug=file_debug)
@@ -341,12 +289,11 @@ def continuous_gradient(client: TCP_client,
heat_source: dict,
n_iterations: int,
precondition: bool = False):
- n = max(1, n_iterations)
- for t_sensor in sensor_times[2:]:
+ for t_sensor in sensor_times:
pprint(cur)
- wait_artss(t_sensor, artss_data_path, artss)
-
- t_artss, t_revert = get_time_step_artss(t_sensor, artss_data_path, dt=artss.get_dt(), time_back=6)
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=artss.get_dt(), time_back=6)
+ wait_artss(t_artss, artss_data_path)
log(f't_sensor: {t_sensor} t_artss: {t_artss} t_revert: {t_revert}', file_debug)
field_reader = FieldReader(t_artss, path=artss_data_path)
@@ -357,7 +304,8 @@ def continuous_gradient(client: TCP_client,
log(f'org: {diff_orig["T"]}', file_debug)
log(f't revert: {t_revert}', file_debug)
- if sum(diff_orig.values()) < 1e-5:
+ if diff_orig['T'] < 1e-5:
+ log(f'skip, difference: {diff_orig["T"]}', file_debug)
continue
if precondition:
@@ -365,27 +313,32 @@ def continuous_gradient(client: TCP_client,
cur['x0'] = minima_x
heat_source['temperature_source']['x0'] = cur['x0']
file_da.write(f'preconditioning;')
+ start_time = time.time()
diff_orig, _ = do_rollback(client=client,
t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
new_para={'x0': cur['x0']}, heat_source=heat_source.values(),
sub_file_name='precondition',
- artss_data_path=artss_data_path, artss=artss,
+ artss_data_path=artss_data_path,
fds_data=fds_data,
devc_info=devc_info,
file_da=file_da, file_debug=file_debug)
-
+ end_time = time.time()
+ file_da.write(f'time: {end_time - start_time}\n')
nabla = {}
for p in keys:
file_da.write(f'calc_nabla;')
+ start_time = time.time()
diff_cur, _ = do_rollback(client=client,
t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
new_para={p: cur[p] + delta[p]}, heat_source=heat_source.values(),
sub_file_name=p,
- artss_data_path=artss_data_path, artss=artss,
+ artss_data_path=artss_data_path,
fds_data=fds_data,
devc_info=devc_info,
file_da=file_da, file_debug=file_debug)
+ end_time = time.time()
+ file_da.write(f'time: {end_time - start_time}\n')
# calc new nabla
nabla[p] = (diff_cur['T'] - diff_orig['T']) / delta[p]
log(f'cur: {diff_cur["T"]}', file_debug)
@@ -432,8 +385,9 @@ def continuous_gradient(client: TCP_client,
alpha_min = min(alpha_min, restrictions[k])
alpha = 0.9 * alpha_min
log(f'mins: 1.0, {restrictions.values()}', file_debug)
- log(f'alpha: {alpha}', file_debug)
+ n = max(1, n_iterations)
while n > 0:
+ log(f'alpha: {alpha}, n: {n}', file_debug)
x = np.asarray([cur[x] for x in keys])
new_para = x + (alpha * d)
new_para = {
@@ -441,23 +395,30 @@ def continuous_gradient(client: TCP_client,
for i, p in enumerate(keys)
}
pprint(new_para)
+ pprint(new_para, stream=file_debug)
file_da.write(f'iterate:{n};')
+ start_time = time.time()
diff_cur, _ = do_rollback(client=client,
t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
new_para=new_para, heat_source=heat_source.values(),
sub_file_name=n,
- artss_data_path=artss_data_path, artss=artss,
+ artss_data_path=artss_data_path,
fds_data=fds_data,
devc_info=devc_info,
file_da=file_da, file_debug=file_debug)
- log(f'conditional statement {diff_orig["T"] + np.dot(alpha * sigma * nabla, d)}', file_debug)
+ end_time = time.time()
+ file_da.write(f'time: {end_time - start_time}\n')
+ log(f'conditional statement {diff_cur["T"]} < {diff_orig["T"] + np.dot(alpha * sigma * nabla, d)}',
+ file_debug)
if diff_cur['T'] < diff_orig['T'] + np.dot(alpha * sigma * nabla, d):
log(f'found alpha: {alpha}', file_debug)
log(f'found x_k+1: {new_para}', file_debug)
log(f"al1: {np.dot(alpha * sigma * nabla, d)}", file_debug)
log(f"als: {diff_cur['T']} < {diff_orig['T'] + np.dot(alpha * sigma * nabla, d)}", file_debug)
cur = copy(new_para)
+ file_da.flush()
+ file_debug.flush()
break
alpha = 0.7 * alpha
@@ -471,167 +432,8 @@ def continuous_gradient(client: TCP_client,
log(f'alpha: {alpha}', file_debug)
log(f'using cur: {cur}', file_debug)
-
-
-def one_time_gradient_parallel(client: TCP_client,
- file_da: TextIO, file_debug: TextIO,
- sensor_times: pandas.Index, devc_info: dict,
- cur: dict, delta: dict, keys: list,
- artss_path: str, artss_data_path: str,
- artss: XML, domain: Domain,
- fds_data: pandas.DataFrame,
- heat_source: dict,
- n_iterations: int,
- precondition: bool = False):
- no_cpus = multiprocessing.cpu_count()
- artss_exe_path = os.path.join(artss_path, 'build', 'bin')
- xml_input_file = os.path.join(artss_data_path, FieldReader.get_xml_file_name(path=artss_data_path))
- for t_sensor in sensor_times[2:]:
- pprint(cur)
- wait_artss(t_sensor, artss_data_path, artss)
-
- t_artss, t_revert = get_time_step_artss(t_sensor, artss_data_path, dt=artss.get_dt(), time_back=6)
-
- log(f't_sensor: {t_sensor} t_artss: {t_artss} t_revert: {t_revert}', file_debug)
- field_reader = FieldReader(t_artss, path=artss_data_path)
- file_da.write(f'original data;time_sensor:{t_sensor};time_artss:{t_artss}\n')
- write_da_data(file_da=file_da, parameters=cur)
- diff_orig, minima_x = comparison_sensor_simulation_data(devc_info, fds_data, field_reader, t_sensor, file_da)
-
- log(f'org: {diff_orig["T"]}', file_debug)
- log(f't revert: {t_revert}', file_debug)
-
- if sum(diff_orig.values()) < 1e-5:
- continue
-
- if precondition:
- precondition = False
- cur['x0'] = minima_x
- heat_source['temperature_source']['x0'] = cur['x0']
- file_da.write(f'preconditioning;')
- diff_orig, _ = do_rollback(client=client,
- t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
- new_para={'x0': cur['x0']}, heat_source=heat_source.values(),
- sub_file_name='precondition',
- artss_data_path=artss_data_path, artss=artss,
- fds_data=fds_data,
- devc_info=devc_info,
- file_da=file_da, file_debug=file_debug)
-
- nabla = {}
-
- jobs = {}
- directories = {}
- counter = 1
- for p in keys:
- file_da.write(f'calc_nabla;')
- output_file, output_dir = create_start_xml(change={p: cur[p] + delta[p]},
- input_file=xml_input_file, output_file=f'config_{p}.xml',
- t_revert=t_revert, t_end=t_artss,
- directory=artss_data_path,
- file=f'{t_revert:.5e}',
- artss_path=artss_path,
- file_debug=file_debug)
- directories[p] = output_dir
- job = multiprocessing.Process(target=start_new_instance, args=(
- output_file, output_dir, artss_exe_path, 'artss_data_assimilation_serial'))
- job.start()
- counter += 1
- print(f'after start job {p} with name {job.name}, {job}')
- jobs[p] = job
- print(f'added job {p}')
- if counter >= no_cpus:
- # wait if you want to start more multi processes than cpus available
- j = jobs.pop(list(jobs.keys())[0])
- j.join()
- print('calc nabla parallel')
- new_para = {}
- for p in keys:
- if p in jobs.keys():
- job = jobs[p]
- print(f'wait for {job.name} [{job}]')
- job.join()
- diff_cur, _ = calc_diff(t_artss=t_artss, t_sensor=t_sensor,
- data_path=directories[p], fds_data=fds_data, devc_info=devc_info,
- file_da=file_da)
- print(p, diff_cur['T'])
- nabla[p] = (diff_cur['T'] - diff_orig['T']) / delta[p]
- log(f'nabla[{p}]={nabla[p]}', file_debug)
- if nabla[p] > 0:
- new_para[p] = cur[p] + delta[p]
-
- log(f'nabla without restrictions: {nabla}', file_debug)
- restrictions = {}
- if 'HRR' in keys:
- restrictions['HRR'] = 1
-
- if 'x0' in keys:
- if nabla['x0'] < 0:
- restrictions['x0'] = (domain.domain_param['X2'] - cur['x0']) / -nabla['x0']
- elif nabla['x0'] > 0:
- restrictions['x0'] = (cur['x0'] - domain.domain_param['X1']) / nabla['x0']
- else:
- restrictions['x0'] = 0
- if 'y0' in keys:
- if nabla['y0'] < 0:
- restrictions['y0'] = (domain.domain_param['Y2'] - cur['y0']) / -nabla['y0']
- elif nabla['y0'] > 0:
- restrictions['y0'] = (cur['y0'] - domain.domain_param['Y1']) / nabla['y0']
- else:
- restrictions['y0'] = 0
- if 'z0' in keys:
- if nabla['z0'] < 0:
- restrictions['z0'] = (domain.domain_param['Z2'] - cur['z0']) / -nabla['z0']
- elif nabla['z0'] > 0:
- restrictions['z0'] = (cur['z0'] - domain.domain_param['Z1']) / nabla['z0']
- else:
- restrictions['z0'] = 0
- log(f'restrictions: {restrictions}', file_debug)
- # search direction
- nabla = np.asarray([nabla[x] for x in keys])
- D = np.eye(len(keys)) # -> hesse matrix
- d = np.dot(-D, nabla)
-
- # calc alpha
- sigma = 0.01
-
- alpha_min = 1
- for k in keys:
- alpha_min = min(alpha_min, restrictions[k])
- alpha = 0.9 * alpha_min
- log(f'mins: 1.0, {restrictions.values()}', file_debug)
- log(f'alpha: {alpha}', file_debug)
- n = n_iterations
- while n > 0:
- x = np.asarray([cur[x] for x in keys])
- new_para = x + (alpha * d)
- new_para = {
- p: new_para[i]
- for i, p in enumerate(keys)
- }
- pprint(new_para)
-
- file_da.write(f'iterate:{n};')
- diff_cur, _ = do_rollback(client=client,
- t_sensor=t_sensor, t_artss=t_artss, t_revert=t_revert,
- new_para=new_para, heat_source=heat_source.values(),
- sub_file_name=n,
- artss_data_path=artss_data_path, artss=artss,
- fds_data=fds_data,
- devc_info=devc_info,
- file_da=file_da, file_debug=file_debug)
- log(f'conditional statement {diff_orig["T"] + np.dot(alpha * sigma * nabla, d)}', file_debug)
- if diff_cur['T'] < diff_orig['T'] + np.dot(alpha * sigma * nabla, d):
- log(f'found alpha: {alpha}', file_debug)
- log(f'found x_k+1: {new_para}', file_debug)
- log(f"al1: {np.dot(alpha * sigma * nabla, d)}", file_debug)
- log(f"als: {diff_cur['T']} < {diff_orig['T'] + np.dot(alpha * sigma * nabla, d)}", file_debug)
- cur = copy(new_para)
- break
-
- alpha = 0.7 * alpha
- n = n - 1
- log(f'ls: {diff_cur["T"]}', file_debug)
+ file_da.flush()
+ file_debug.flush()
def start(fds_data_path: str, fds_input_file_name: str, artss_data_path: str, artss_path: str, parallel=False):
@@ -643,22 +445,25 @@ def start(fds_data_path: str, fds_input_file_name: str, artss_data_path: str, ar
xml = XML(FieldReader.get_xml_file_name(artss_data_path), path=artss_data_path)
xml.read_xml()
- domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles, enable_computational_domain=xml.computational_domain)
+ domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles,
+ enable_computational_domain=xml.computational_domain)
domain.print_info()
domain.print_debug()
- devc_info_temperature, devc_info_thermocouple, fds_data = read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devices, fds_data = fds_utility.read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devc_info_temperature = devices['temperature']
print(devc_info_temperature)
- sensor_times = fds_data.index
+ sensor_times = fds_data.index[3:]
+ print('sensor times:', fds_data.index)
heat_source = xml.get_temperature_source()
- file_da = open('da_details.csv', 'w')
- file_debug = open('da_debug_details.dat', 'w')
+ file_da = open(os.path.join(artss_data_path, 'da_details.csv'), 'w')
+ file_debug = open(os.path.join(artss_data_path, 'da_debug_details.dat'), 'w')
delta = {
- 'HRR': float(heat_source['temperature_source']['HRR']) * 0.5,
+ 'HRR': float(heat_source['temperature_source']['HRR']) * 0.2,
'x0': domain.domain_param['dx'] * 10,
'y0': domain.domain_param['dy'] * 1,
'z0': domain.domain_param['dz'] * 1
@@ -672,6 +477,7 @@ def start(fds_data_path: str, fds_input_file_name: str, artss_data_path: str, ar
}
keys = ['HRR', 'x0']
+
if parallel:
continuous_gradient_parallel(client=client, file_da=file_da, file_debug=file_debug,
sensor_times=sensor_times,
@@ -696,25 +502,6 @@ def start(fds_data_path: str, fds_input_file_name: str, artss_data_path: str, ar
file_da.close()
-def wait_artss(t_sensor, artss_data_path, artss: ARTSS.XML):
- time_step = artss.get_output_time_step()
- if t_sensor % time_step == 0:
- t = t_sensor
- else:
- t = (t_sensor // time_step + 1) * time_step
- # time.sleep(30) # needed for artss to rewrite meta file
- t_cur = FieldReader.get_t_current(path=artss_data_path)
- pbar = tqdm(total=t)
- pbar.update(t_cur)
- while t_cur <= t:
- time.sleep(1)
- t_new = FieldReader.get_t_current(path=artss_data_path)
- pbar.update(t_new - t_cur)
- t_cur = t_new
-
- pbar.close()
-
-
def write_changes_xml(change: dict, source: list, file_name: str, file_da: TextIO, path='.') -> [str, str]:
source_type, temperature_source, random = data_assimilation.change_heat_source(*source,
changes={'source_type': {},
@@ -722,7 +509,6 @@ def write_changes_xml(change: dict, source: list, file_name: str, file_da: TextI
'random': {}})
write_da_data(file_da=file_da, parameters=temperature_source)
da = DAFile()
- da.create_config({'u': False, 'v': False, 'w': False, 'p': False, 'T': False, 'C': False})
da.create_temperature_source_changes(
source_type=source_type,
temperature_source=temperature_source,
@@ -734,138 +520,8 @@ def write_changes_xml(change: dict, source: list, file_name: str, file_da: TextI
return config_file_name, config_file_path
-def get_time_step_artss(t_sensor: float, artss_data_path: str, dt: float, time_back: float = 10) -> [float, float]:
- files = FieldReader.get_all_time_steps(path=artss_data_path)
- times = list(filter(lambda x: x >= t_sensor, files))
- t_revert = ([dt] + list(filter(lambda x: x <= max(0., t_sensor - time_back), files)))[-1]
- return times[0], t_revert
-
-
-def read_fds_data(input_path: str, input_file_name: str, artss: Domain) -> [dict, dict, pd.DataFrame]:
- full_name = os.path.join(input_path, input_file_name)
- str_devc = read_devc_from_input_file(full_name + '.fds')
- dict_devc = parse_devc(str_devc)
- chid_file_name = get_chid_from_input_file(full_name + '.fds')
-
- fds_data = read_devc_from_csv_file(os.path.join(input_path, chid_file_name + '_devc.csv'))
- # fds_data = read_fds_file(input_path, artss)
-
- devc_temperature = {}
- devc_thermocouple = {}
- for d in dict_devc:
- dict_devc[d]['type']: str = 'T'
-
- ijk = artss.get_ijk_from_xyz(dict_devc[d]['XYZ'][0], dict_devc[d]['XYZ'][2], dict_devc[d]['XYZ'][1])
- dict_devc[d]['index']: int = artss.get_index(*ijk)
- if d.startswith('Temperatur'):
- devc_temperature[d] = dict_devc[d]
- elif d.startswith('Thermocouple'):
- devc_thermocouple[d] = dict_devc[d]
-
- return devc_temperature, devc_thermocouple, fds_data
-
-
-def read_devc_from_csv_file(file_name: str) -> pd.DataFrame:
- df = pd.read_csv(file_name, skiprows=1)
- df.index = df.iloc[:, 0]
- df.drop('Time', axis=1, inplace=True)
- return df
-
-
-def get_chid_from_input_file(input_file: str) -> str:
- with open(input_file, 'r') as inp:
- for line in inp:
- if line.startswith('&HEAD'):
- line = line[len('&HEAD'):].strip()
- tmp_arr = line.split('=')
- index = tmp_arr.index('CHID')
- chid_str = tmp_arr[index + 1].strip()
- index_end = chid_str[1:].index("'")
- chid = chid_str[1:index_end + 1]
- break
- return chid
-
-
-def parse_devc(str_devc: list, start='&DEVC') -> dict:
- devc = {}
- for line in str_devc:
- line = line[len(start):-1]
- line.strip()
- tmp_arr = line.split(',')
- parameters = []
- for i in tmp_arr:
- if type(i) == str and '=' in i:
- parameters.append(i)
- else:
- parameters[-1] += f',{i}'
-
- param_dict = {}
- for p in parameters:
- tmp_arr = p.split('=')
- param_dict[tmp_arr[0].strip()] = tmp_arr[1].strip().strip("'").strip('"')
- id = param_dict.pop('ID')
- if 'XYZ' in param_dict.keys():
- param_dict['XYZ'] = [float(i) for i in param_dict['XYZ'].split(',')]
- devc[id] = param_dict
- return devc
-
-
-def read_devc_from_input_file(fds_file_name: str, start="&DEVC", end="/") -> list:
- """
- Reads FDS input out of a list of string and collects namelist definitions
- of a given type, e.g. "&DEVC ".
- Reads all definitions line by line, without processing.
-
- :param fds_file_name: List of string, of a single namelist definition.
- :param start: String, indicate the namelist, default "&DEVC ".
- :param end: String, indicate the end of the namelist, default "/".
- :return: list with lines
- inspired by Tristan Hehnen
- """
-
- namelist_lines = []
-
- collect_line = False
- fds_input_content = open(fds_file_name, 'r')
- for line in fds_input_content:
- line = line.strip()
- if collect_line is True:
- namelist_lines[-1] += f' {line}'
- elif line.startswith(start):
- namelist_lines.append(line)
- collect_line = True
- collect_line = collect_line and (end not in line)
- fds_input_content.close()
- return namelist_lines
-
-
-def read_fds_file(data_path: str, artss: Domain) -> pd.DataFrame:
- data = fdsreader.Simulation(data_path)
- devc = data.devices
- return_val: pd.DataFrame = pd.DataFrame([], index=devc['Time'].data)
- return_val.index.names = ['Time']
- return_val.columns.names = ['Position of Thermocouple']
- for key in devc:
- if key.startswith('Thermocouple'):
- if devc[key] is not None:
- print(key, devc[key].xyz, *devc[key].xyz, devc[key].xyz[1])
- # be careful, in FDS the third parameter indicates the height, but in ARTSS it is the second
- i, j, k = artss.get_ijk_from_xyz(devc[key].xyz[0], devc[key].xyz[2], devc[key].xyz[1])
- index = artss.get_index(i, j, k)
- print(f'index: {index} of {i}|{j}|{k}')
- if index in return_val:
- print(f'index {index} already exists, take average')
- tmp = (devc[key].data + return_val[index]) / 2.
- return_val[index] = tmp
- else:
- return_val[index] = devc[key].data
- else:
- print(f'{key} is None, skipped')
- return return_val
-
-
def comparison_sensor_simulation_data(devc_info: dict, sensor_data: pd.DataFrame, field_reader: FieldReader,
- t_sensor: float, file_da: TextIO) -> (Dict[str, float], float):
+ t_sensor: float, file_da: TextIO) -> Tuple[Dict[str, float], float]:
diff: dict = {'T': [], 'C': []}
fields_sim = field_reader.read_field_data()
min_pos_x: float = 0
@@ -878,10 +534,10 @@ def comparison_sensor_simulation_data(devc_info: dict, sensor_data: pd.DataFrame
value_sim = fields_sim[type_sensor][index_sensor]
difference = kelvin_to_celsius(value_sim) - value_sensor
- if abs(difference) > 0.5:
+ if abs(difference) > 2:
diff[type_sensor].append(difference)
file_da.write(
- f'sensor:{key};time_sensor:{t_sensor};time_artss:{field_reader.t};sensor_val:{value_sensor};artss_val:{value_sim};diff:{difference};considered:{abs(difference) > 0.5};position:{devc_info[key]["XYZ"]}\n')
+ f'sensor:{key};time_sensor:{t_sensor};time_artss:{field_reader.t};sensor_val:{value_sensor};artss_val:{value_sim};diff:{difference};considered:{abs(difference) > 2};position:{devc_info[key]["XYZ"]}\n')
file_da.flush()
if difference < min_sensor_val:
@@ -891,6 +547,8 @@ def comparison_sensor_simulation_data(devc_info: dict, sensor_data: pd.DataFrame
result = {}
for key in diff:
if len(diff[key]) == 0:
+ result[key] = 0
+ file_da.write(f'differences:{key}:{result[key]};no_of_sensors:{len(diff[key])}\n')
continue
result[key] = np.sqrt(sum(np.array(diff[key]) ** 2))
file_da.write(f'differences:{key}:{result[key]};no_of_sensors:{len(diff[key])}\n')
@@ -925,19 +583,17 @@ def calc_single_differences(devc_info: dict, sensor_data: pd.DataFrame, field_re
return diff
-def kelvin_to_celsius(kelvin):
- return kelvin - 273.5
-
-
def plot_differences(fds_data_path: str, fds_input_file_name: str, artss_data_path: str):
a_path = os.path.join(artss_data_path, '30')
xml = XML(FieldReader.get_xml_file_name(path=a_path), path=a_path)
xml.read_xml()
- domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles, enable_computational_domain=xml.computational_domain)
+ domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles,
+ enable_computational_domain=xml.computational_domain)
domain.print_info()
domain.print_debug()
- devc_info_temperature, devc_info_thermocouple, fds_data = read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devices, fds_data = fds_utility.read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devc_info_temperature = devices['temperature']
print(devc_info_temperature)
sensor_times = fds_data.index
@@ -946,9 +602,9 @@ def plot_differences(fds_data_path: str, fds_input_file_name: str, artss_data_pa
t_sensor = sensor_times[t]
fig, axs = plt.subplots(2, 1)
for pos in [20, 30, 40, 50, 60, 70, 80]:
- t_artss, t_revert = get_time_step_artss(t_sensor=t_sensor,
- artss_data_path=os.path.join(artss_data_path, str(pos)),
- dt=xml.get_dt())
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor=t_sensor,
+ artss_data_path=os.path.join(artss_data_path, str(pos)),
+ dt=xml.get_dt())
field_reader = FieldReader(t_artss, path=os.path.join(artss_data_path, str(pos)))
fields = field_reader.read_field_data()
diff = calc_single_differences(devc_info_temperature, fds_data, field_reader, t_artss, t_sensor)
@@ -1017,11 +673,13 @@ def plot_comparison_da(fds_data_path: str, fds_input_file_name: str, artss_data_
a_path = os.path.join(artss_data_path, 'with_da')
xml = XML(FieldReader.get_xml_file_name(path=a_path), path=a_path)
xml.read_xml()
- domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles, enable_computational_domain=xml.computational_domain)
+ domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles,
+ enable_computational_domain=xml.computational_domain)
domain.print_info()
domain.print_debug()
- devc_info_temperature, devc_info_thermocouple, fds_data = read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devices, fds_data = fds_utility.read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devc_info_temperature = devices['temperature']
print(devc_info_temperature)
sensor_times = fds_data.index[:31]
@@ -1062,11 +720,13 @@ def plot_sensor_data(fds_data_path: str, fds_input_file_name: str, artss_data_pa
a_path = os.path.join(artss_data_path, '30')
xml = XML(FieldReader.get_xml_file_name(path=a_path), path=a_path)
xml.read_xml()
- domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles, enable_computational_domain=xml.computational_domain)
+ domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles,
+ enable_computational_domain=xml.computational_domain)
domain.print_info()
domain.print_debug()
- devc_info_temperature, devc_info_thermocouple, fds_data = read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devices, fds_data = fds_utility.read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devc_info_temperature = devices['temperature']
print(devc_info_temperature)
sensor_times = fds_data.index[:13]
@@ -1143,11 +803,13 @@ def compare_distance(fds_data_path: str, fds_input_file_name: str, a_data_path:
xml = XML(FieldReader.get_xml_file_name(path=artss_data_path), path=artss_data_path)
xml.read_xml()
- domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles, enable_computational_domain=xml.computational_domain)
+ domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles,
+ enable_computational_domain=xml.computational_domain)
domain.print_info()
domain.print_debug()
- devc_info_temperature, devc_info_thermocouple, fds_data = read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devices, fds_data = fds_utility.read_fds_data(fds_data_path, fds_input_file_name, domain)
+ devc_info_temperature = devices['temperature']
log(f'devc info: {devc_info_temperature}', file_debug)
x_pos_orig = 52.5
@@ -1173,7 +835,7 @@ def compare_distance(fds_data_path: str, fds_input_file_name: str, a_data_path:
devc_info = devc_info_temperature
for t_sensor in sensor_times:
- t_artss, t_revert = get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(), time_back=6)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(), time_back=6)
log(f't_sensor: {t_sensor} t_artss: {t_artss} t_revert: {t_revert}', file_debug)
field_reader = FieldReader(t_artss, path=artss_data_path)
diff --git a/data_assimilation/main.py b/data_assimilation/main.py
index b7de4b52..7726db62 100644
--- a/data_assimilation/main.py
+++ b/data_assimilation/main.py
@@ -1,15 +1,22 @@
+import multiprocessing
import os
import time
+from pprint import pprint
+from typing import List, Dict, Tuple
import numpy as np
+import pandas as pd
from numpy import ndarray
+import ARTSS
import TCP_client
import data_assimilation
-import gradient_based_optimisation
+import downhill_simplex
+import fds_utility
from ARTSS import XML, Domain
from data_assimilation import FieldReader, DAFile
from obstacle import Obstacle
+from utility import log
def change_something(domain: Domain, field: list) -> list:
@@ -86,7 +93,7 @@ def main(dry_run=False):
'x0': float(source[1]['x0']) + 10 * index},
'random': {}})
da = DAFile()
- da.create_config({'u': False, 'v': False, 'w': False, 'p': False, 'T': True, 'C': False}, field_file_name)
+ da.create_field_changes({'T': True}, field_file_name)
da.create_temperature_source_changes(
source_type=source_type,
temperature_source=temperature_source,
@@ -109,14 +116,160 @@ def tmp():
obstacle.add_boundary(fields=['p'], patches=['back'], boundary_condition='neumann', value=10)
obstacle.add_boundary(fields=['p'], patches=['left'], boundary_condition='dirichlet', value=11)
da = DAFile()
- da.create_config({'u': False, 'v': False, 'w': False, 'p': False, 'T': False, 'C': False})
da.create_obstacle_changes([obstacle], True)
da.write_xml('change_obstacle.xml', pretty_print=True)
+def initial_guess_heat_source_position(time: float, fds_data: pd.DataFrame, temp_devices: Dict[str, any], threshold=22):
+ """
+ determine position of heat source based on the average of all sensors which values are higher than threshold
+ :param time:
+ :param fds_data:
+ :param temp_devices:
+ :param threshold:
+ :return:
+ """
+ positions = {'x': [], 'y': [], 'z': []}
+ for sensor in temp_devices:
+ if fds_data[sensor][time] > threshold:
+ x, y, z = temp_devices[sensor]['XYZ']
+ positions['x'].append(x)
+ positions['y'].append(z) # in ARTSS and FDS y and z are swapped
+ positions['z'].append(y)
+
+ x = sum(positions['x']) / len(positions['x'])
+ y = sum(positions['y']) / len(positions['y'])
+ z = sum(positions['z']) / len(positions['z'])
+ return x, y, z
+
+
+def start_simulation(fds_data_path: str, fds_input_file_name: str, artss_data_path: str, artss_input_file_name: str, artss_root_path: str, port=7777, time_difference: float = 0):
+ domain, xml = parse_artss_input_file(artss_data_path, artss_input_file_name=artss_input_file_name)
+
+ fds_data, devices, sensor_times = parse_fds_data(fds_data_path=fds_data_path, fds_input_file_name=fds_input_file_name, time_difference=time_difference, domain=domain)
+ pos_heat_source = initial_guess_heat_source_position(sensor_times[0], fds_data, devices['temperature'])
+
+ file_da = open(os.path.join(artss_data_path, 'da_details.csv'), 'w')
+ file_debug = open(os.path.join(artss_data_path, 'da_debug_details.dat'), 'w')
+
+ heat_source = xml.get_temperature_source()
+ delta = {
+ 'HRR': float(heat_source['temperature_source']['HRR']) * 0.05,
+ 'x0': 1.5,
+ 'y0': domain.domain_param['dy'] * 1,
+ 'z0': 0.5
+ }
+
+ cur = {
+ 'HRR': float(heat_source['temperature_source']['HRR']),
+ 'x0': pos_heat_source[0],
+ 'y0': pos_heat_source[1],
+ 'z0': pos_heat_source[2]
+ }
+
+ start_file_name = artss_input_file_name[:-4] + '_initial_guess.xml'
+ ARTSS.change_xml(change={'x0': cur['x0'], 'z0': cur['z0']},
+ input_file=os.path.join(artss_data_path, artss_input_file_name), output_file=start_file_name,
+ artss_data_path=artss_data_path, artss_root_path=artss_root_path,
+ file_debug=file_debug)
+ job = multiprocessing.Process(target=ARTSS.start_new_instance, args=(start_file_name, artss_data_path, os.path.join(artss_root_path, 'build', 'bin'), 'artss_data_assimilation_serial'))
+ job.start()
+ time.sleep(2) # wait for ARTSS to start
+ client = set_up_client(port=port)
+
+ keys = ['HRR', 'x0', 'z0']
+
+ log(f'sensor times: {sensor_times}', file_debug)
+ log(f'cur: {cur}', file_debug)
+ log(f'delta: {delta}', file_debug)
+ log(f'keys: {keys}', file_debug)
+
+ devc_info_temperature = devices['temperature']
+ print('devices:')
+ pprint(devc_info_temperature)
+
+ downhill_simplex.opt_scipy(client=client, file_da=file_da, file_debug=file_debug,
+ sensor_times=sensor_times,
+ devc_info=devc_info_temperature, fds_data=fds_data,
+ artss_data_path=artss_data_path,
+ domain=domain, heat_source=heat_source,
+ cur=cur, delta=delta, keys=keys, n_iterations=5, artss=xml)
+
+
+def set_up_client(ip_address: str = 'localhost', port: int = 7777) -> TCP_client:
+ """
+ set up TCP client and connect to ARTSS
+ :param ip_address: IP address for connecting to the running simulation (default: localhost)
+ :param port: port for connecting to the running simulation, can be found in input file (default: 7777)
+ :return: connected TCP client
+ """
+ client = TCP_client.TCPClient()
+ client.set_server_address(ip_address, port)
+ client.connect()
+ return client
+
+
+def parse_artss_input_file(artss_data_path: str, artss_input_file_name: str = None, quiet=False) -> Tuple[Domain, XML]:
+ artss_data_path = os.path.abspath(artss_data_path)
+ if artss_input_file_name is None:
+ artss_input_file_name = FieldReader.get_xml_file_name(artss_data_path)
+ xml = XML(artss_input_file_name, path=artss_data_path)
+ xml.read_xml()
+ domain = Domain(domain_param=xml.domain, obstacles=xml.obstacles,
+ enable_computational_domain=xml.computational_domain)
+ if not quiet:
+ domain.print_info()
+ domain.print_debug()
+ return domain, xml
+
+
+def parse_fds_data(fds_data_path: str, fds_input_file_name: str, domain: Domain, time_difference: float = 2, threshold=22, keys=['temperature']) -> [pd.DataFrame, Dict[str, Dict[str, any]],
+ List[float]]:
+ devices, fds_data = fds_utility.read_fds_data(fds_data_path, fds_input_file_name, domain)
+ starting_time = fds_utility.get_starting_time(fds_data, threshold, keys)
+ starting_time_index = list(fds_data.index).index(starting_time)
+ sensor_times = fds_data.index[starting_time_index:]
+ res = [sensor_times[0]]
+ for i in sensor_times[1:]:
+ if i > res[-1] + time_difference:
+ res.append(i)
+ return fds_data, devices, res
+
+
+def set_up(fds_data_path: str, fds_input_file_name: str, artss_data_path: str, port: int = 7777, time_difference: float = 0) \
+ -> [TCP_client, Domain, XML, List[float], Dict[str, Dict[str, any]], pd.DataFrame, float]:
+ """
+ set up necessary objects for data assimilation
+ :param fds_data_path: path to directory where the FDS and devc file are
+ :param fds_input_file_name: name of FDS file (without ending .fds)
+ :param artss_data_path: path to where the .vis directory of the simulation is
+ :param port: port defined in the XML used by the simulation
+ :param time_difference: restrict sensor times, two times have to be at least time_difference apart (in seconds)
+ :return: TCPClient: connecting to ARTSS
+ Domain: information about running simulation
+ XML: parsed ARTSS input file
+ sensor times
+ devices: data about devices from FDS input file combined with location data from ARTSS
+ fds_data: content of FDS device file
+ offset: time when the first (temperature) sensor in FDS rises about a certain threshold (22C)
+ """
+ client = set_up_client(port=port)
+ domain, xml = parse_artss_input_file(artss_data_path)
+ fds_data, devices, sensor_times = parse_fds_data(fds_data_path=fds_data_path, fds_input_file_name=fds_input_file_name, time_difference=time_difference, domain=domain)
+ print('sensor times:', fds_data.index)
+ return client, domain, xml, sensor_times, devices, fds_data
+
+
if __name__ == '__main__':
- gradient_based_optimisation.start(artss_data_path='example',
- fds_data_path='example/fds_data', fds_input_file_name='tunnel',
- artss_path=os.path.join(os.getcwd(), '..'),
- parallel=True)
+ start_simulation(fds_data_path='example/fds_data', fds_input_file_name='tunnel',
+ artss_data_path='../thesis', artss_input_file_name='tunnel.xml', artss_root_path='..',
+ time_difference=2)
+ # gradient_based_optimisation.start(artss_data_path='example',
+ # fds_data_path='example/fds_data', fds_input_file_name='tunnel',
+ # artss_path=os.path.join(os.getcwd(), '..'),
+ # parallel=True)
+ # downhill_simplex.start(artss_data_path='../thesis/run3',
+ # fds_data_path='example/fds_data', fds_input_file_name='tunnel',
+ # parallel=True,
+ # port=7777)
# main(dry_run=False)
diff --git a/data_assimilation/obstacle.py b/data_assimilation/obstacle.py
index aa6d798e..efcc6933 100644
--- a/data_assimilation/obstacle.py
+++ b/data_assimilation/obstacle.py
@@ -1,5 +1,5 @@
#!/usr/bin/env python3
-from typing import List, Dict
+from typing import List, Dict, Set
STATE: Dict[str, int] = {'XML': 0, 'unmodified': 1, 'modified': 2, 'new': 3, 'deleted': 4}
FIELD_TYPES: Dict[str, int] = {'u': 0, 'v': 1, 'w': 2, 'p': 3, 'T': 4, 'C': 5}
@@ -28,6 +28,8 @@ def __init__(self, name: str, state: str = 'unmodified'):
self.geometry: Dict[str, float] = {}
self.boundary: List[BoundaryData] = [BoundaryData(f) for f in FIELD_TYPES.keys()]
self.state = state
+ self.index: Dict[str, int] = {}
+ self.cells: Dict[str, Set[int]] = {}
def add_boundary(self, fields: List[str], patches: List[str], boundary_condition: str, value: float):
for f in fields:
@@ -36,7 +38,8 @@ def add_boundary(self, fields: List[str], patches: List[str], boundary_condition
def add_boundary_line(self, boundary: Dict[str, str]):
fields = boundary['field'].split(",")
patches = boundary['patch'].split(",")
- self.add_boundary(fields=fields, patches=patches, boundary_condition=boundary['type'], value=float(boundary['value']))
+ self.add_boundary(fields=fields, patches=patches, boundary_condition=boundary['type'],
+ value=float(boundary['value']))
def add_geometry(self, geometry: List[float]):
keys = ['ox1', 'ox2', 'oy1', 'oy2', 'oz1', 'oz2']
diff --git a/data_assimilation/obstacle_changer.py b/data_assimilation/obstacle_changer.py
index a8ddb77d..ee83a273 100644
--- a/data_assimilation/obstacle_changer.py
+++ b/data_assimilation/obstacle_changer.py
@@ -1,16 +1,14 @@
#!/usr/bin/env python3
import os
-from tqdm import tqdm
-import time
-from typing import List, Dict
+from typing import List, Dict, Tuple
import numpy as np
import TCP_client
from ARTSS import XML, Domain
from data_assimilation import FieldReader, create_message, DAFile
-from gradient_based_optimisation import get_time_step_artss
from obstacle import Obstacle
+from utility import wait_artss
def create_obstacles() -> List[Obstacle]:
@@ -67,11 +65,10 @@ def obstacle_wonder(artss_data_path: str):
for t_sensor in sensor_times[:8]:
wait_artss(t_sensor, artss_data_path)
- t_artss, t_revert = get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
- time_back=time_back * xml.get_dt())
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
da = DAFile()
- da.create_config({'u': False, 'v': False, 'w': False, 'p': False, 'T': False, 'C': False})
da.create_obstacle_changes([obstacles[counter]], True)
counter = (counter + 1) % 4
config_file_name = f'change_obstacle_{int(t_sensor * 10)}.xml'
@@ -83,11 +80,10 @@ def obstacle_wonder(artss_data_path: str):
for t_sensor in sensor_times[8:]:
wait_artss(t_sensor, artss_data_path)
- t_artss, t_revert = get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
- time_back=time_back * xml.get_dt())
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
da = DAFile()
- da.create_config({'u': False, 'v': False, 'w': False, 'p': False, 'T': False, 'C': False})
da.create_obstacle_changes([obstacles[counter], obstacles[counter + 2]], True)
counter = (counter + 1) % 2
config_file_name = f'change_obstacle_{int(t_sensor * 10)}.xml'
@@ -96,8 +92,165 @@ def obstacle_wonder(artss_data_path: str):
client.send_message(create_message(t_revert, config_file_path))
+def full_corridor_doors(artss_data_path: str):
+ client = TCP_client.TCPClient()
+ client.connect()
+
+ xml = XML(FieldReader.get_xml_file_name(artss_data_path), path=artss_data_path)
+ xml.read_xml()
+
+ obstacles = xml.obstacles
+ door1, neighbours1 = create_door(obstacles, door_identifier='Room1', name='door1')
+ door2, neighbours2 = create_door(obstacles, door_identifier='Room2', name='door2')
+ door3, neighbours3 = create_door(obstacles, door_identifier='Room3', name='door3')
+ door4, neighbours4 = create_door(obstacles, door_identifier='Room4', name='door4')
+ door5, neighbours5 = create_door(obstacles, door_identifier='Room5', name='door5')
+ door6, neighbours6 = create_door(obstacles, door_identifier='Room6', name='door6')
+
+ doors = [door1, door2, door3, door4, door5, door6]
+ neighbours = [neighbours1, neighbours2, neighbours3, neighbours4, neighbours5, neighbours6]
+
+ time_back = 1
+ domain = Domain(xml.domain, xml.computational_domain, obstacles)
+ for counter in range(len(doors)):
+ print('door', counter)
+ # close door
+ t_sensor = 11 * xml.get_dt() + time_back * xml.get_dt() + counter * xml.get_dt() * 4
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
+ domain.add_obstacle(doors[counter])
+ domain.print_info()
+ da = DAFile()
+ reader = FieldReader(t_artss, path=artss_data_path)
+ field_changes, fields = set_value(fields=reader.read_field_data(), domain=domain,
+ obstacle_names=[doors[counter].name],
+ neighbouring_obstacle_patches={doors[counter].name: neighbours[counter]})
+ field_key_changes: List[str] = merge_field_keys(field_changes, field_key_changes=[])
+ field_file_name = f'fields_{t_artss:.5e}.hdf'
+ FieldReader.write_field_data_keys(file_name=field_file_name, data=fields, field_keys=field_key_changes,
+ path=artss_data_path)
+ da.create_field_changes(field_changes, field_file_name=field_file_name)
+ da.write_obstacle_changes(domain.get_obstacles())
+ config_file_name = f'full_corridor_obstacle_{int(t_sensor * 10)}.xml'
+ config_file_path = os.path.join(artss_data_path, config_file_name)
+ da.write_xml(config_file_path, pretty_print=True)
+ client.send_message(create_message(t_revert, config_file_name))
+
+ # open door
+ t_sensor = 13 * xml.get_dt() + time_back * xml.get_dt() + counter * xml.get_dt() * 4
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
+ door = domain.remove_obstacle(doors[counter].name)
+ domain.print_info()
+ da = DAFile()
+ reader = FieldReader(t_artss, path=artss_data_path)
+ field_changes, fields = set_ambient_temperature(domain=domain, obstacles=[door], value=299.14,
+ fields=reader.read_field_data())
+ field_key_changes = merge_field_keys(field_changes, field_key_changes=[])
+ field_file_name = f'fields_{t_artss:.5e}.hdf'
+ FieldReader.write_field_data_keys(file_name=field_file_name, data=fields, field_keys=field_key_changes,
+ path=artss_data_path)
+ da.create_field_changes(field_changes, field_file_name=field_file_name)
+ da.write_obstacle_changes(domain.get_obstacles() + [door])
+ config_file_name = f'full_corridor_obstacle_{int(t_sensor * 10)}.xml'
+ config_file_path = os.path.join(artss_data_path, config_file_name)
+ da.write_xml(config_file_path, pretty_print=True)
+ client.send_message(create_message(t_revert, config_file_name))
+
+
+def merge_field_keys(field_changes: Dict[str, bool], field_key_changes: List[str]):
+ for f in field_changes:
+ if field_changes[f] and f not in field_key_changes:
+ field_key_changes.append(f)
+ return field_key_changes
+
+
+def full_corridor_rooms(artss_data_path: str):
+ client = TCP_client.TCPClient()
+ client.connect()
+
+ xml = XML(FieldReader.get_xml_file_name(artss_data_path), path=artss_data_path)
+ xml.read_xml()
+
+ obstacles = xml.obstacles
+
+ domain = Domain(xml.domain, xml.computational_domain, obstacles)
+ rooms = [replace_room2(domain)]
+
+ time_back = 1
+ for counter in range(len(rooms)):
+ print('room', counter + 2)
+ # close door
+ t_sensor = 11 * xml.get_dt() + time_back * xml.get_dt() + counter * xml.get_dt() * 4
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
+ removed = []
+ for o in rooms[counter][0]:
+ domain.add_obstacle(o)
+ for o_name in rooms[counter][1]:
+ removed.append(domain.remove_obstacle(o_name))
+
+ o_names = [x.name for x in rooms[counter][0]]
+ domain.print_info()
+ reader = FieldReader(t_artss, path=artss_data_path)
+ field_changes, fields = set_value(fields=reader.read_field_data(), domain=domain,
+ obstacle_names=o_names,
+ neighbouring_obstacle_patches={})
+ field_key_changes: List[str] = merge_field_keys(field_changes, field_key_changes=[])
+ field_file_name = f'fields_{t_artss:.5e}.hdf'
+ FieldReader.write_field_data_keys(file_name=field_file_name, data=fields, field_keys=field_key_changes,
+ path=artss_data_path)
+ da = DAFile()
+ da.create_field_changes(field_changes, field_file_name=field_file_name)
+ da.write_obstacle_changes(domain.get_obstacles() + removed)
+ config_file_name = f'full_corridor_obstacle_{int(t_sensor * 10)}.xml'
+ config_file_path = os.path.join(artss_data_path, config_file_name)
+ da.write_xml(config_file_path, pretty_print=True)
+ client.send_message(create_message(t_revert, config_file_name))
+
+ # open door
+ t_sensor = 13 * xml.get_dt() + time_back * xml.get_dt() + counter * xml.get_dt() * 4
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
+ for o in removed:
+ domain.add_obstacle(o)
+ o_names = [x.name for x in removed]
+ removed = []
+ for o in rooms[counter][0]:
+ removed.append(domain.remove_obstacle(o.name))
+ domain.print_info()
+ da = DAFile()
+ reader = FieldReader(t_artss, path=artss_data_path)
+ field_changes, fields = set_ambient_temperature(fields=reader.read_field_data(), domain=domain,
+ obstacles=removed,
+ value=293.15)
+ field_key_changes = merge_field_keys(field_changes, field_key_changes=[])
+ field_changes, fields = set_value(fields=fields, domain=domain,
+ obstacle_names=o_names,
+ neighbouring_obstacle_patches={})
+ field_key_changes = merge_field_keys(field_changes, field_key_changes=field_key_changes)
+
+ field_file_name = f'fields_{t_artss:.5e}.hdf'
+ FieldReader.write_field_data_keys(file_name=field_file_name, data=fields, field_keys=field_key_changes,
+ path=artss_data_path)
+ da.create_field_changes(field_changes, field_file_name=field_file_name)
+ da.write_obstacle_changes(domain.get_obstacles() + removed)
+ config_file_name = f'change_obstacle_{int(t_sensor * 10)}.xml'
+ config_file_path = os.path.join(artss_data_path, config_file_name)
+ da.write_xml(config_file_path, pretty_print=True)
+ client.send_message(create_message(t_revert, config_file_name))
+
+
def steckler_door(artss_data_path: str):
- artss_data_path = os.path.abspath(artss_data_path)
+ """
+ experimental setup to close and open the door
+ :param artss_data_path: path to where ARTSS was started
+ """
+ # artss_data_path = os.path.abspath(artss_data_path)
client = TCP_client.TCPClient()
client.connect()
@@ -110,48 +263,135 @@ def steckler_door(artss_data_path: str):
door, neighbours = create_door(obstacles)
time_back = 1
- t_sensor = 81 * xml.get_dt() + time_back * xml.get_dt()
+ # close door
+ t_sensor = 81 * xml.get_dt() + time_back * xml.get_dt()
wait_artss(t_sensor, artss_data_path)
-
- t_artss, t_revert = get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
- time_back=time_back * xml.get_dt())
-
- # da = DAFile()
- # da.create_config({'u': False, 'v': False, 'w': False, 'p': False, 'T': False, 'C': False})
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
domain = Domain(xml.domain, xml.computational_domain, obstacles)
domain.add_obstacle(door)
- domain.print_info()
- print()
domain.print_debug()
- da = set_zero(t_artss, t_revert, domain=domain, obstacle_name=door.name, neighbouring_obstacle_patches=neighbours)
- da.write_obstacle_changes(obstacles + [door], True)
- # da.remove_obstacle(obstacles, door)
+ da = DAFile()
+ reader = FieldReader(t_artss, path=artss_data_path)
+ field_changes, fields = set_value(domain=domain, obstacle_names=[door.name],
+ neighbouring_obstacle_patches={door.name: neighbours},
+ fields=reader.read_field_data())
+ field_key_changes: List[str] = []
+ for f in field_changes:
+ if field_changes[f] and f not in field_key_changes:
+ field_key_changes.append(f)
+ field_file_name = f'fields_{t_artss:.5e}.hdf'
+ FieldReader.write_field_data_keys(file_name=field_file_name, data=fields, field_keys=field_key_changes,
+ path=artss_data_path)
+ da.create_field_changes(field_changes, field_file_name=field_file_name)
+ da.write_obstacle_changes(domain.get_obstacles())
config_file_name = f'change_obstacle_{int(t_sensor * 10)}.xml'
config_file_path = os.path.join(artss_data_path, config_file_name)
da.write_xml(config_file_path, pretty_print=True)
- client.send_message(create_message(t_revert, config_file_path))
+ client.send_message(create_message(t_revert, config_file_name))
+ # open door
+ t_sensor = 101 * xml.get_dt() + time_back * xml.get_dt() # 101
+ wait_artss(t_sensor, artss_data_path)
+ t_artss, t_revert = FieldReader.get_time_step_artss(t_sensor, artss_data_path, dt=xml.get_dt(),
+ time_back=time_back * xml.get_dt())
+ door = domain.remove_obstacle(door.name)
+ domain.print_debug()
+ da = DAFile()
+ [field_changes, fields] = set_ambient_temperature(domain=domain, obstacles=[door], value=299.14,
+ fields=reader.read_field_data())
+ # [field_changes, field_file_name] = set_gradient_x(t_artss=t_artss, domain=domain, obstacle=door,
+ # artss_data_path=artss_data_path)
+ field_key_changes: List[str] = []
+ for f in field_changes:
+ if field_changes[f] and f not in field_key_changes:
+ field_key_changes.append(f)
+ field_file_name = f'fields_{t_artss:.5e}.hdf'
+ FieldReader.write_field_data_keys(file_name=field_file_name, data=fields, field_keys=field_key_changes,
+ path=artss_data_path)
+ da.create_field_changes(field_changes, field_file_name=field_file_name)
+ da.write_obstacle_changes(domain.get_obstacles() + [door])
+ config_file_name = f'change_obstacle_{int(t_sensor * 10)}.xml'
+ config_file_path = os.path.join(artss_data_path, config_file_name)
+ da.write_xml(config_file_path, pretty_print=True)
+ client.send_message(create_message(t_revert, config_file_name))
+
+
+def set_ambient_temperature(fields: Dict[str, np.ndarray], obstacles: List[Obstacle], domain: Domain, value: float) -> [
+ Dict[str, bool], str]:
+ for obstacle in obstacles:
+ for j in range(obstacle.index['y1'], obstacle.index['y2'] + 1):
+ for k in range(obstacle.index['z1'], obstacle.index['z2'] + 1):
+ for i in range(obstacle.index['x1'], obstacle.index['x2'] + 1):
+ index = domain.get_index(i, j, k)
+ fields['T'][index] = value
+
+ return {'T': True}, fields
-def set_zero(t_artss: float, t_revert: float, obstacle_name: str, domain: Domain, neighbouring_obstacle_patches: Dict[str, List[str]]) -> DAFile:
- reader = FieldReader(t_artss, path='example')
+
+def set_gradient_x(t_artss: float, obstacle: Obstacle, domain: Domain, artss_data_path: str) -> [Dict[str, bool], str]:
+ """
+ set gradient in x-direction. interpolate cell values from left neighbour cell to right neighbour cell
+ :param t_artss: time of field to be changed
+ :param obstacle: obstacle which will be removed
+ :param domain: domain of ARTSS
+ :param artss_data_path: path to where ARTSS was executed
+ :return: A dictionary containing the changed fields and the name of changed field file
+ """
+ reader = FieldReader(t_artss, path=artss_data_path)
field_T = reader.read_field_data()['T']
- domain.set_value_of_obstacle_cells(value=-1, field=field_T, obstacle_name=obstacle_name)
- for o_name in neighbouring_obstacle_patches:
- for patch in neighbouring_obstacle_patches[o_name]:
- domain.set_value_of_obstacle_patch(value=-1, field=field_T, obstacle_name=o_name, patch=patch)
- field_file_name = f'temperature_{t_revert:.5e}'
+ field_file_name = f'temperature_{t_artss:.5e}'
+
+ for j in range(obstacle.index['y1'], obstacle.index['y2'] + 1):
+ for k in range(obstacle.index['z1'], obstacle.index['z2'] + 1):
+ start_val = field_T[domain.get_index(obstacle.index['x1'] - 1, j, k)]
+ end_val = field_T[domain.get_index(obstacle.index['x2'] + 1, j, k)]
+ delta = (end_val - start_val) / (obstacle.index['x2'] - obstacle.index['x1'])
+ counter = 1
+ print(
+ f"start val: {field_T[domain.get_index(obstacle.index['x1'] - 1, j, k)]}, index: {domain.get_index(obstacle.index['x1'] - 1, j, k)}")
+ print(
+ f"end val: {field_T[domain.get_index(obstacle.index['x2'] + 1, j, k)]}, index: {domain.get_index(obstacle.index['x2'] + 1, j, k)}")
+ print(f"delta: {(end_val - start_val) / (obstacle.index['x2'] - obstacle.index['x1'])}")
+ for i in range(obstacle.index['x1'], obstacle.index['x2'] + 1):
+ index = domain.get_index(i, j, k)
+ field_T[index] = start_val + delta * counter
+ counter += 1
FieldReader.write_field_data_keys(file_name=field_file_name,
data={'T': field_T},
- field_keys=['T'])
- da = DAFile()
- da.create_config({'u': False, 'v': False, 'w': False, 'p': False, 'T': True, 'C': False},
- field_file_name=os.path.abspath(field_file_name))
- return da
+ field_keys=['T'],
+ path=artss_data_path)
+ return {'T': True}, field_file_name
+
+
+def set_value(fields: Dict[str, np.ndarray], obstacle_names: List[str], domain: Domain,
+ neighbouring_obstacle_patches: Dict[str, Dict[str, List[str]]], value=0) -> [Dict[str, bool],
+ Dict[str, np.ndarray]]:
+ """
+ set all cells given obstacle to zero in order to see the obstacle in the simulation
+ :param fields: fields to be changed
+ :param obstacle_names: name of obstacles which cell values should be changed
+ :param domain: domain of ARTSS
+ :param neighbouring_obstacle_patches: specify additional cells which should be changed.
+ E.g. the patch of a neighbour obstacle
+ :return:
+ """
+ for obstacle_name in obstacle_names:
+ for field in fields:
+ domain.set_value_of_obstacle_cells(value=value, field=fields[field], obstacle_name=obstacle_name)
+ if obstacle_name in neighbouring_obstacle_patches:
+ for o_name in neighbouring_obstacle_patches[obstacle_name]:
+ for patch in neighbouring_obstacle_patches[obstacle_name][o_name]:
+ domain.set_value_of_obstacle_patch(value=value, field=fields[field], obstacle_name=o_name,
+ patch=patch)
+
+ return dict(zip(fields.keys(), [True] * len(fields))), fields
-def create_door(obstacles: List[Obstacle], door_identifier: str = None, ground_coordinate: float = 0) -> [Obstacle, Dict[str, List[int]]]:
+def create_door(obstacles: List[Obstacle], door_identifier: str = None, ground_coordinate: float = 0, name='door') \
+ -> [Obstacle, Dict[str, List[int]]]:
"""
door measurements are calculated based on the obstacle names "left from door", "right from door", "above door" +
door identifier. E.g. door identifier = "Room1" results in looking for obstacle names which include
@@ -168,10 +408,8 @@ def create_door(obstacles: List[Obstacle], door_identifier: str = None, ground_c
neighbours[obstacle.name] = ['bottom']
elif names[1] in obstacle.name:
obstacle_right = obstacle
- neighbours[obstacle.name] = ['left']
elif names[0] in obstacle.name:
obstacle_left = obstacle
- neighbours[obstacle.name] = ['right']
if obstacle_left.geometry['ox2'] < obstacle_right.geometry['ox1']:
# door faces z-direction
@@ -181,6 +419,8 @@ def create_door(obstacles: List[Obstacle], door_identifier: str = None, ground_c
door_coordinates['oy2'] = obstacle_above.geometry['oy1']
door_coordinates['oz1'] = max(obstacle_left.geometry['oz1'], obstacle_right.geometry['oz1'])
door_coordinates['oz2'] = min(obstacle_left.geometry['oz2'], obstacle_right.geometry['oz2'])
+ neighbours[obstacle_left.name] = ['right']
+ neighbours[obstacle_right.name] = ['left']
else:
# door faces x-direction
door_coordinates['ox1'] = max(obstacle_left.geometry['ox1'], obstacle_right.geometry['ox1'])
@@ -189,26 +429,37 @@ def create_door(obstacles: List[Obstacle], door_identifier: str = None, ground_c
door_coordinates['oy2'] = obstacle_above.geometry['oy1']
door_coordinates['oz1'] = obstacle_left.geometry['oz2']
door_coordinates['oz2'] = obstacle_right.geometry['oz1']
+ neighbours[obstacle_left.name] = ['back']
+ neighbours[obstacle_right.name] = ['front']
- door = Obstacle(name="door", state='new')
+ door = Obstacle(name=name, state='new')
door.geometry = door_coordinates
door.boundary = obstacle_above.boundary
return door, neighbours
-def wait_artss(t_sensor: float, artss_data_path: str):
- t_cur = FieldReader.get_t_current(path=artss_data_path)
- pbar = tqdm(total=t_sensor)
- pbar.update(t_cur)
- while t_cur <= t_sensor:
- time.sleep(1)
- t_new = FieldReader.get_t_current(path=artss_data_path)
- pbar.update(t_new - t_cur)
- t_cur = t_new
+def replace_room2(domain: Domain) -> Tuple[List[Obstacle], List[str]]:
+ created: List[Obstacle] = []
+ room_blocker = Obstacle(name='room2', state='new')
+ room_blocker.add_geometry([domain.obstacles['wall between 1 and 2'].geometry['ox1'],
+ domain.obstacles['wall between 2 and 3'].geometry['ox2'],
+ domain.obstacles['wall between 1 and 2'].geometry['oy1'],
+ domain.obstacles['wall between 1 and 2'].geometry['oy2'],
+ domain.obstacles['wall between 1 and 2'].geometry['oz1'],
+ domain.obstacles['wall between 1 and 2'].geometry['oz2'],
+ ])
+ created.append(room_blocker)
- pbar.close()
+ deleted: List[str] = ['wall between 1 and 2',
+ 'wall between 2 and 3',
+ 'room 2 right from doorRoom2',
+ 'room 2 above doorRoom2',
+ 'room 2 left from doorRoom2']
+ return created, deleted
if __name__ == '__main__':
# obstacle_wonder(artss_data_path='example')
- steckler_door(artss_data_path='example')
+ # steckler_door(artss_data_path='example')
+ # full_corridor_doors(artss_data_path='full_corridor')
+ full_corridor_rooms(artss_data_path='full_corridor')
diff --git a/data_assimilation/requirements.txt b/data_assimilation/requirements.txt
index d066c09b..8e4a0832 100644
--- a/data_assimilation/requirements.txt
+++ b/data_assimilation/requirements.txt
@@ -4,4 +4,5 @@ h5py
fdsreader
pandas
numpy
-matplotlib
\ No newline at end of file
+matplotlib
+scipy
\ No newline at end of file
diff --git a/data_assimilation/utility.py b/data_assimilation/utility.py
new file mode 100644
index 00000000..80db1a13
--- /dev/null
+++ b/data_assimilation/utility.py
@@ -0,0 +1,39 @@
+import time
+from typing import TextIO
+
+from tqdm import tqdm
+
+from data_assimilation import FieldReader
+
+
+def wait_artss(wait_time: float, artss_data_path: str):
+ t_cur = FieldReader.get_t_current(path=artss_data_path)
+ pbar = tqdm(total=wait_time)
+ pbar.update(t_cur)
+ while t_cur <= wait_time:
+ time.sleep(1)
+ t_new = FieldReader.get_t_current(path=artss_data_path)
+ pbar.update(t_new - t_cur)
+ t_cur = t_new
+
+ pbar.close()
+
+
+def write_da_data(file_da: TextIO, parameters: dict):
+ """
+ write sensor data and ARTSS data to file
+ :param file_da: name of file to write to
+ :param parameters:
+ """
+ for key in parameters:
+ file_da.write(f';{key}:{parameters[key]}')
+ file_da.write('\n')
+
+
+def log(message: str, file_debug: TextIO):
+ print(message)
+ file_debug.write(f'{message}\n')
+
+
+def kelvin_to_celsius(kelvin):
+ return kelvin - 273.5
diff --git a/gtests/randomField/UniformRandom.cpp b/gtests/randomField/UniformRandom.cpp
index 9da3dd05..b856cc12 100644
--- a/gtests/randomField/UniformRandom.cpp
+++ b/gtests/randomField/UniformRandom.cpp
@@ -2,6 +2,9 @@
#include
#include "src/randomField/UniformRandom.h"
+#include "src/source/GaussFunction.h"
+#include "src/domain/DomainData.h"
+#include "src/domain/DomainController.h"
#define EPS 10e-5
@@ -216,7 +219,7 @@ TEST_F(UniformRandomFieldTest, range_100_at_least_one) {
real no = -range + j * step_size;
bool found = false;
for (size_t i = 0; i < size; ++i) {
- if (abs(a[i] == no) < EPS) {
+ if (abs(a[i] - no) < EPS) {
found = true;
break;
}
@@ -410,3 +413,33 @@ TEST_F(UniformRandomFieldTest, relative) {
EXPECT_DOUBLE_EQ(a[i], result[i]);
}
}
+
+TEST_F(UniformRandomFieldTest, gauss) {
+ Settings::domain_parameters domain_params{ };
+ domain_params.start_coords_PD.set_coordinate(0, 0, 0);
+ domain_params.start_coords_CD.copy(domain_params.start_coords_PD);
+ domain_params.end_coords_PD.set_coordinate(5, 10, 15);
+ domain_params.end_coords_CD.copy(domain_params.end_coords_PD);
+ domain_params.enable_computational_domain = false;
+ domain_params.number_of_inner_cells.set_coordinate(5, 20, 10);
+ DomainData::init({ }, domain_params, 0);
+ Settings::Settings settings;
+ settings.domain_parameters = domain_params;
+ DomainController::init(settings);
+
+ Coordinate position(10, 9, 8);
+ Coordinate dimension(1, 2, 3);
+ Settings::solver::sources::gauss settings_gauss = {100, 1.5, position, dimension, 0.5};
+
+ Field a1(FieldType::UNKNOWN_FIELD, 100);
+ auto gauss1 = GaussFunction(settings_gauss);
+ gauss1.update_source(a1, 0);
+
+ Field a2(FieldType::UNKNOWN_FIELD, 100);
+ auto gauss2 = GaussFunction(settings_gauss);
+ gauss2.update_source(a2, 0);
+
+ for (size_t i = 0; i < DomainData::getInstance()->get_size(); i++){
+ EXPECT_DOUBLE_EQ(a1[i], a2[i]);
+ }
+}
diff --git a/gtests/utility/Settings.cpp b/gtests/utility/Settings.cpp
index 85acfca8..dd8db1d6 100644
--- a/gtests/utility/Settings.cpp
+++ b/gtests/utility/Settings.cpp
@@ -423,8 +423,8 @@ TEST(SettingsTest, requiredObstaclesParameters) {
)";
tinyxml2::XMLDocument doc;
doc.Parse(xml.c_str());
- Settings::obstacles_parameters obstacles_parameters = Settings::parse_obstacles_parameters(doc.RootElement());
- EXPECT_FALSE(obstacles_parameters.enabled);
+ Settings::obstacles_parameters obstacle_parameters = Settings::parse_obstacles_parameters(doc.RootElement());
+ EXPECT_FALSE(obstacle_parameters.enabled);
}
TEST(SettingsTest, obstacles) {
@@ -445,12 +445,13 @@ TEST(SettingsTest, obstacles) {
)";
tinyxml2::XMLDocument doc;
doc.Parse(xml.c_str());
- Settings::obstacles_parameters obstacles_parameters = Settings::parse_obstacles_parameters(doc.RootElement());
- EXPECT_TRUE(obstacles_parameters.enabled);
- EXPECT_EQ(obstacles_parameters.obstacles.size(), 2);
+ Settings::obstacles_parameters obstacle_parameters = Settings::parse_obstacles_parameters(doc.RootElement());
+ EXPECT_TRUE(obstacle_parameters.enabled);
+ EXPECT_EQ(obstacle_parameters.obstacles.size(), 2);
- const auto &obstacle1 = obstacles_parameters.obstacles[0];
+ const auto &obstacle1 = obstacle_parameters.obstacles[0];
EXPECT_EQ(obstacle1.name, "ceiling");
+ EXPECT_EQ(obstacle1.state, State::XML);
EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::X], 0);
EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::Y], 2.18);
EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::Z], -1.4);
@@ -474,7 +475,7 @@ TEST(SettingsTest, obstacles) {
o1_b1.field_type.end(),
FieldType::W),
o1_b1.field_type.end());
- for (Patch patch : {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
EXPECT_NE(std::find(o1_b1.patch.begin(),
o1_b1.patch.end(),
patch),
@@ -491,15 +492,16 @@ TEST(SettingsTest, obstacles) {
o1_b2.field_type.end(),
FieldType::T),
o1_b2.field_type.end());
- for (Patch patch : {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
EXPECT_NE(std::find(o1_b2.patch.begin(),
o1_b2.patch.end(),
patch),
o1_b2.patch.end());
}
- const auto &obstacle2 = obstacles_parameters.obstacles[1];
+ const auto &obstacle2 = obstacle_parameters.obstacles[1];
EXPECT_EQ(obstacle2.name, "right wall left from door");
+ EXPECT_EQ(obstacle2.state, State::XML);
EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::X], 2.8);
EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::Y], 0);
EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::Z], -1.4);
@@ -522,7 +524,7 @@ TEST(SettingsTest, obstacles) {
o2_b1.field_type.end(),
FieldType::P),
o2_b1.field_type.end());
- for (Patch patch : {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
EXPECT_NE(std::find(o2_b1.patch.begin(),
o2_b1.patch.end(),
patch),
@@ -539,13 +541,80 @@ TEST(SettingsTest, obstacles) {
o2_b2.field_type.end(),
FieldType::T),
o2_b2.field_type.end());
- for (Patch patch : {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
EXPECT_NE(std::find(o2_b2.patch.begin(),
o2_b2.patch.end(),
patch),
o2_b2.patch.end());
}
}
+TEST(SettingsTest, obstacles2) {
+ std::string xml = R"(
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+)";
+ tinyxml2::XMLDocument doc;
+ doc.Parse(xml.c_str());
+ Settings::obstacles_parameters obstacle_parameters = Settings::parse_obstacles_parameters(doc.RootElement());
+ EXPECT_TRUE(obstacle_parameters.enabled);
+ EXPECT_EQ(obstacle_parameters.obstacles.size(), 7);
+
+ const auto &obstacle1 = obstacle_parameters.obstacles[1];
+ EXPECT_EQ(obstacle1.name, "ceiling");
+ EXPECT_EQ(obstacle1.state, State::XML);
+ EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::X], -1.6625);
+ EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::Y], 2.13);
+ EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::Z], -1.6625);
+ EXPECT_DOUBLE_EQ(obstacle1.end_coords[CoordinateAxis::X], 1.6625);
+ EXPECT_DOUBLE_EQ(obstacle1.end_coords[CoordinateAxis::Y], 2.3296875);
+ EXPECT_DOUBLE_EQ(obstacle1.end_coords[CoordinateAxis::Z], 1.6625);
+
+ const auto &obstacle2 = obstacle_parameters.obstacles[4];
+ EXPECT_EQ(obstacle2.name, "right wall left from door");
+ EXPECT_EQ(obstacle2.state, State::XML);
+ EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::X], 1.4);
+ EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::Y], 0);
+ EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::Z], -1.4);
+ EXPECT_DOUBLE_EQ(obstacle2.end_coords[CoordinateAxis::X], 1.6625);
+ EXPECT_DOUBLE_EQ(obstacle2.end_coords[CoordinateAxis::Y], 2.13);
+ EXPECT_DOUBLE_EQ(obstacle2.end_coords[CoordinateAxis::Z], -0.4375);
+}
TEST(SettingsTest, requiredSurfacesParameters) {
std::string xml = R"(
@@ -689,7 +758,7 @@ TEST(SettingsTest, boundaries) {
patch),
boundaries_parameters.boundaries[1].patch.end());
}
- for (Patch patch : {Patch::LEFT, Patch::RIGHT, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::TOP, Patch::FRONT, Patch::BACK}) {
EXPECT_NE(std::find(boundaries_parameters.boundaries[2].patch.begin(),
boundaries_parameters.boundaries[2].patch.end(),
patch),
@@ -1766,13 +1835,12 @@ TEST(SettingsTest, assimilation) {
TEST(SettingsTest, assimilation2) {
std::string xml = R"(
-
+
)";
tinyxml2::XMLDocument doc;
doc.Parse(xml.c_str());
Settings::data_assimilation_parameters da = Settings::parse_assimilation_parameters(doc.RootElement());
EXPECT_TRUE(da.enabled);
- EXPECT_EQ("TemperatureSource", da.class_name);
EXPECT_FALSE(da.load_data);
EXPECT_EQ(da.output_dir, ".vis");
}
@@ -1780,29 +1848,7 @@ TEST(SettingsTest, assimilation2) {
TEST(SettingsTest, assimilation3) {
std::string xml = R"(
-
- 1
-
- 10
-
-)";
- tinyxml2::XMLDocument doc;
- doc.Parse(xml.c_str());
- Settings::data_assimilation_parameters da = Settings::parse_assimilation_parameters(doc.RootElement());
- EXPECT_TRUE(da.enabled);
- EXPECT_EQ("TemperatureSourceChanger", da.class_name);
- EXPECT_TRUE(da.load_data);
- EXPECT_EQ(da.file, "1.10000e+01");
- EXPECT_EQ(da.output_dir, ".vis");
- EXPECT_EQ(da.output_time_interval, 1);
- EXPECT_EQ(da.time, 11);
- EXPECT_EQ(da.port, 10);
-}
-
-TEST(SettingsTest, assimilation4) {
- std::string xml = R"(
-
-
+
1
10
@@ -1812,7 +1858,6 @@ TEST(SettingsTest, assimilation4) {
doc.Parse(xml.c_str());
Settings::data_assimilation_parameters da = Settings::parse_assimilation_parameters(doc.RootElement());
EXPECT_TRUE(da.enabled);
- EXPECT_EQ("ObstacleChanger", da.class_name);
EXPECT_TRUE(da.load_data);
EXPECT_EQ(da.file, "1.10000e+01");
EXPECT_EQ(da.output_dir, ".vis");
@@ -1839,7 +1884,7 @@ TEST(SettingsTest, assimilationHeatSourceChanges) {
)";
tinyxml2::XMLDocument doc;
doc.Parse(xml.c_str());
- Settings::data_assimilation::field_changes field_changes = Settings::parse_field_changes(doc.RootElement(), "temperature_source");
+ Settings::data_assimilation::field_changes field_changes = Settings::parse_field_changes(doc.RootElement());
EXPECT_FALSE(field_changes.u_changed);
EXPECT_FALSE(field_changes.v_changed);
EXPECT_FALSE(field_changes.w_changed);
@@ -1868,11 +1913,21 @@ TEST(SettingsTest, assimilationHeatSourceChanges) {
TEST(SettingsTest, assimilationFieldChanges) {
std::string xml = R"(
-
+
+
+
+
+
+
)";
tinyxml2::XMLDocument doc;
doc.Parse(xml.c_str());
- Settings::data_assimilation::field_changes field_changes = Settings::parse_field_changes(doc.RootElement(), "temperature");
+ auto methods = Settings::parse_data_assimilation_methods(doc.RootElement());
+ EXPECT_EQ(std::get<0>(methods[0]), "FieldChanger");
+ std::string tag = std::get<1>(methods[0]);
+ EXPECT_EQ(tag, "field_changer");
+ auto subsection = Settings::get_subsection(tag, doc.RootElement());
+ Settings::data_assimilation::field_changes field_changes = Settings::parse_field_changes(subsection);
EXPECT_FALSE(field_changes.u_changed);
EXPECT_FALSE(field_changes.v_changed);
EXPECT_FALSE(field_changes.w_changed);
@@ -1880,3 +1935,251 @@ TEST(SettingsTest, assimilationFieldChanges) {
EXPECT_FALSE(field_changes.T_changed);
EXPECT_FALSE(field_changes.C_changed);
}
+
+TEST(SettingsTest, assimilationChanges) {
+ std::string xml = R"(
+
+
+
+
+
+
+
+)";
+ tinyxml2::XMLDocument doc;
+ doc.Parse(xml.c_str());
+ auto methods = Settings::parse_data_assimilation_methods(doc.RootElement());
+ EXPECT_EQ(std::get<0>(methods[0]), "ObstacleChanger");
+ EXPECT_EQ(std::get<1>(methods[0]), "obstacle_changer");
+ EXPECT_EQ(std::get<0>(methods[1]), "TemperatureSourceChanger");
+ EXPECT_EQ(std::get<1>(methods[1]), "temperature_changer");
+}
+
+TEST(SettingsTest, assimilationChanges2) {
+ std::string xml = R"(
+
+
+
+
+
+
+
+
+ 50.3
+ 1.
+ 1.4
+ 0.02
+ 0.
+ 0.15
+ 0.6
+ 0.1
+ 5.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+)";
+ tinyxml2::XMLDocument doc;
+ doc.Parse(xml.c_str());
+ auto field_changes = Settings::parse_field_changes(doc.RootElement());
+ auto da_methods = Settings::parse_data_assimilation_methods(doc.RootElement());
+ Settings::solver::temperature_source temperature_source;
+ Settings::obstacles_parameters obstacle_parameters{ };
+ std::vector names_of_deleted_obstacles;
+ for (std::tuple tuple: da_methods) {
+ std::string name = std::get<0>(tuple);
+ std::string tag = std::get<1>(tuple);
+ if (name == "TemperatureSourceChanger") {
+ auto subsection = Settings::get_subsection(tag, doc.RootElement());
+ temperature_source = Settings::solver::parse_temperature_source(subsection, tag);
+ } else if (name == "ObstacleChanger") {
+ auto subsection = Settings::get_subsection(tag, doc.RootElement());
+ for (const auto *i = subsection->FirstChildElement(); i; i = i->NextSiblingElement()) {
+ Settings::obstacle o = Settings::parse_obstacle(i, tag);
+ if (o.state != State::DELETED) {
+ obstacle_parameters.obstacles.push_back(o);
+ } else {
+ names_of_deleted_obstacles.push_back(o.name);
+ }
+ }
+ int counter_unmodified = 0;
+ int counter_deleted = 0;
+ int counter_new = 0;
+ int counter_modified = 0;
+ int counter_xml = 0;
+ for (const auto &obstacle: obstacle_parameters.obstacles) {
+ if (obstacle.state == State::UNMODIFIED) {
+ counter_unmodified++;
+ } else if (obstacle.state == State::DELETED) {
+ counter_deleted++;
+ } else if (obstacle.state == State::MODIFIED) {
+ counter_modified++;
+ } else if (obstacle.state == State::NEW) {
+ counter_new++;
+ } else if (obstacle.state == State::XML) {
+ counter_xml++;
+ }
+ }
+ bool parameter_changes = counter_deleted + counter_new + counter_modified > 0;
+ EXPECT_EQ(counter_new, 1);
+ EXPECT_EQ(counter_modified, 1);
+ EXPECT_EQ(counter_unmodified, 0);
+ EXPECT_EQ(counter_deleted, 0);
+ EXPECT_EQ(counter_xml, 1);
+ EXPECT_TRUE(parameter_changes);
+ obstacle_parameters.enabled = counter_new + counter_modified + counter_unmodified + counter_xml > 0;
+ }
+ }
+ auto methods = Settings::parse_data_assimilation_methods(doc.RootElement());
+ EXPECT_EQ(std::get<0>(methods[0]), "ObstacleChanger");
+ EXPECT_EQ(std::get<1>(methods[0]), "obstacle_changer");
+ EXPECT_EQ(std::get<0>(methods[1]), "TemperatureSourceChanger");
+ EXPECT_EQ(std::get<1>(methods[1]), "temperature_changer");
+
+ // temperature source
+ EXPECT_FALSE(temperature_source.dissipation);
+ EXPECT_EQ(temperature_source.temp_fct, SourceMethods::Gauss);
+ const auto &gauss_temp = std::get(temperature_source.temp_function);
+ EXPECT_DOUBLE_EQ(gauss_temp.tau, 5);
+ EXPECT_DOUBLE_EQ(gauss_temp.heat_release_rate, 50.3);
+ EXPECT_DOUBLE_EQ(gauss_temp.heat_capacity, 1);
+ EXPECT_DOUBLE_EQ(gauss_temp.position[CoordinateAxis::X], 1.4);
+ EXPECT_DOUBLE_EQ(gauss_temp.position[CoordinateAxis::Y], 0.02);
+ EXPECT_DOUBLE_EQ(gauss_temp.position[CoordinateAxis::Z], 0);
+ EXPECT_DOUBLE_EQ(gauss_temp.dimension[CoordinateAxis::X], 0.15);
+ EXPECT_DOUBLE_EQ(gauss_temp.dimension[CoordinateAxis::Y], 0.6);
+ EXPECT_DOUBLE_EQ(gauss_temp.dimension[CoordinateAxis::Z], 0.1);
+
+ // obstacle
+ EXPECT_TRUE(obstacle_parameters.enabled);
+ EXPECT_EQ(obstacle_parameters.obstacles.size(), 3);
+
+ const auto &obstacle1 = obstacle_parameters.obstacles[0];
+ EXPECT_EQ(obstacle1.name, "ceiling");
+ EXPECT_EQ(obstacle1.state, State::MODIFIED);
+ EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::X], 0);
+ EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::Y], 2.18);
+ EXPECT_DOUBLE_EQ(obstacle1.start_coords[CoordinateAxis::Z], -1.4);
+ EXPECT_DOUBLE_EQ(obstacle1.end_coords[CoordinateAxis::X], 2.8);
+ EXPECT_DOUBLE_EQ(obstacle1.end_coords[CoordinateAxis::Y], 2.38);
+ EXPECT_DOUBLE_EQ(obstacle1.end_coords[CoordinateAxis::Z], 1.4);
+
+ EXPECT_EQ(obstacle1.boundaries.size(), 2);
+ const auto &o1_b1 = obstacle1.boundaries[0];
+ EXPECT_DOUBLE_EQ(o1_b1.value.value(), -1);
+ EXPECT_EQ(o1_b1.boundary_condition, BoundaryCondition::DIRICHLET);
+ EXPECT_NE(std::find(o1_b1.field_type.begin(),
+ o1_b1.field_type.end(),
+ FieldType::U),
+ o1_b1.field_type.end());
+ EXPECT_NE(std::find(o1_b1.field_type.begin(),
+ o1_b1.field_type.end(),
+ FieldType::V),
+ o1_b1.field_type.end());
+ EXPECT_NE(std::find(o1_b1.field_type.begin(),
+ o1_b1.field_type.end(),
+ FieldType::W),
+ o1_b1.field_type.end());
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ EXPECT_NE(std::find(o1_b1.patch.begin(),
+ o1_b1.patch.end(),
+ patch),
+ o1_b1.patch.end());
+ }
+ const auto &o1_b2 = obstacle1.boundaries[1];
+ EXPECT_DOUBLE_EQ(o1_b2.value.value(), 0);
+ EXPECT_EQ(o1_b2.boundary_condition, BoundaryCondition::NEUMANN);
+ EXPECT_NE(std::find(o1_b2.field_type.begin(),
+ o1_b2.field_type.end(),
+ FieldType::P),
+ o1_b2.field_type.end());
+ EXPECT_NE(std::find(o1_b2.field_type.begin(),
+ o1_b2.field_type.end(),
+ FieldType::T),
+ o1_b2.field_type.end());
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ EXPECT_NE(std::find(o1_b2.patch.begin(),
+ o1_b2.patch.end(),
+ patch),
+ o1_b2.patch.end());
+ }
+
+ const auto &obstacle2 = obstacle_parameters.obstacles[1];
+ EXPECT_EQ(obstacle2.name, "right wall left from door");
+ EXPECT_EQ(obstacle2.state, State::NEW);
+ EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::X], 2.8);
+ EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::Y], 0);
+ EXPECT_DOUBLE_EQ(obstacle2.start_coords[CoordinateAxis::Z], -1.4);
+ EXPECT_DOUBLE_EQ(obstacle2.end_coords[CoordinateAxis::X], 3.2);
+ EXPECT_DOUBLE_EQ(obstacle2.end_coords[CoordinateAxis::Y], 2.38);
+ EXPECT_DOUBLE_EQ(obstacle2.end_coords[CoordinateAxis::Z], -0.43);
+
+ EXPECT_EQ(obstacle2.boundaries.size(), 2);
+ const auto &o2_b1 = obstacle2.boundaries[0];
+ EXPECT_EQ(o2_b1.boundary_condition, BoundaryCondition::PERIODIC);
+ EXPECT_NE(std::find(o2_b1.field_type.begin(),
+ o2_b1.field_type.end(),
+ FieldType::U),
+ o2_b1.field_type.end());
+ EXPECT_NE(std::find(o2_b1.field_type.begin(),
+ o2_b1.field_type.end(),
+ FieldType::V),
+ o2_b1.field_type.end());
+ EXPECT_NE(std::find(o2_b1.field_type.begin(),
+ o2_b1.field_type.end(),
+ FieldType::P),
+ o2_b1.field_type.end());
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ EXPECT_NE(std::find(o2_b1.patch.begin(),
+ o2_b1.patch.end(),
+ patch),
+ o2_b1.patch.end());
+ }
+ const auto &o2_b2 = obstacle2.boundaries[1];
+ EXPECT_DOUBLE_EQ(o2_b2.value.value(), 10);
+ EXPECT_EQ(o2_b2.boundary_condition, BoundaryCondition::DIRICHLET);
+ EXPECT_NE(std::find(o2_b2.field_type.begin(),
+ o2_b2.field_type.end(),
+ FieldType::W),
+ o2_b2.field_type.end());
+ EXPECT_NE(std::find(o2_b2.field_type.begin(),
+ o2_b2.field_type.end(),
+ FieldType::T),
+ o2_b2.field_type.end());
+ for (Patch patch: {Patch::LEFT, Patch::RIGHT, Patch::BOTTOM, Patch::TOP, Patch::FRONT, Patch::BACK}) {
+ EXPECT_NE(std::find(o2_b2.patch.begin(),
+ o2_b2.patch.end(),
+ patch),
+ o2_b2.patch.end());
+ }
+
+ const auto &obstacle3 = obstacle_parameters.obstacles[2];
+ EXPECT_EQ(obstacle3.name, "test of behaviour");
+ EXPECT_EQ(obstacle3.state, State::XML);
+ EXPECT_DOUBLE_EQ(obstacle3.start_coords[CoordinateAxis::X], 2.8);
+ EXPECT_DOUBLE_EQ(obstacle3.start_coords[CoordinateAxis::Y], 0);
+ EXPECT_DOUBLE_EQ(obstacle3.start_coords[CoordinateAxis::Z], -1.4);
+ EXPECT_DOUBLE_EQ(obstacle3.end_coords[CoordinateAxis::X], 3.2);
+ EXPECT_DOUBLE_EQ(obstacle3.end_coords[CoordinateAxis::Y], 2.38);
+ EXPECT_DOUBLE_EQ(obstacle3.end_coords[CoordinateAxis::Z], -0.43);
+}
diff --git a/src/dataAssimilation/CMakeLists.txt b/src/dataAssimilation/CMakeLists.txt
index f163d832..27fe6937 100644
--- a/src/dataAssimilation/CMakeLists.txt
+++ b/src/dataAssimilation/CMakeLists.txt
@@ -4,10 +4,6 @@ target_sources(artss
${CMAKE_CURRENT_LIST_DIR}/DataAssimilation.h
${CMAKE_CURRENT_LIST_DIR}/FieldIO.cpp
${CMAKE_CURRENT_LIST_DIR}/FieldIO.h
- ${CMAKE_CURRENT_LIST_DIR}/ObstacleChanger.cpp
- ${CMAKE_CURRENT_LIST_DIR}/ObstacleChanger.h
- ${CMAKE_CURRENT_LIST_DIR}/TemperatureSourceChanger.cpp
- ${CMAKE_CURRENT_LIST_DIR}/TemperatureSourceChanger.h
${CMAKE_CURRENT_LIST_DIR}/ParameterReader.cpp
${CMAKE_CURRENT_LIST_DIR}/ParameterReader.h
)
diff --git a/src/dataAssimilation/DataAssimilation.cpp b/src/dataAssimilation/DataAssimilation.cpp
index 3c5aeb08..6a6ff85b 100644
--- a/src/dataAssimilation/DataAssimilation.cpp
+++ b/src/dataAssimilation/DataAssimilation.cpp
@@ -6,17 +6,14 @@
#include "DataAssimilation.h"
#include "../domain/DomainData.h"
-#include "../TCP/TCPServer.h"
#include "mpi.h"
-#include "TemperatureSourceChanger.h"
-#include "ObstacleChanger.h"
#include "../domain/DomainController.h"
DataAssimilation::DataAssimilation(const SolverController &solver_controller,
FieldController *field_controller,
const Settings::Settings &settings) :
m_logger(Utility::create_logger(typeid(this).name())),
- m_settings(settings),
+ m_settings(settings.assimilation_parameters),
m_field_controller(field_controller),
m_solver_controller(solver_controller),
m_new_field_u(Field(FieldType::U)),
@@ -26,23 +23,13 @@ DataAssimilation::DataAssimilation(const SolverController &solver_controller,
m_new_field_T(Field(FieldType::T)),
m_new_field_C(Field(FieldType::RHO)) {
m_time_interval_counter = 0;
- m_output_time_interval = m_settings.assimilation_parameters.output_time_interval;
+ m_output_time_interval = m_settings.output_time_interval;
m_field_IO_handler = new FieldIO(settings.filename,
- m_settings.assimilation_parameters.output_dir);
+ m_settings.output_dir);
- if (m_settings.assimilation_parameters.enabled) {
- if (m_settings.assimilation_parameters.class_name == AssimilationMethods::standard) {
- m_parameter_handler = new ParameterReader();
- } else if (m_settings.assimilation_parameters.class_name == AssimilationMethods::temperature_source) {
- m_parameter_handler = new TemperatureSourceChanger(m_solver_controller,
- m_settings.solver_parameters.temperature.source);
- } else if (m_settings.assimilation_parameters.class_name == AssimilationMethods::obstacle_changer) {
- m_parameter_handler = new ObstacleChanger(m_solver_controller, m_settings.obstacles_parameters);
- } else {
- m_logger->error("assimilation method {} not known", m_settings.assimilation_parameters.class_name);
- std::exit(1);
- }
+ if (m_settings.enabled) {
+ m_parameter_handler = new ParameterReader(m_solver_controller);
}
}
@@ -54,7 +41,7 @@ void DataAssimilation::initiate_time_skip(const real t_cur) {
}
void DataAssimilation::save_data(real t_cur) {
- if (t_cur >= m_output_time_interval * m_time_interval_counter) {
+ if (t_cur >= m_output_time_interval * m_time_interval_counter - DomainData::getInstance()->get_physical_parameters().dt/2) {
m_logger->debug("save data for {} with interval {} counter {}",
t_cur, m_output_time_interval,
m_time_interval_counter);
@@ -107,16 +94,30 @@ bool DataAssimilation::config_rollback(const char *msg) {
m_new_field_u, m_new_field_v, m_new_field_w,
m_new_field_p, m_new_field_T, m_new_field_C);
auto domain_controller = DomainController::getInstance();
+ if (field_changes.u_changed) {
+ domain_controller->apply_boundary(m_new_field_u);
+ }
+ if (field_changes.v_changed) {
+ domain_controller->apply_boundary(m_new_field_v);
+ }
+ if (field_changes.w_changed) {
+ domain_controller->apply_boundary(m_new_field_w);
+ }
+ if (field_changes.p_changed) {
+ domain_controller->apply_boundary(m_new_field_p);
+ }
if (field_changes.T_changed) {
domain_controller->apply_boundary(m_new_field_T);
}
- //TODO
+ if (field_changes.C_changed) {
+ domain_controller->apply_boundary(m_new_field_C);
+ }
return changes || field_changes.changed;
}
}
bool DataAssimilation::requires_rollback(const real t_cur) {
- if (!m_settings.assimilation_parameters.enabled) {
+ if (!m_settings.enabled) {
return false;
}
m_t_cur = t_cur;
@@ -146,13 +147,13 @@ bool DataAssimilation::requires_rollback(const real t_cur) {
}
bool DataAssimilation::load_data() {
- bool load_data = m_settings.assimilation_parameters.load_data;
+ bool load_data = m_settings.load_data;
if (!load_data) {
return load_data;
}
- m_t_cur = m_settings.assimilation_parameters.time;
+ m_t_cur = m_settings.time;
- m_field_IO_handler->read_fields(m_settings.assimilation_parameters.file,
+ m_field_IO_handler->read_fields(m_settings.file,
m_t_cur,
m_new_field_u, m_new_field_v, m_new_field_w,
m_new_field_p, m_new_field_T, m_new_field_C);
diff --git a/src/dataAssimilation/DataAssimilation.h b/src/dataAssimilation/DataAssimilation.h
index e5aa3544..c6a84055 100644
--- a/src/dataAssimilation/DataAssimilation.h
+++ b/src/dataAssimilation/DataAssimilation.h
@@ -15,10 +15,9 @@
#include "FieldIO.h"
#include "../solver/SolverController.h"
#include "ParameterReader.h"
-#include "../interfaces/IParameterReader.h"
struct AssimilationMethods {
- inline static const std::string standard = "default";
+ inline static const std::string field_changer = "FieldChanger";
inline static const std::string temperature_source = "TemperatureSourceChanger";
inline static const std::string obstacle_changer = "ObstacleChanger";
};
@@ -44,12 +43,12 @@ class DataAssimilation {
private:
std::shared_ptr m_logger;
- const Settings::Settings &m_settings;
+ const Settings::data_assimilation_parameters &m_settings;
FieldController *m_field_controller;
const SolverController &m_solver_controller;
FieldIO *m_field_IO_handler;
- IParameterReader *m_parameter_handler;
+ ParameterReader *m_parameter_handler;
real m_t_cur = 0;
real m_output_time_interval;
diff --git a/src/dataAssimilation/ObstacleChanger.cpp b/src/dataAssimilation/ObstacleChanger.cpp
deleted file mode 100644
index 89645ecb..00000000
--- a/src/dataAssimilation/ObstacleChanger.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-/// \file ObstacleChanger.cpp
-/// \brief
-/// \date Nov 23, 2022
-/// \author My Linh Wuerzburger
-/// \copyright <2015-2022> Forschungszentrum Juelich. All rights reserved
-
-#include
-#include "ObstacleChanger.h"
-#include "../domain/DomainController.h"
-
-return_parameter_reader ObstacleChanger::read_config(const std::string &filename) {
- try {
- m_logger->debug("parse file to string {}", filename);
- auto file_content = Settings::parse_settings_from_file(filename);
- m_logger->debug("parse document from {} to XMLTree {}", filename, file_content);
- tinyxml2::XMLDocument doc;
- doc.Parse(file_content.c_str());
- m_logger->debug("parse obstacles {}", static_cast(doc.RootElement()));
- auto obstacle_parameters = Settings::parse_obstacles_parameters(doc.RootElement());
- bool parameter_changes = false;
- for (const auto& obstacle: obstacle_parameters.obstacles) {
- if (obstacle.state == State::MODIFIED || obstacle.state == State::NEW || obstacle.state == State::DELETED) {
- parameter_changes = true;
- break;
- }
- }
- if (parameter_changes) {
- m_logger->debug("apply obstacle changes");
- std::chrono::time_point start, end;
- start = std::chrono::system_clock::now();
- DomainController::getInstance()->replace_obstacles(obstacle_parameters);
- end = std::chrono::system_clock::now();
- long ms = std::chrono::duration_cast < std::chrono::milliseconds > (end - start).count();
- m_logger->debug("replace obstacles time: {}", ms);
-#ifndef BENCHMARKING
- int counter_unmodified = 0;
- int counter_deleted = 0;
- int counter_rest = 0;
- for (const auto& obstacle: obstacle_parameters.obstacles) {
- if (obstacle.state == State::UNMODIFIED) {
- counter_unmodified++;
- } else if (obstacle.state == State::DELETED) {
- counter_deleted++;
- } else if (obstacle.state == State::MODIFIED || obstacle.state == State::NEW) {
- counter_rest++;
- } else {
- m_logger->warn("wrong state: {} ({})", obstacle.state, obstacle.name);
- }
- }
- m_logger->debug("changed obstacles, unmodified: {}, modified/new: {}, deleted {}", counter_unmodified, counter_deleted, counter_rest);
-#endif
- m_solver_controller.update_sight();
- start = std::chrono::system_clock::now();
- m_solver_controller.m_solver->update_obstacle_change();
- end = std::chrono::system_clock::now();
- ms = std::chrono::duration_cast < std::chrono::milliseconds > (end - start).count();
- m_logger->debug("update source: {}", ms);
-#ifndef BENCHMARKING
- } else {
- m_logger->debug("no obstacle changes");
-#endif
- }
- m_logger->debug("parse field changes");
-
- auto field_changes = Settings::parse_field_changes(doc.RootElement(), "TemperatureSourceChanger");
- return {parameter_changes, field_changes};
- } catch (const std::exception &ex) {
- std::cerr << ex.what() << std::endl;
- }
- Settings::data_assimilation::field_changes field_changes;
- field_changes.changed = false;
- return {false, field_changes};
-}
diff --git a/src/dataAssimilation/ObstacleChanger.h b/src/dataAssimilation/ObstacleChanger.h
deleted file mode 100644
index 17d03da9..00000000
--- a/src/dataAssimilation/ObstacleChanger.h
+++ /dev/null
@@ -1,33 +0,0 @@
-/// \file ObstacleChanger.h
-/// \brief
-/// \date Nov 23, 2022
-/// \author My Linh Wuerzburger
-/// \copyright <2015-2022> Forschungszentrum Juelich. All rights reserved
-
-#ifndef ARTSS_DATAASSIMILATION_OBSTACLECHANGER_H
-#define ARTSS_DATAASSIMILATION_OBSTACLECHANGER_H
-
-
-#include "../interfaces/IParameterReader.h"
-#include "../utility/settings/Settings.h"
-#include "../utility/Utility.h"
-#include "../solver/SolverController.h"
-
-class ObstacleChanger : public IParameterReader {
-public:
- ObstacleChanger(const SolverController &solver_controller,
- const Settings::obstacles_parameters &obstacle_parameters) :
- m_solver_controller(solver_controller),
- m_obstacles(obstacle_parameters),
- m_logger(Utility::create_logger(typeid(this).name())) { }
-
- return_parameter_reader read_config(const std::string &filename) override;
-
-private:
- const SolverController &m_solver_controller;
- const Settings::obstacles_parameters &m_obstacles;
- std::shared_ptr m_logger;
-};
-
-
-#endif /* ARTSS_DATAASSIMILATION_OBSTACLECHANGER_H */
diff --git a/src/dataAssimilation/ParameterReader.cpp b/src/dataAssimilation/ParameterReader.cpp
index ed6ef1eb..4169a1f2 100644
--- a/src/dataAssimilation/ParameterReader.cpp
+++ b/src/dataAssimilation/ParameterReader.cpp
@@ -4,22 +4,120 @@
/// \author My Linh Wuerzburger
/// \copyright <2015-2021> Forschungszentrum Juelich. All rights reserved
+#include
#include "ParameterReader.h"
+#include "DataAssimilation.h"
+#include "../domain/DomainController.h"
return_parameter_reader ParameterReader::read_config(const std::string &file_name) {
+ bool parameter_changes = false;
+ Settings::data_assimilation::field_changes field_changes{ };
+ field_changes.changed = false;
try {
- m_logger->debug("parse file to string");
+ m_logger->debug("parse file to string {}", file_name);
auto file_content = Settings::parse_settings_from_file(file_name);
- m_logger->debug("parse document to XMLTree");
+ m_logger->debug("parse document to XMLTree\n{}", file_content);
tinyxml2::XMLDocument doc;
doc.Parse(file_content.c_str());
- m_logger->debug("parse field changes");
- auto field_changes = Settings::parse_field_changes(doc.RootElement(), "field_changes");
- return {true, field_changes};
+ auto da_methods = Settings::parse_data_assimilation_methods(doc.RootElement());
+ for (std::tuple tuple: da_methods) {
+ std::string name = std::get<0>(tuple);
+ std::string tag = std::get<1>(tuple);
+ if (name == AssimilationMethods::field_changer) {
+ auto subsection = Settings::get_subsection(tag, doc.RootElement());
+ field_changes = Settings::parse_field_changes(subsection);
+ } else if (name == AssimilationMethods::temperature_source) {
+ auto subsection = Settings::get_subsection(tag, doc.RootElement());
+ parameter_changes = parameter_changes || temperature_source_changer(subsection, tag);
+ } else if (name == AssimilationMethods::obstacle_changer) {
+ auto subsection = Settings::get_subsection(tag, doc.RootElement());
+ parameter_changes = parameter_changes || obstacle_changer(subsection, tag);
+ } else {
+ m_logger->debug("Unknown Assimilation method: {}", name);
+ }
+ }
} catch (const std::exception &ex) {
std::cerr << ex.what() << std::endl;
}
- Settings::data_assimilation::field_changes field_changes;
- field_changes.changed = false;
- return {false, field_changes} ;
+ return {parameter_changes, field_changes};
+}
+
+bool ParameterReader::temperature_source_changer(const tinyxml2::XMLElement *head, const std::string &context) const {
+ auto temperature_source = Settings::solver::parse_temperature_source(head, context);
+ bool parameter_changes = true;
+ // TODO (c++20)
+ // bool parameter_changes = temperature_source != m_temperature_source;
+ if (parameter_changes) {
+ m_logger->debug("apply heat source changes");
+ m_solver_controller.m_solver->replace_heat_source(temperature_source);
+ }
+ return parameter_changes;
+}
+
+bool ParameterReader::obstacle_changer(const tinyxml2::XMLElement *head, const std::string &context) const {
+ Settings::obstacles_parameters obstacle_parameters{ };
+ std::vector names_of_deleted_obstacles;
+ for (const auto *i = head->FirstChildElement(); i; i = i->NextSiblingElement()) {
+ Settings::obstacle o = Settings::parse_obstacle(i, context);
+ if (o.state != State::DELETED) {
+ obstacle_parameters.obstacles.push_back(o);
+ } else {
+ names_of_deleted_obstacles.push_back(o.name);
+ }
+ }
+ size_t counter_deleted = names_of_deleted_obstacles.size();
+ size_t counter_unmodified = 0;
+ size_t counter_new = 0;
+ size_t counter_modified = 0;
+ size_t counter_xml = 0;
+ for (const auto &obstacle: obstacle_parameters.obstacles) {
+ if (obstacle.state == State::UNMODIFIED) {
+ counter_unmodified++;
+ } else if (obstacle.state == State::MODIFIED) {
+ counter_modified++;
+ } else if (obstacle.state == State::NEW) {
+ counter_new++;
+ } else if (obstacle.state == State::XML) {
+ counter_xml++;
+ }
+ }
+#ifndef BENCHMARKING
+ m_logger->debug("obstacles read in:\n "
+ "unmodified: {}\n "
+ "modified: {}\n "
+ "deleted: {}\n "
+ "new: {}\n "
+ "XML: {}",
+ counter_unmodified,
+ counter_modified,
+ counter_deleted,
+ counter_new,
+ counter_xml);
+#endif
+ bool parameter_changes = counter_deleted + counter_new + counter_modified > 0;
+ obstacle_parameters.enabled = counter_new + counter_modified + counter_unmodified + counter_xml > 0;
+ if (parameter_changes) {
+#ifndef BENCHMARKING
+ m_logger->debug("apply obstacle changes");
+ std::chrono::time_point start, end;
+ start = std::chrono::system_clock::now();
+#endif
+ DomainController::getInstance()->replace_obstacles(obstacle_parameters);
+#ifndef BENCHMARKING
+ end = std::chrono::system_clock::now();
+ long ms = std::chrono::duration_cast(end - start).count();
+ m_logger->debug("replace obstacles time: {}", ms);
+#endif
+ m_solver_controller.update_sight();
+#ifndef BENCHMARKING
+ start = std::chrono::system_clock::now();
+ m_solver_controller.m_solver->update_obstacle_change();
+ end = std::chrono::system_clock::now();
+ ms = std::chrono::duration_cast(end - start).count();
+ m_logger->debug("update source: {}", ms);
+ } else {
+ m_logger->debug("no obstacle changes");
+#endif
+ }
+ return parameter_changes;
}
diff --git a/src/dataAssimilation/ParameterReader.h b/src/dataAssimilation/ParameterReader.h
index d780b81d..0cd0929e 100644
--- a/src/dataAssimilation/ParameterReader.h
+++ b/src/dataAssimilation/ParameterReader.h
@@ -9,18 +9,27 @@
#include
-#include "../interfaces/IParameterReader.h"
#include "../utility/settings/Settings.h"
#include "../utility/Utility.h"
+#include "../solver/SolverController.h"
+
+using return_parameter_reader = std::tuple;
+
+class ParameterReader {
+public:
+ explicit ParameterReader(const SolverController &solver_controller) : m_solver_controller(solver_controller), m_logger(Utility::create_logger(typeid(this).name())) { }
-class ParameterReader : public IParameterReader {
- public:
- ParameterReader() : m_logger(Utility::create_logger(typeid(this).name())) {}
~ParameterReader() = default;
- return_parameter_reader read_config(const std::string &file_name) override;
- private:
+ return_parameter_reader read_config(const std::string &file_name);
+
+private:
+ const SolverController &m_solver_controller;
std::shared_ptr m_logger;
+
+ bool temperature_source_changer(const tinyxml2::XMLElement *doc, const std::string &context) const;
+
+ bool obstacle_changer(const tinyxml2::XMLElement *head, const std::string &context) const;
};
#endif /* ARTSS_DATAASSIMILATION_PARAMETERREADER_H */
diff --git a/src/dataAssimilation/TemperatureSourceChanger.cpp b/src/dataAssimilation/TemperatureSourceChanger.cpp
deleted file mode 100644
index d6807287..00000000
--- a/src/dataAssimilation/TemperatureSourceChanger.cpp
+++ /dev/null
@@ -1,34 +0,0 @@
-/// \file TemperatureSourceChanger.cpp
-/// \brief Class for changing temperature source
-/// \date Jan 05, 2021
-/// \author My Linh Wuerzburger
-/// \copyright <2015-2022> Forschungszentrum Juelich All rights reserved.
-
-#include "TemperatureSourceChanger.h"
-
-return_parameter_reader TemperatureSourceChanger::read_config(const std::string &filename) {
- try {
- m_logger->debug("parse file to string {}", filename);
- auto file_content = Settings::parse_settings_from_file(filename);
- m_logger->debug("parse document from {} to XMLTree {}", filename, file_content);
- tinyxml2::XMLDocument doc;
- doc.Parse(file_content.c_str());
- m_logger->debug("parse heat source changes {}", static_cast(doc.RootElement()));
- auto temperature_source = Settings::solver::parse_temperature_source(doc.RootElement(), "temperature_source");
- bool parameter_changes = true;
- // TODO (c++20)
- // bool parameter_changes = temperature_source != m_temperature_source;
- if (parameter_changes) {
- m_logger->debug("apply heat source changes");
- m_solver_controller.m_solver->replace_heat_source(temperature_source);
- }
- m_logger->debug("parse field changes");
- auto field_changes = Settings::parse_field_changes(doc.RootElement(), "TemperatureSourceChanger");
- return {parameter_changes, field_changes};
- } catch (const std::exception &ex) {
- std::cerr << ex.what() << std::endl;
- }
- Settings::data_assimilation::field_changes field_changes;
- field_changes.changed = false;
- return {false, field_changes};
-}
diff --git a/src/dataAssimilation/TemperatureSourceChanger.h b/src/dataAssimilation/TemperatureSourceChanger.h
deleted file mode 100644
index 4525f529..00000000
--- a/src/dataAssimilation/TemperatureSourceChanger.h
+++ /dev/null
@@ -1,31 +0,0 @@
-/// \file TemperatureSourceChanger.h
-/// \brief Class for changing temperature source
-/// \date Jan 05, 2021
-/// \author My Linh Wuerzburger
-/// \copyright <2015-2022> Forschungszentrum Juelich All rights reserved.
-
-#ifndef ARTSS_DATAASSIMILATION_TEMPERATURESOURCECHANGER_H
-#define ARTSS_DATAASSIMILATION_TEMPERATURESOURCECHANGER_H
-
-#include
-
-#include "../solver/SolverController.h"
-#include "../interfaces/IParameterReader.h"
-#include "../utility/settings/Settings.h"
-#include "../utility/Utility.h"
-
-class TemperatureSourceChanger : public IParameterReader {
- public:
- TemperatureSourceChanger(const SolverController &solver_controller,
- const Settings::solver::temperature_source &temperature_source) :
- m_solver_controller(solver_controller),
- m_temperature_source(temperature_source),
- m_logger(Utility::create_logger(typeid(this).name())) { }
- return_parameter_reader read_config(const std::string &filename) override;
- private:
- const SolverController &m_solver_controller;
- const Settings::solver::temperature_source &m_temperature_source;
- std::shared_ptr m_logger;
-};
-
-#endif /* ARTSS_DATAASSIMILATION_TEMPERATURESOURCECHANGER_H */
diff --git a/src/domain/DomainController.cpp b/src/domain/DomainController.cpp
index a4a06a99..e16bcd47 100644
--- a/src/domain/DomainController.cpp
+++ b/src/domain/DomainController.cpp
@@ -7,7 +7,7 @@
#include "DomainController.h"
#include
-std::unique_ptr DomainController::single{};
+std::unique_ptr DomainController::single{ };
DomainController::DomainController(const Settings::Settings &settings) :
m_settings(settings) {
@@ -34,7 +34,7 @@ return_xml_objects DomainController::read_XML() {
m_logger->debug("start parsing XML");
m_logger->debug("start parsing boundary parameter");
#endif
- BoundaryDataController bdc_domain(m_settings.boundary_parameters.boundaries);
+ BoundaryDataController bdc_domain(m_settings.boundary_parameters.boundaries);
#ifndef BENCHMARKING
m_logger->debug("finished parsing boundary parameter");
#endif
@@ -88,17 +88,17 @@ return_obstacle DomainController::parse_obstacle_parameter(const Settings::obsta
std::vector obstacles;
std::vector bdc_obstacles;
if (obstacle_settings.enabled) {
- obstacles.reserve(obstacle_settings.obstacles.size());
- bdc_obstacles.reserve(obstacle_settings.obstacles.size());
- for (const Settings::obstacle &o: obstacle_settings.obstacles) {
+ obstacles.reserve(obstacle_settings.obstacles.size());
+ bdc_obstacles.reserve(obstacle_settings.obstacles.size());
+ for (const Settings::obstacle &o: obstacle_settings.obstacles) {
#ifndef BENCHMARKING
- m_logger->debug("read {}", o.name);
- m_logger->debug("start coords {}", o.start_coords);
- m_logger->debug("end coords {}", o.end_coords);
+ m_logger->debug("read {}", o.name);
+ m_logger->debug("start coords {}", o.start_coords);
+ m_logger->debug("end coords {}", o.end_coords);
#endif
- obstacles.emplace_back(o.start_coords, o.end_coords, o.name);
- bdc_obstacles.emplace_back(o.boundaries);
- }
+ obstacles.emplace_back(o.start_coords, o.end_coords, o.name);
+ bdc_obstacles.emplace_back(o.boundaries);
+ }
}
#ifndef BENCHMARKING
m_logger->debug("finished parsing obstacle parameter");
@@ -119,10 +119,10 @@ void DomainController::print_boundaries(const BoundaryDataController &bdc_domain
m_logger->info("-- Info summary");
DomainData::getInstance()->print();
bdc_domain.print();
- for (const auto & bdc_obstacle : bdc_obstacles) {
+ for (const auto &bdc_obstacle: bdc_obstacles) {
bdc_obstacle.print();
}
- for (const auto & bdc_surface : bdc_surfaces) {
+ for (const auto &bdc_surface: bdc_surfaces) {
bdc_surface.print();
}
#endif
@@ -181,41 +181,10 @@ bool DomainController::is_blocked_by_obstacle(const Coordinate &start, c
}
void DomainController::replace_obstacles(const Settings::obstacles_parameters &obstacle_parameters) {
-//#ifndef BENCHMARKING
-// m_logger->debug("start parsing obstacle parameter");
-//#endif
-// std::vector unmodified;
-// std::vector deleted;
-// std::vector obstacles;
-// std::vector bdc_obstacles;
-// if (obstacle_parameters.enabled) {
-// obstacles.reserve(obstacle_parameters.obstacles.size());
-// bdc_obstacles.reserve(obstacle_parameters.obstacles.size());
-// for (const Settings::obstacle &o: obstacle_parameters.obstacles) {
-// switch (o.state) {
-// case State::DELETED:
-// deleted.push_back(o.name);
-// case State::UNMODIFIED:
-// unmodified.push_back(o.name);
-// break;
-// case State::NEW:
-// case State::MODIFIED:
-// obstacles.emplace_back(o.start_coords, o.end_coords, o.name);
-// bdc_obstacles.emplace_back(o.boundaries);
-// break;
-// default:
-// m_logger->warn("obstacle ({}) with unknown state: {}", o.name, o.state);
-// }
-// }
-// }
-//#ifndef BENCHMARKING
-// m_logger->debug("finished parsing obstacle parameter");
-//#endif
-
auto [bdc_domain, surfaces, bdc_surfaces, obstacles2, bdc_obstacles2] = read_XML();
auto [obstacles, bdc_obstacles] = parse_obstacle_parameter(obstacle_parameters);
detect_neighbouring_obstacles(obstacles);
- //m_multigrid->replace_obstacles(obstacles, bdc_obstacles);
+ // m_multigrid->replace_obstacles();
size_t multigrid_level = DomainData::getInstance()->get_levels();
delete m_multigrid;
m_multigrid = new Multigrid(surfaces, bdc_surfaces,
diff --git a/src/domain/Obstacle.h b/src/domain/Obstacle.h
index 18fa0818..3e1d8d84 100644
--- a/src/domain/Obstacle.h
+++ b/src/domain/Obstacle.h
@@ -66,7 +66,7 @@ class Obstacle {
[[nodiscard]] const Coordinate &get_end_coordinates() const { return m_end; }
[[nodiscard]] const Coordinate &get_strides() const { return m_strides; }
- std::string get_name() { return m_name; }
+ std::string get_name() const { return m_name; }
void replace_patch(size_t *indices, size_t size, Patch p);
void control();
diff --git a/src/interfaces/IParameterReader.h b/src/interfaces/IParameterReader.h
deleted file mode 100644
index 57fc3fa7..00000000
--- a/src/interfaces/IParameterReader.h
+++ /dev/null
@@ -1,19 +0,0 @@
-/// \file IParameterReader.h
-/// \brief
-/// \date Jan 05, 2022
-/// \author My Linh Wuerzburger
-/// \copyright <2015-2021> Forschungszentrum Juelich. All rights reserved
-
-#ifndef ARTSS_INTERFACES_IPARAMETERREADER_H
-#define ARTSS_INTERFACES_IPARAMETERREADER_H
-
-#include
-#include "../utility/settings/Settings.h"
-using return_parameter_reader = std::tuple;
-
-class IParameterReader {
-public:
- virtual return_parameter_reader read_config(const std::string &filename) = 0;
-};
-
-#endif /* ARTSS_INTERFACES_IPARAMETERREADER_H */
diff --git a/src/utility/Mapping.h b/src/utility/Mapping.h
index 3755fb1b..b3adfb4c 100644
--- a/src/utility/Mapping.h
+++ b/src/utility/Mapping.h
@@ -15,6 +15,22 @@
inline static const std::vector state_name = {"XML", "unmodified", "modified", "new", "deleted"};
constexpr size_t number_of_states = 5;
+/// \brief 5 different states for obstacles, negligible if no data assimilation is used. State
+/// always refers to the difference between the current ARTSS state and the config file read in.
+/// \details
+/// XML: state of an obstacle from the original XML. (replaces current obstacle with the data from
+/// the XML, not meant for later usage)\n
+/// UNMODIFIED: obstacle already exist in ARTSS, no changes necessary. Usage especially for
+/// obstacles which were at least once modified/newly created\n
+/// MODIFIED: obstacle is to be changed\n
+/// NEW: obstacle doesn't exist in ARTSS and has to be newly created. This also counts for
+/// obstacles which were defined in the XML, got deleted and now need to be newly created\n
+/// DELETED: obstacle exists in ARTSS and has to be deleted\n
+/// \example At the start of ARTSS, obstacles have the state XML. The first obstacle changes happen:
+/// an obstacle is deleted (state: deleted) and another one changed (state: modified). With the
+/// second obstacle change the deleted obstacle is to be restored. The deleted obstacle has to be
+/// newly created (state:new) and and the previously changed obstacle gets the status unmodified
+/// because it is not changed any further.
enum State : int {
UNKNOWN_STATE = -1,
XML = 0,
diff --git a/src/utility/settings/Settings.cpp b/src/utility/settings/Settings.cpp
index f61b5d15..9f2127eb 100644
--- a/src/utility/settings/Settings.cpp
+++ b/src/utility/settings/Settings.cpp
@@ -15,6 +15,7 @@
#include "../../solver/SolverSelection.h"
#include
#include
+#include
namespace Settings {
static const std::string xml_true = "Yes";
@@ -412,17 +413,18 @@ namespace Settings {
obstacle.name = get_required_string(values, "name", context);
obstacle.state = Mapping::match_state(get_optional_string(values, "state", Mapping::get_state_name(State::XML)));
- if (obstacle.state != State::UNMODIFIED) {
- std::string context_geometry = create_context(context, fmt::format("geometry ({})", obstacle.name));
+ std::string sub_context = create_context(context, obstacle.name);
+ if (obstacle.state == State::NEW || obstacle.state == State::MODIFIED || obstacle.state == State::XML) {
+ std::string context_geometry = create_context(sub_context, fmt::format("geometry ({})", obstacle.name));
auto [head_geometry, values_geometry] = map_parameter_section(head, "geometry");
for (size_t a = 0; a < number_of_axes; a++) {
auto axis = CoordinateAxis(a);
std::string axis_name = Utility::to_lower(Mapping::get_axis_name(axis));
- obstacle.start_coords[axis] = get_required_real(values_geometry, "o" + axis_name + "1", context);
- obstacle.end_coords[axis] = get_required_real(values_geometry, "o" + axis_name + "2", context);
+ obstacle.start_coords[axis] = get_required_real(values_geometry, "o" + axis_name + "1", sub_context);
+ obstacle.end_coords[axis] = get_required_real(values_geometry, "o" + axis_name + "2", sub_context);
}
- std::string context_boundaries = create_context(context, fmt::format("boundaries ({})", obstacle.name));
+ std::string context_boundaries = create_context(sub_context, fmt::format("boundaries ({})", obstacle.name));
for (const auto *i = head->FirstChildElement(); i; i = i->NextSiblingElement()) {
if (i->Name() == std::string("boundary")) {
obstacle.boundaries.emplace_back(parse_boundary(i, context_boundaries));
@@ -443,6 +445,17 @@ namespace Settings {
}
}
op.obstacles.shrink_to_fit();
+ std::unordered_set obstacle_names;
+ for (const obstacle &o: op.obstacles) {
+ if (obstacle_names.find(o.name) == obstacle_names.end()) {
+ obstacle_names.insert(o.name);
+ } else {
+ throw config_error(fmt::format("obstacle names have to be unique. Duplicate: {}", o.name));
+ }
+ if (o.state == State::UNKNOWN_STATE) {
+ throw config_error(fmt::format("obstacle state of {} is unknown", o.name));
+ }
+ }
return op;
}
@@ -976,7 +989,6 @@ namespace Settings {
ap.output_dir = get_optional_string(values, "output_dir", ".vis");
if (ap.enabled) {
- ap.class_name = get_required_string(values, "class_name", context);
ap.output_time_interval = get_optional_real(values, "write_output", 1);
ap.load_data = get_optional_bool(values, "load_data", false);
if (ap.load_data) {
@@ -1028,10 +1040,9 @@ namespace Settings {
return settings;
}
- data_assimilation::field_changes parse_field_changes(const tinyxml2::XMLElement *head,
- const std::string &parent_context) {
+ data_assimilation::field_changes parse_field_changes(const tinyxml2::XMLElement *head) {
std::string own_context = "fields_changed";
- std::string context = create_context(parent_context, own_context);
+ std::string context = create_context("assimilation_changes", own_context);
auto [subsection, values] = map_parameter_section(head, own_context);
data_assimilation::field_changes field_changes{ };
field_changes.u_changed = get_required_bool(values, "u", context);
@@ -1041,13 +1052,36 @@ namespace Settings {
field_changes.T_changed = get_required_bool(values, "T", context);
field_changes.C_changed = get_required_bool(values, "concentration", context);
field_changes.changed = field_changes.u_changed || field_changes.v_changed
- || field_changes.w_changed || field_changes.p_changed
- || field_changes.T_changed || field_changes.C_changed;
+ || field_changes.w_changed || field_changes.p_changed
+ || field_changes.T_changed || field_changes.C_changed;
if (field_changes.changed) {
field_changes.file_name = get_required_string(values, "filename", context);
}
return field_changes;
}
+ std::vector> parse_data_assimilation_methods(const tinyxml2::XMLElement *head) {
+ std::vector> methods;
+ std::string own_context = "data_assimilation";
+ std::string context = create_context("assimilation_changes", own_context);
+ auto [subsection, empty_values] = map_parameter_section(head, own_context);
+ for (auto i = subsection->FirstChildElement(); i; i = i->NextSiblingElement()) {
+ auto values = map_parameter_line(i);
+ std::string class_name = get_required_string(values, "name", context);
+ std::string tag_name = get_required_string(values, "tag", context);
+ methods.emplace_back(class_name, tag_name);
+ }
+ return methods;
+ }
+
+}
+
+const tinyxml2::XMLElement *Settings::get_subsection(const std::string &context, const tinyxml2::XMLElement *head) {
+ for (const auto *i = head->FirstChildElement(); i; i = i->NextSiblingElement()) {
+ if (i->Name() == context) {
+ return i;
+ }
+ }
+ throw config_error(fmt::format("No section named '{}' found", context));
}
diff --git a/src/utility/settings/Settings.h b/src/utility/settings/Settings.h
index 75fee1f6..96864f4c 100644
--- a/src/utility/settings/Settings.h
+++ b/src/utility/settings/Settings.h
@@ -334,7 +334,6 @@ namespace Settings {
}
struct data_assimilation_parameters {
bool enabled;
- std::string class_name;
real output_time_interval;
std::string output_dir;
bool load_data;
@@ -359,6 +358,7 @@ namespace Settings {
random_parameters parse_random_parameters(const tinyxml2::XMLElement *head, const std::string &parent_context);
solver_parameters parse_solver_parameters(const tinyxml2::XMLElement *root);
surfaces_parameters parse_surfaces_parameters(const tinyxml2::XMLElement *root);
+ obstacle parse_obstacle(const tinyxml2::XMLElement *root, const std::string &context);
obstacles_parameters parse_obstacles_parameters(const tinyxml2::XMLElement *root);
adaption_parameters parse_adaption_parameters(const tinyxml2::XMLElement *root);
data_assimilation_parameters parse_assimilation_parameters(const tinyxml2::XMLElement *root);
@@ -370,6 +370,9 @@ namespace Settings {
physical_parameters parse_physical_parameters(const tinyxml2::XMLElement *root, const std::string &solver_description);
Settings parse_settings(const std::filesystem::path &path);
std::string parse_settings_from_file(const std::filesystem::path &path);
- data_assimilation::field_changes parse_field_changes(const tinyxml2::XMLElement *head, const std::string &parent_context);
+ data_assimilation::field_changes parse_field_changes(const tinyxml2::XMLElement *head);
+ std::vector> parse_data_assimilation_methods(const tinyxml2::XMLElement *head);
+ const tinyxml2::XMLElement *get_subsection(const std::string &context, const tinyxml2::XMLElement *head);
+
}
#endif
diff --git a/tools/change_xml.sh b/tools/change_xml.sh
index 7bedd6ab..f3f7759f 100755
--- a/tools/change_xml.sh
+++ b/tools/change_xml.sh
@@ -63,6 +63,11 @@ do
shift
shift
;;
+ --z0)
+ z0=$2
+ shift
+ shift
+ ;;
*)
echo "unknown parameter $1"
shift
@@ -108,6 +113,13 @@ then
cp $OUTPUT $TMP
fi
+if [ ! -z $z0 ]
+then
+ echo "change z0 to $z0"
+ sed 's/\s*[0-9]*\.*[0-9]*\s*<\/z0>/ '${z0}' <\/z0>/g' "${TMP}" > "${OUTPUT}"
+ cp $OUTPUT $TMP
+fi
+
if [ ! -z $LOGLEVEL ]
then
echo "change log level to $LOGLEVEL"