USEPA / WNTR

An EPANET compatible python package to simulate and analyze water distribution networks under disaster scenarios.
Other
315 stars 183 forks source link

Skeletonization #128

Closed bahadurr closed 4 years ago

bahadurr commented 4 years ago

I am using WNTRMorphSkel.py script for Skeletonization. The issue is with the aggregation of nodal base demands. Total aggregated base demand for all nodes should remain same even after Skeletonization. But in my test example, the script captures only 10% of the demand. The demands after Skeletonization do not get aggregated and assigned to nodes correctly. I I am using 16 inch pipe diameter for Skeletonization.

Total Demand MGD (aggregate of all nodal base demands) - All pipe model 271.49 MGD Skeltonized Model 30.29 MGD

Thanks for your help

kaklise commented 4 years ago

Are you able to share your script and network model? You can email the files to kaklise@sandia.gov if you would rather not post the files here.

bahadurr commented 4 years ago

Katherine

I can share the script (attached) but the network model I cannot because of obvious reasons. The model input file was given to me by Robert Janke from EPA Cincinnati.

I tried the same script with NET6.inp file that is packaged with WNTR and getting the same results. The script can be tested on NET6.inp file if it helps

Rakesh Bahadur, Ph.D., PMP | Leidos MS FS 1-5-1 11951 Freedom Drive Reston, VA 20190 phone: 571.526.7120 bahadurr@leidos.commailto:bahadurr@leidos.com | leidos.com

[cid:image001.png@01CF9B6D.2D0979D0] This email and any attachments to it are intended only for the identified recipients. It may contain proprietary or otherwise legally protected information of Leidos. Any unauthorized use or disclosure of this communication is strictly prohibited. If you have received this communication in error, please notify the sender and delete or otherwise destroy the email and all attachments immediately.

From: Katherine Klise notifications@github.com Sent: Thursday, March 5, 2020 10:13 AM To: USEPA/WNTR WNTR@noreply.github.com Cc: Bahadur, Rakesh [US-US] RAKESH.BAHADUR@leidos.com; Author author@noreply.github.com Subject: EXTERNAL: Re: [USEPA/WNTR] Skeletonization (#128)

Are you able to share your script and network model? You can email the files to kaklise@sandia.govmailto:kaklise@sandia.gov if you would rather not post the files here.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/USEPA/WNTR/issues/128?email_source=notifications&email_token=ACNPI4HEYZVH6SNNZUISVPTRF66RPA5CNFSM4LCLXZP2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN5UKOI#issuecomment-595281209, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACNPI4EALC76OLUWVI74GMDRF66RPANCNFSM4LCLXZPQ.

bahadurr commented 4 years ago

Katherine

You may not get the attachment as Outlook sometimes strips the python scripts. If that is case then let me know if I can upload the script for you

Rakesh Bahadur, Ph.D., PMP | Leidos MS FS 1-5-1 11951 Freedom Drive Reston, VA 20190 phone: 571.526.7120 bahadurr@leidos.commailto:bahadurr@leidos.com | leidos.com

[cid:image001.png@01CF9B6D.2D0979D0] This email and any attachments to it are intended only for the identified recipients. It may contain proprietary or otherwise legally protected information of Leidos. Any unauthorized use or disclosure of this communication is strictly prohibited. If you have received this communication in error, please notify the sender and delete or otherwise destroy the email and all attachments immediately.

From: Katherine Klise notifications@github.com Sent: Thursday, March 5, 2020 10:13 AM To: USEPA/WNTR WNTR@noreply.github.com Cc: Bahadur, Rakesh [US-US] RAKESH.BAHADUR@leidos.com; Author author@noreply.github.com Subject: EXTERNAL: Re: [USEPA/WNTR] Skeletonization (#128)

Are you able to share your script and network model? You can email the files to kaklise@sandia.govmailto:kaklise@sandia.gov if you would rather not post the files here.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/USEPA/WNTR/issues/128?email_source=notifications&email_token=ACNPI4HEYZVH6SNNZUISVPTRF66RPA5CNFSM4LCLXZP2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN5UKOI#issuecomment-595281209, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACNPI4EALC76OLUWVI74GMDRF66RPANCNFSM4LCLXZPQ.

kaklise commented 4 years ago

I don't see the attached file, but I can provide an example using Net6 to start the discussion.

The issue might be explained by how demands are extracted from the original and skeletonized model. Within the skeletonization process, base demands and demand patterns are moved to remaining nodes. We can compute total demand at each junction and time using base demands, demand patterns, and the demand multiplier to make sure demands are conserved (using wntr.metrics.expected_demand). If we just extract the base demand (using query_node_attribute), the results will be different because this does not account for the multiple demands and patterns on each junction.

In the following example, Net 6 is skeletonized using a 16 inch pipe diameter threshold. The total expected demand is conserved (130.42 m3/s) but the base demands differ. Does this explain the issue you are seeing?

import wntr

# Create a water network model
orig_wn = wntr.network.WaterNetworkModel('networks/Net6.inp')
print(orig_wn.describe())

# Skeletonize the model using a 16 inch pipe diameter threshold
skel_wn = wntr.morph.skeletonize(orig_wn, 16*0.0254)
print(skel_wn.describe())

# Extract the expected demand for each water network model
# This uses base demand, demand pattern, and the demand multiplier
# Expected demand for both models should be the same
orig_expected_demand = wntr.metrics.expected_demand(orig_wn)
skel_expected_demand = wntr.metrics.expected_demand(skel_wn)

orig_expected_demand.sum(axis=1).plot() # plot timeseries of total demand
skel_expected_demand.sum(axis=1).plot()

print(orig_expected_demand.sum().sum()) # print total demand
print(skel_expected_demand.sum().sum())

# Extract the base demand for each water network model
# The base demand will not be the same
orig_base_demand = orig_wn.query_node_attribute('base_demand')
skel_base_demand = skel_wn.query_node_attribute('base_demand')

print(orig_base_demand.sum()) # print total base demand
print(skel_base_demand.sum())
bahadurr commented 4 years ago

Katherine

I changed the extension of the script from .py to .txt and resending. Hopefully this file won’t be stripped by Outlook. You may want to change the extension back to .y to run it.

My understanding is that the Base demand for the whole model (not individual nodes) should remain same before and after Skeletonization.

Rakesh Bahadur, Ph.D., PMP | Leidos MS FS 1-5-1 11951 Freedom Drive Reston, VA 20190 phone: 571.526.7120 bahadurr@leidos.commailto:bahadurr@leidos.com | leidos.com

[cid:image001.png@01CF9B6D.2D0979D0] This email and any attachments to it are intended only for the identified recipients. It may contain proprietary or otherwise legally protected information of Leidos. Any unauthorized use or disclosure of this communication is strictly prohibited. If you have received this communication in error, please notify the sender and delete or otherwise destroy the email and all attachments immediately.

From: Katherine Klise notifications@github.com Sent: Thursday, March 5, 2020 10:59 AM To: USEPA/WNTR WNTR@noreply.github.com Cc: Bahadur, Rakesh [US-US] RAKESH.BAHADUR@leidos.com; Author author@noreply.github.com Subject: EXTERNAL: Re: [USEPA/WNTR] Skeletonization (#128)

I don't see the attached file, but I can provide an example using Net6 to start the discussion.

The issue might be explained by how demands are extracted from the original and skeletonized model. Within the skeletonization process, base demands and demand patterns are moved to remaining nodes. We can compute total demand at each junction and time using base demands, demand patterns, and the demand multiplier to make sure demands are conserved (using wntr.metrics.expected_demand). If we just extract the base demand (using query_node_attribute), the results will be different because this does not account for the multiple demands and patterns on each junction.

In the following example, Net 6 is skeletonized using a 16 inch pipe diameter threshold. The total expected demand is conserved (130.42 m3/s) but the base demands differ. Does this explain the issue you are seeing?

import wntr

Create a water network model

orig_wn = wntr.network.WaterNetworkModel('networks/Net6.inp')

print(orig_wn.describe())

Skeletonize the model using a 16 inch pipe diameter threshold

skel_wn = wntr.morph.skeletonize(orig_wn, 16*0.0254)

print(skel_wn.describe())

Extract the expected demand for each water network model

This uses base demand, demand pattern, and the demand multiplier

Expected demand for both models should be the same

orig_expected_demand = wntr.metrics.expected_demand(orig_wn)

skel_expected_demand = wntr.metrics.expected_demand(skel_wn)

orig_expected_demand.sum(axis=1).plot() # plot timeseries of total demand

skel_expected_demand.sum(axis=1).plot()

print(orig_expected_demand.sum().sum()) # print total demand

print(skel_expected_demand.sum().sum())

Extract the base demand for each water network model

The base demand will not be the same

orig_base_demand = orig_wn.query_node_attribute('base_demand')

skel_base_demand = skel_wn.query_node_attribute('base_demand')

print(orig_base_demand.sum()) # print total base demand

print(skel_base_demand.sum())

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/USEPA/WNTR/issues/128?email_source=notifications&email_token=ACNPI4GVDJEJRYXPR4SXINDRF7D2TA5CNFSM4LCLXZP2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN52CRY#issuecomment-595304775, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACNPI4CC7INRRK7HEDLYOQTRF7D2TANCNFSM4LCLXZPQ.

import numpy as np import matplotlib.pyplot as plt import pickle import wntr """ The wntr.morph.skel module contains functions to skeletonize water network models. """ import logging import itertools import networkx as nx

from wntr.network.elements import Pipe, Junction from wntr.sim.core import WNTRSimulator from wntr.sim import EpanetSimulator from wntr.morph.node import _deepcopy_wn

logger = logging.getLogger(name)

Create a water network model

inp_file = 'C:/Users/bahadurr/Documents/Projects 2019/eCOE/Test Cases/Skeleonization/Texas_Example2_Testing_FOUO.inp' wn = wntr.network.WaterNetworkModel(inp_file)

def skeletonize(wn, pipe_diameter_threshold, branch_trim=True, series_pipe_merge=True, parallel_pipe_merge=True, max_cycles=None, use_epanet=True, return_map=False):

"""
Perform network skeletonization using branch trimming, series pipe merge, 
and parallel pipe merge operations. Candidate pipes for removal is based 
on a pipe diameter threshold.  

Parameters
-------------
wn: wntr WaterNetworkModel
    A WaterNetworkModel object

pipe_diameter_threshold: float 
    Pipe diameter threshold used to determine candidate pipes for skeletonization

branch_trim: bool (optional, default = True)
    Include branch trimming in skeletonization

series_pipe_merge: bool (optional, default = True)
    Include series pipe merge in skeletonization

parallel_pipe_merge: bool (optional, default = True)
    Include parallel pipe merge in skeletonization

max_cycles: int or None (optional, default = None)
    Defines the maximum number of cycles in the skeletonization process. 
    One cycle performs branch trimming for all candidate pipes, followed
    by series pipe merging for all candidate pipes, followed by parallel 
    pipe merging for all candidate pipes. If max_cycles is set to None, 
    skeletonization will run until the network can no longer be reduced.

use_epanet: bool (optional)
    If True, use the EpanetSimulator to compute headloss in pipes.  If False, 
    use the WNTRSimulator to compute headloss in pipes

return_map: bool (optional, default = False)
    Return a skeletonization map.   The map is a dictionary 
    that includes original nodes as keys and a list of skeletonized nodes 
    that were merged into each original node as values.

Returns
--------
A skeletonized WaterNetworkModel object and (if return_map = True) a skeletonization map.
"""
skel = _Skeletonize(wn, use_epanet)

skel.run(pipe_diameter_threshold, branch_trim, series_pipe_merge, 
         parallel_pipe_merge, max_cycles)

if return_map:
    return skel.wn, skel.skeleton_map
else:
    return skel.wn

class _Skeletonize(object):

def __init__(self, wn, use_epanet):

    # Get a copy of the WaterNetworkModel
    self.wn = _deepcopy_wn(wn)

    # Get the WaterNetworkModel graph
    G = self.wn.get_graph()
    G = G.to_undirected()
    self.G = G

    # Create a map of original nodes to skeletonized nodes
    skel_map = {}
    for node_name in self.wn.node_name_list:
        skel_map[node_name] = [node_name]
    self.skeleton_map = skel_map

    # Get a list of junction and pipe names that are associated with controls
    junc_with_controls = []
    pipe_with_controls = []
    for name, control in self.wn.controls():
        for req in control.requires():
            if isinstance(req, Junction):
                junc_with_controls.append(req.name)
            elif isinstance(req, Pipe):
                pipe_with_controls.append(req.name)
    self.junc_with_controls = list(set(junc_with_controls))
    self.pipe_with_controls = list(set(pipe_with_controls))

    # Calculate pipe headloss using a single period EPANET simulation
    duration = self.wn.options.time.duration
    if use_epanet:
        sim = EpanetSimulator(self.wn)
    else:
        sim = WNTRSimulator(self.wn)
    self.wn.options.time.duration = 0
    results = sim.run_sim()
    head = results.node['head']
    headloss = {}
    for link_name, link in self.wn.links():
        headloss[link_name] = float(abs(head[link.start_node_name] - head[link.end_node_name]))
    self.headloss = headloss
    self.wn.options.time.duration = duration

    self.num_branch_trim = 0
    self.num_series_merge = 0
    self.num_parallel_merge = 0

def run(self, pipe_threshold, branch_trim=True, series_pipe_merge=True, 
            parallel_pipe_merge=True, max_cycles=None):
    """
    Run iterative branch trim, series pipe merge, and parallel pipe merge 
    operations based on a pipe diameter threshold.  
    """
    num_junctions = self.wn.num_junctions
    iteration = 0
    flag = True

    while flag:
        if branch_trim:
            self.branch_trim(pipe_threshold)
        if series_pipe_merge:
            self.series_pipe_merge(pipe_threshold)
        if parallel_pipe_merge:
            self.parallel_pipe_merge(pipe_threshold)

        iteration = iteration + 1

        if (max_cycles is not None) and (iteration > max_cycles):
            flag = False
        if num_junctions == self.wn.num_junctions:
            flag = False
        else:
            num_junctions = self.wn.num_junctions

    return self.wn, self.skeleton_map

def branch_trim(self, pipe_threshold):
    """
    Run a single branch trim operation based on a pipe diameter threshold.
    Branch trimming removes dead-end pipes smaller than the pipe 
    diameter threshold and redistributes demands (and associated demand 
    patterns) to the neighboring junction.
    """
    for junc_name in self.wn.junction_name_list:
        if junc_name in self.junc_with_controls:
            continue
        neighbors = list(nx.neighbors(self.G,junc_name))
        if len(neighbors) > 1:
            continue
        if len(neighbors) == 0:
            continue
        neigh_junc_name = neighbors[0] # only one neighbor
        nPipes = len(self.G.adj[junc_name][neigh_junc_name])
        if nPipes > 1:
            continue
        neigh_junc = self.wn.get_node(neigh_junc_name)
        if not (isinstance(neigh_junc, Junction)):
            continue
        pipe_name = list(self.G.adj[junc_name][neigh_junc_name].keys())[0] # only one pipe
        pipe = self.wn.get_link(pipe_name)
        if not ((isinstance(pipe, Pipe)) and \
            (pipe.diameter <= pipe_threshold) and \
            pipe_name not in self.pipe_with_controls):
            continue

        logger.info('Branch trim: '+ str(junc_name) + str(neighbors))

        # Update skeleton map        
        self.skeleton_map[neigh_junc_name].extend(self.skeleton_map[junc_name])
        self.skeleton_map[junc_name] = []

        # Move demand
        junc = self.wn.get_node(junc_name)
        for demand in junc.demand_timeseries_list:
            neigh_junc.demand_timeseries_list.append(demand)
        junc.demand_timeseries_list.clear()

        # Remove node and links from wn and G
        self.wn.remove_link(pipe_name)
        self.wn.remove_node(junc_name)
        self.G.remove_node(junc_name)

        self.num_branch_trim +=1

    return self.wn, self.skeleton_map

def series_pipe_merge(self, pipe_threshold):
    """
    Run a single series pipe merge operation based on a pipe diameter 
    threshold.  This operation combines pipes in series if both pipes are 
    smaller than the pipe diameter threshold. The larger diameter pipe is 
    retained, demands (and associated demand patterns) are redistributed 
    to the nearest junction.
    """
    for junc_name in self.wn.junction_name_list:
        if junc_name in self.junc_with_controls:
            continue
        neighbors = list(nx.neighbors(self.G,junc_name))
        if not (len(neighbors) == 2):
            continue
        neigh_junc_name0 = neighbors[0]
        neigh_junc_name1 = neighbors[1]
        neigh_junc0 = self.wn.get_node(neigh_junc_name0)
        neigh_junc1 = self.wn.get_node(neigh_junc_name1)
        if not ((isinstance(neigh_junc0, Junction)) or \
           (isinstance(neigh_junc1, Junction))):
            continue
        pipe_name0 = list(self.G.adj[junc_name][neigh_junc_name0].keys())
        pipe_name1 = list(self.G.adj[junc_name][neigh_junc_name1].keys())
        if (len(pipe_name0) > 1) or (len(pipe_name1) > 1):
            continue
        pipe_name0 = pipe_name0[0] # only one pipe
        pipe_name1 = pipe_name1[0] # only one pipe
        pipe0 = self.wn.get_link(pipe_name0)
        pipe1 = self.wn.get_link(pipe_name1)
        if not ((isinstance(pipe0, Pipe)) and \
            (isinstance(pipe1, Pipe)) and \
            ((pipe0.diameter <= pipe_threshold) and \
            (pipe1.diameter <= pipe_threshold)) and \
            pipe_name0 not in self.pipe_with_controls and \
            pipe_name1 not in self.pipe_with_controls):
            continue
        # Find closest neighbor junction
        if (isinstance(neigh_junc0, Junction)) and \
           (isinstance(neigh_junc1, Junction)):
            if pipe0.length < pipe1.length:
                closest_junc = neigh_junc0
            else:
                closest_junc = neigh_junc1
        elif (isinstance(neigh_junc0, Junction)):
            closest_junc = neigh_junc0
        elif (isinstance(neigh_junc1, Junction)):
            closest_junc = neigh_junc1
        else:
            continue

        logger.info('Series pipe merge: ' + str(junc_name) + str(neighbors))

        # Update skeleton map    
        self.skeleton_map[closest_junc.name].extend(self.skeleton_map[junc_name])
        self.skeleton_map[junc_name] = []

        # Move demand
        junc = self.wn.get_node(junc_name)
        for demand in junc.demand_timeseries_list:
            closest_junc.demand_timeseries_list.append(demand)
        junc.demand_timeseries_list.clear()

        # Remove node and links from wn and G
        self.wn.remove_link(pipe_name0)
        self.wn.remove_link(pipe_name1)
        self.wn.remove_node(junc_name)
        self.G.remove_node(junc_name)

        # Compute new pipe properties
        props = self._series_merge_properties(pipe0, pipe1)

        # Add new pipe to wn and G
        dominant_pipe = self._select_dominant_pipe(pipe0, pipe1)
        self.wn.add_pipe(dominant_pipe.name, 
                         start_node_name=neigh_junc_name0, 
                         end_node_name=neigh_junc_name1, 
                         length=props['length'], 
                         diameter=props['diameter'], 
                         roughness=props['roughness'], 
                         minor_loss=props['minorloss'],
                         status=props['status']) 
        self.G.add_edge(neigh_junc_name0, 
                        neigh_junc_name1, 
                        dominant_pipe.name)

        self.num_series_merge +=1

    return self.wn, self.skeleton_map

def parallel_pipe_merge(self, pipe_threshold):
    """
    Run a single parallel pipe merge operation based on a pipe diameter 
    threshold.  This operation combines pipes in parallel if both pipes are 
    smaller than the pipe diameter threshold. The larger diameter pipe is 
    retained.
    """

    for junc_name in self.wn.junction_name_list:
        if junc_name in self.junc_with_controls:
            continue
        neighbors = nx.neighbors(self.G,junc_name)
        for neighbor in neighbors:
            parallel_pipe_names = list(self.G.adj[junc_name][neighbor].keys())
            if len(parallel_pipe_names) == 1:
                continue
            for (pipe_name0, pipe_name1) in itertools.combinations(parallel_pipe_names, 2):
                try:
                    pipe0 = self.wn.get_link(pipe_name0)
                    pipe1 = self.wn.get_link(pipe_name1)
                except:
                    continue # one of the pipes removed in previous loop
                if not ((isinstance(pipe0, Pipe)) and \
                   (isinstance(pipe1, Pipe)) and \
                    ((pipe0.diameter <= pipe_threshold) and \
                    (pipe1.diameter <= pipe_threshold)) and \
                    pipe_name0 not in self.pipe_with_controls and \
                    pipe_name1 not in self.pipe_with_controls):
                    continue

                logger.info('Parallel pipe merge: '+ str(junc_name) + str((pipe_name0, pipe_name1)))

                # Remove links from wn and G   
                self.wn.remove_link(pipe_name0)
                self.wn.remove_link(pipe_name1)
                self.G.remove_edge(neighbor, junc_name, pipe_name0) 
                self.G.remove_edge(junc_name, neighbor, pipe_name1)

                # Compute new pipe properties
                props = self._parallel_merge_properties(pipe0, pipe1)

                # Add a new pipe to wn and G
                dominant_pipe = self._select_dominant_pipe(pipe0, pipe1)
                self.wn.add_pipe(dominant_pipe.name, 
                                 start_node_name=dominant_pipe.start_node_name, 
                                 end_node_name=dominant_pipe.end_node_name,
                                 length=props['length'], 
                                 diameter=props['diameter'], 
                                 roughness=props['roughness'], 
                                 minor_loss=props['minorloss'],
                                 status=props['status']) 
                self.G.add_edge(dominant_pipe.start_node_name, 
                                dominant_pipe.end_node_name, 
                                dominant_pipe.name)

                self.num_parallel_merge +=1

    return self.wn, self.skeleton_map

def _select_dominant_pipe(self, pipe0, pipe1):

    # Dominant pipe = larger diameter
    if pipe0.diameter >= pipe1.diameter:
        dominant_pipe = pipe0
    else:
        dominant_pipe = pipe1

    return dominant_pipe

def _series_merge_properties(self, pipe0, pipe1):

    props = {}
    dominant_pipe = self._select_dominant_pipe(pipe0, pipe1)

    props['length'] = pipe0.length + pipe1.length
    props['diameter'] = dominant_pipe.diameter
    props['minorloss'] = dominant_pipe.minor_loss
    props['status'] = dominant_pipe.status

    props['roughness'] = (props['length']/(props['diameter']**4.87))**0.54 * \
        ((pipe0.length/((pipe0.diameter**4.87)*(pipe0.roughness**1.85))) + \
         (pipe1.length/((pipe1.diameter**4.87)*(pipe1.roughness**1.85))))**-0.54

    return props

def _parallel_merge_properties(self, pipe0, pipe1):

    props = {}
    dominant_pipe = self._select_dominant_pipe(pipe0, pipe1)

    props['length'] = dominant_pipe.length
    props['diameter'] = dominant_pipe.diameter
    props['minorloss'] = dominant_pipe.minor_loss
    props['status'] = dominant_pipe.status

    props['roughness'] = ((props['length']**0.54)/(props['diameter']**2.63)) * \
        ((pipe0.roughness*(pipe0.diameter**2.63))/(pipe0.length**0.54) + \
         (pipe1.roughness*(pipe1.diameter**2.63))/(pipe1.length**0.54))

    return props

wn_model = skeletonize(wn, 16, use_epanet = True) wn_model.write_inpfile('Example2_outwithRB')

kaklise commented 4 years ago

We preserve the original base demand values and demand patterns when combining junctions in the skeletonization process. We should clarify this in the user manual. You could aggregate these into new base demands and demand patterns for each junction, but this could potentially create many new patterns, and makes debugging difficult.

bahadurr commented 4 years ago

Katherine

I think I know the issue. The small script you have sent to me converts pipe diameter from INCHES to METERS. I was using simple inch as there was no mention of the conversion in the manual. It is working fine and giving me reasonable results. I will now compare the original model with the skeletonized model

skel_wn = wntr.morph.skeletonize(orig_wn, 16*0.0254)

Thanks again for all your help

Rakesh Bahadur, Ph.D., PMP | Leidos MS FS 1-5-1 11951 Freedom Drive Reston, VA 20190 phone: 571.526.7120 bahadurr@leidos.commailto:bahadurr@leidos.com | leidos.com

[cid:image001.png@01CF9B6D.2D0979D0] This email and any attachments to it are intended only for the identified recipients. It may contain proprietary or otherwise legally protected information of Leidos. Any unauthorized use or disclosure of this communication is strictly prohibited. If you have received this communication in error, please notify the sender and delete or otherwise destroy the email and all attachments immediately.

From: Katherine Klise notifications@github.com Sent: Thursday, March 5, 2020 10:59 AM To: USEPA/WNTR WNTR@noreply.github.com Cc: Bahadur, Rakesh [US-US] RAKESH.BAHADUR@leidos.com; Author author@noreply.github.com Subject: EXTERNAL: Re: [USEPA/WNTR] Skeletonization (#128)

I don't see the attached file, but I can provide an example using Net6 to start the discussion.

The issue might be explained by how demands are extracted from the original and skeletonized model. Within the skeletonization process, base demands and demand patterns are moved to remaining nodes. We can compute total demand at each junction and time using base demands, demand patterns, and the demand multiplier to make sure demands are conserved (using wntr.metrics.expected_demand). If we just extract the base demand (using query_node_attribute), the results will be different because this does not account for the multiple demands and patterns on each junction.

In the following example, Net 6 is skeletonized using a 16 inch pipe diameter threshold. The total expected demand is conserved (130.42 m3/s) but the base demands differ. Does this explain the issue you are seeing?

import wntr

Create a water network model

orig_wn = wntr.network.WaterNetworkModel('networks/Net6.inp')

print(orig_wn.describe())

Skeletonize the model using a 16 inch pipe diameter threshold

skel_wn = wntr.morph.skeletonize(orig_wn, 16*0.0254)

print(skel_wn.describe())

Extract the expected demand for each water network model

This uses base demand, demand pattern, and the demand multiplier

Expected demand for both models should be the same

orig_expected_demand = wntr.metrics.expected_demand(orig_wn)

skel_expected_demand = wntr.metrics.expected_demand(skel_wn)

orig_expected_demand.sum(axis=1).plot() # plot timeseries of total demand

skel_expected_demand.sum(axis=1).plot()

print(orig_expected_demand.sum().sum()) # print total demand

print(skel_expected_demand.sum().sum())

Extract the base demand for each water network model

The base demand will not be the same

orig_base_demand = orig_wn.query_node_attribute('base_demand')

skel_base_demand = skel_wn.query_node_attribute('base_demand')

print(orig_base_demand.sum()) # print total base demand

print(skel_base_demand.sum())

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/USEPA/WNTR/issues/128?email_source=notifications&email_token=ACNPI4GVDJEJRYXPR4SXINDRF7D2TA5CNFSM4LCLXZP2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN52CRY#issuecomment-595304775, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACNPI4CC7INRRK7HEDLYOQTRF7D2TANCNFSM4LCLXZPQ.

kaklise commented 4 years ago

I'm glad to hear that it's working. WNTR uses SI units throughout the software. https://wntr.readthedocs.io/en/latest/units.html

bahadurr commented 4 years ago

Katherine

I am getting the following Waning:

UserWarning:

Not all curves were used in "C:\WNTR\Testing_Script/Texas_Example2_Testing_FOUO.inp"; added with type None, units conversion left to user

What does it mean?

Thanks again for all your help today

Rakesh Bahadur, Ph.D., PMP | Leidos MS FS 1-5-1 11951 Freedom Drive Reston, VA 20190 phone: 571.526.7120 bahadurr@leidos.commailto:bahadurr@leidos.com | leidos.com

[cid:image001.png@01CF9B6D.2D0979D0] This email and any attachments to it are intended only for the identified recipients. It may contain proprietary or otherwise legally protected information of Leidos. Any unauthorized use or disclosure of this communication is strictly prohibited. If you have received this communication in error, please notify the sender and delete or otherwise destroy the email and all attachments immediately.

From: Katherine Klise notifications@github.com Sent: Thursday, March 5, 2020 12:39 PM To: USEPA/WNTR WNTR@noreply.github.com Cc: Bahadur, Rakesh [US-US] RAKESH.BAHADUR@leidos.com; Author author@noreply.github.com Subject: EXTERNAL: Re: [USEPA/WNTR] Skeletonization (#128)

I'm glad to hear that it's working. WNTR uses SI units throughout the software. https://wntr.readthedocs.io/en/latest/units.html

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/USEPA/WNTR/issues/128?email_source=notifications&email_token=ACNPI4E3XSB45GDZC2UGITDRF7PSPA5CNFSM4LCLXZP2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN6GO6Q#issuecomment-595355514, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACNPI4F34PJU6XNXWV6Z5YTRF7PSPANCNFSM4LCLXZPQ.

dbhart commented 4 years ago

Hi Rakesh,

This warning means that for things like energy curves, we don’t use them in WNTR so we leave them alone. We assign curves their properties (pump, energy, etc.) based on what elements they get assigned to (pump-head, tank-volumes). So if there is a curve that never gets used, or which is an energy curve, we put this warning in there so the user knows that there may be a curve we didn’t identify the units for, so we left them alone. Hopefully that makes sense.

Dave H.