OPAL++ Example - 2d axisymmetric Tesla Pump

This is an example of an Optimization performed with OPAL++.

An axisymmetric 2D representation of a Tesla Pump has been investigated.

It was run on Karman (2014).

Optimization Details

The objective of the optimization was to increase the hydraulic efficiency of a simplified tesla pump. The pump was modeled 2d-axisymmetric. The mesh was parametrically created in ICEM CFD and the simulations have been run in Starccm. Because the highest efficiency in a specific operating point was searched for, each simulation had to find its own optimal rotation speed to pump 5L/min with a pressure rise of 100mmHg.

Simulation Flow

The simulation script for opal (here: Tesla_Simulation.o2script) does the following things.

  1. Substitute the defined parameter into the file in_parameter.tcl
  2. Borrow an ANSYS License for ICEM CFD
    1. remove previous indicator (version.txt)
    2. run script which writes version.txt when license is ready (et_icem_license.sh)
    3. wait for version.txt
    4. retain version to the output variables
  3. run ICEM CFD with available license version
  4. retain output of ICEM CFD and check for errors
  5. run in_ioHelper.tcl to prepare StarCCM+ macro, and retain output (can be done completely in OPAL++, too. Used tcl script because i was more familiar with it.)
  6. run StarCCM+ simulation with in_tesrot.java macro and retain output
  7. check for errors during simulation
  8. postprocess the simulation results with in_ioHelper
  9. retain output values
  10. save simulation in zip file and retain it
  11. Check for validity (whether there are any errors in the simulation), mark valid if criteria met

OPAL++ Input Files

  • Tesla_Optimization.o2script
  • Tesla_Simulation.o2script
    • in_icem.tcl
    • in_ioHelper.tcl
    • in_parameter.tcl
    • in_tesrot.java
    • in_extractvalues.sh
    • get_icem_license.sh

Results

How to run

  • Create Checkstring
~/OPAL/OPAL -opal.mode=NUMERIC_TEST -opal.file={Tesla_Optimization.o2script}
  • submit job with pbs script

Files

Tesla_Optimization.o2script

Tesla_Optimization.o2script
#include "stdafx.h"
#include "OPAL2_MasterScript_V22.h"
 
OPAL2_MASTERSCRIPT_BEGIN
 
  OPAL2.SETUP.PROBLEM.NAME(STR/*Tesla Pump 2d Axisymmetric Model Main Case Optimization*/);
	OPAL2.SETUP.PROBLEM.DESCRIPTION(STR/*This module is responsible for the Tesla pump optimization*/);
 
	// Setup Scripts
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Input/Tesla_simulation.o2script*/,STR/*simulation.o2script*/,OPAL2SCRIPT);
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Results_AxiSym_Main3_Detail2*/,STR/**/,WORKDIRECTORY);
 
	// Input scripts
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Input/get_icem_license.sh*/,STR/*get_icem_license.sh*/,TEXT_DATAFILE);
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Input/in_icem.tcl*/,STR/*in_icem.tcl*/,TEXT_DATAFILE);	
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Input/in_parameter.tcl*/,STR/*in_parameter.tcl*/,TEXT_DATAFILE);
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Input/in_tesrot.java*/,STR/*in_tesrot.java*/,TEXT_DATAFILE);
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Input/in_ioHelper.tcl*/,STR/*in_ioHelper.tcl*/,TEXT_DATAFILE);	
	OPAL2.SETUP.FILES.ADD(STR/*/home/engel/Tesla/AxiSym/Optimization/Input/in_extractvalues.sh*/,STR/*in_extractvalues.sh*/,TEXT_DATAFILE);
 
	// Files to Retrieve
	OPAL2.SETUP.FILES.ADD(STR/*case.zip*/,STR/*case.zip*/,RETAINFILE);	
 
	//Solver method
	OPAL2.SETUP.SOLVER.METHOD=GENETIC1;
 
	//Variables
	OPAL2.SETUP.PROBLEM.VARIABLES.ADD_REAL(STR/*RDOUT*/,45.0,50.0,1.0);	
	OPAL2.SETUP.PROBLEM.VARIABLES.ADD_REAL(STR/*RDIN*/,7.0,8.0,1.0);	
	OPAL2.SETUP.PROBLEM.VARIABLES.ADD_REAL(STR/*GAP*/,0.75,1.0,1.0);	
	OPAL2.SETUP.PROBLEM.VARIABLES.ADD_REAL(STR/*WALLSPACE*/,1.2,1.7,1.0);	
 
 
	OPAL2.SETUP.PROBLEM.VARIABLES.FINALIZED=TRUE;
 
	//Derived variables
	OPAL2.SETUP.PROBLEM.DERIVED_VARIABLES.ADD(STR/*ICEM_VERSION*/,STR/*0*/);	
	OPAL2.SETUP.PROBLEM.DERIVED_VARIABLES.ADD(STR/*REPORT_SPEED*/,STR/*0*/);	
	OPAL2.SETUP.PROBLEM.DERIVED_VARIABLES.ADD(STR/*REPORT_PRESDROP*/,STR/*0*/);
	OPAL2.SETUP.PROBLEM.DERIVED_VARIABLES.ADD(STR/*REPORT_TORQUE*/,STR/*0*/);
	OPAL2.SETUP.PROBLEM.DERIVED_VARIABLES.ADD(STR/*REPORT_EFFICIENCY*/,STR/*0*/);
	OPAL2.SETUP.PROBLEM.DERIVED_VARIABLES.FINALIZED=TRUE;
 
	//Maximize power coefficient (not direct evaluation)
	OPAL2.SETUP.PROBLEM.OBJECTIVES.ADD_INTERNAL(STR/*OBJECTIVE_CP*/,MAXIMIZE,STR/*REPORT_EFFICIENCY*/,false,false,0,1,NO_NORMALIZATION,1);
	OPAL2.SETUP.PROBLEM.OBJECTIVES.FINALIZED=TRUE;
 
        OPAL2.SETUP.PROBLEM.CONSTRAINS.ADD(STR/*RDIN+5.0-RDOUT*/,DESIGN,1.0);
	OPAL2.SETUP.PROBLEM.CONSTRAINS.FINALIZED=TRUE;
	OPAL2.SETUP.PROBLEM.FINALIZED=TRUE;
 
	//Separator: 'space'
	OPAL2.RUNTIME.DATAFILES.SEPARATOR=STR/* */;
	//Print information of all individuals
	OPAL2.RUNTIME.INFORMATION.REFRESHRATE=INDIVIDUAL;
	//Solver settings
	OPAL2.SETUP SOLVER.CHANGE_SETTINGS(N,100);
	OPAL2.SETUP SOLVER.CHANGE_SETTINGS(GENERATION_MAX,3)	
 	//Dual-mode license
        OPAL2.RUNTIME.LICENSE_SHIFTING(3.0);
        OPAL2.RUNTIME.LICENSING.ADD({ANSYS},1,1);
    	//DOE
    	OPAL2.SETUP.SOLVER.INITIALIZATION.TYPE=SOBOL;
    	OPAL2.SETUP.SOLVER.INITIALIZATION.CONSTRAINED=TRUE;
 
OPAL2_MASTERSCRIPT_END

Tesla_Simulation.o2script

Tesla_Simulation.o2script
#include "stdafx.h"
#include "OPAL2_SimulationScript_V22.h"
 
OPAL2_SIMULATIONSCRIPT_BEGIN
	//Substitute values into the files
	OPAL2S.TOOLS.SUBSTITUTE(STR/*in_parameter.tcl*/);
 
 
	// -------------------------------------------------
	// Make Mesh
	// -------------------------------------------------
	OPAL2S.LICENSING.BORROW({ANSYS});	
		OPAL2S.TOOLS.COMMAND(STR/*rm version.txt*/);
		OPAL2S.TOOLS.COMMAND(STR/*sh get_icem_license.sh &*/);
		OPAL2S.TOOLS.WAITFOR(STR/*version.txt*/);
		OPAL2S.PROCESSION.PROCESSFILE(STR/*version.txt*/,SINGLE_VALUE,STR/*ICEM_VERSION*/);
	OPAL2S.LICENSING.RETURN({ANSYS});
 
	OPAL2S.CONTROL.IF(STR/*ICEM_VERSION>2.5 && ICEM_VERSION<3.5*/);		
		OPAL2S.TOOLS.COMMAND(STR/*export LM_LICENSE_FILE=1055@liclux.urz.uni-magdeburg.de;cp /home/engel/.ansys/v140/licensing/license.preferences.xml.default /home/engel/.ansys/v140/licensing/license.preferences.xml;/home/software/ansys_inc/v140/icemcfd/linux64_amd/bin/icemcfd -batch -script in_icem.tcl > icemlog.txt 2>&1;cp /home/engel/.ansys/v140/licensing/license.preferences.xml.default /home/engel/.ansys/v140/licensing/license.preferences.xml*/,TRUE);
	OPAL2S.CONTROL.END;				
 
	OPAL2S.CONTROL.IF(STR/*ICEM_VERSION>3.5 && ICEM_VERSION<4.5*/);		
		OPAL2S.TOOLS.COMMAND(STR/*export LM_LICENSE_FILE=1055@liclux.urz.uni-magdeburg.de;cp /home/engel/.ansys/v140/licensing/license.preferences.xml.default /home/engel/.ansys/v140/licensing/license.preferences.xml;/home/software/ansys_inc/v140/icemcfd/linux64_amd/bin/icemcfd -batch -script in_icem.tcl > icemlog.txt 2>&1;cp /home/engel/.ansys/v140/licensing/license.preferences.xml.default /home/engel/.ansys/v140/licensing/license.preferences.xml*/,TRUE);
	OPAL2S.CONTROL.END;	
 
	OPAL2S.TOOLS.APPEND_TO_LOG(STR/*icemlog.txt*/,STR/*ICEM output:*/);	
 
	// Error detection	
	OPAL2S.PROCESSION.VALIDITYCHECK(STR/*icemlog.txt*/,STR/*Licensed number of users already reached*/,NOT_PRESENT,STR/*ICEM, no licences left*/);
	OPAL2S.PROCESSION.VALIDITYCHECK(STR/*icemlog.txt*/,STR/*ERROR*/,NOT_PRESENT,STR/*ICEM meshing error*/);
	OPAL2S.PROCESSION.VALIDITYCHECK(STR/*icemlog.txt*/,STR/*Error*/,NOT_PRESENT,STR/*ICEM meshing error*/);
 
 
 
	// -------------------------------------------------
	//  RUN
	// -------------------------------------------------
	OPAL2S.TOOLS.COMMAND(STR/*tclsh in_ioHelper.tcl start > ioHelper.txt 2>&1*/,TRUE);		
	OPAL2S.TOOLS.APPEND_TO_LOG(STR/*ioHelper.txt*/,STR/*ioHelper output:*/);		
 
	OPAL2S.TOOLS.COMMAND(STR/*/home/software/CD-adapco/STAR-CCM+9.02.005-R8/star/bin/starccm+ -new -batch in_tesrot.java -licpath 1999@flex.cd-adapco.com -podkey xKf6MtxYPARpfvJDUP9+SA -rsh ssh -np 8 -power > output.txt 2>&1*/,TRUE);
	OPAL2S.TOOLS.APPEND_TO_LOG(STR/*output.txt*/,STR/*Star CCM+ output:*/);
 
	// Error detection
	OPAL2S.PROCESSION.VALIDITYCHECK(STR/*output.txt*/,STR/*Divergence detected*/,NOT_PRESENT,STR/*Divergence error*/);
	OPAL2S.PROCESSION.VALIDITYCHECK(STR/*output.txt*/,STR/*point exception*/,NOT_PRESENT,STR/*SIGFPE detected*/);
	OPAL2S.PROCESSION.VALIDITYCHECK(STR/*output.txt*/,STR/*Error*/,NOT_PRESENT,STR/*Error detected*/);
 
	OPAL2S.TOOLS.COMMAND(STR/*tclsh in_ioHelper.tcl end > ioHelper.txt 2>&1*/,TRUE);		
	OPAL2S.TOOLS.APPEND_TO_LOG(STR/*ioHelper.txt*/,STR/*ioHelper output:*/);		
 
 
	// -------------------------------------------------
	// Post
	// -------------------------------------------------	
 
	OPAL2S.PROCESSION.PROCESSFILE(STR/*speed2.val*/,SINGLE_VALUE,STR/*REPORT_SPEED*/);		
	OPAL2S.PROCESSION.PROCESSFILE(STR/*presDrop.val*/,SINGLE_VALUE,STR/*REPORT_PRESDROP*/);	
	OPAL2S.PROCESSION.PROCESSFILE(STR/*torque.val*/,SINGLE_VALUE,STR/*REPORT_TORQUE*/);		
	OPAL2S.PROCESSION.PROCESSFILE(STR/*eff.val*/,SINGLE_VALUE,STR/*REPORT_EFFICIENCY*/);		
 
	OPAL2S.TOOLS.COMMAND(STR/*zip -r9 case.zip . > ziplog.txt 2>&1*/,TRUE);			
	OPAL2S.TOOLS.APPEND_TO_LOG(STR/*ziplog.txt*/,STR/*zip output:*/);
 
        // Error detection
        OPAL2S.PROCESSION.VALIDITYCHECK(STR/*output.txt*/,STR/*Divergence detected*/,NOT_PRESENT,STR/*Divergence error*/);
        OPAL2S.PROCESSION.VALIDITYCHECK(STR/*output.txt*/,STR/*point exception*/,NOT_PRESENT,STR/*SIGFPE detected*/);
        OPAL2S.PROCESSION.VALIDITYCHECK(STR/*output.txt*/,STR/*Error*/,NOT_PRESENT,STR/*Error detected*/)
        OPAL2S.PROCESSION.VALIDITYCHECK(STR/*output.txt*/,STR/*Stopping criteria X-momentum Criterion and Y-momentum Criterion and Z-momentum Criterion and Continuity Criterion and Tke Criterion and Sdr Criterion satisfied*/,PRESENT,STR/*Not Converged*/);
 
	//Check, if any dynamic library based contrains have been violated
	OPAL2S.NATIVE.CONTROL.IF(STR/*constrain>0*/);
		OPAL2S.NATIVE.PROCESSION.EQUATION(STR/*validity*/,STR/*1*/);
	OPAL2S.NATIVE.CONTROL.END;
 
OPAL2_SIMULATIONSCRIPT_END

get_icem_license.sh

get_icem_license.sh
#!/bin/bash
CENTRAL_ALL=30
CENTRAL_BORROWED=0
CENTRAL_PROTECTED=0
CENTRAL_AVAILABLE=0
TEACHING_ALL=25
TEACHING_AVAILABLE=0
TEACHING_BORROWED=0
TEACHING_PROTECTED=0
SAFETY_MARGIN=1
SUCCESS=0
 
sleep 4
 
export LM_LICENSE_FILE=1055@liclux.urz.uni-magdeburg.de
 
while [ $SUCCESS -ne 1 ]; do
        #
	TEACHING_PROTECTED=2
	CENTRAL_PROTECTED=0
	CENTRAL_BORROWED=$(/home/software/ansys_inc/v140/fluent/bin/lmstat -f aa_r_cfd | grep "Total of" | awk '{print $11}')
	TEACHING_BORROWED=$(/home/software/ansys_inc/v140/fluent/bin/lmstat -f aa_t_a | grep "Total of" | awk '{print $11}')
	let CENTRAL_AVAILABLE=$CENTRAL_ALL-$CENTRAL_BORROWED
	let TEACHING_AVAILABLE=$TEACHING_ALL-$TEACHING_BORROWED
 
	#Preference List: Parallel (1); Lehrstuhl (2); Teaching (3); CFD (4)
 
 
if [ $TEACHING_AVAILABLE -gt $TEACHING_PROTECTED ]; then
	echo 3 > version.txt
				SUCCESS=1
	else
 
		if [ $CENTRAL_AVAILABLE -gt $CENTRAL_PROTECTED ]; then
							echo 4 > version.txt
							SUCCESS=1
			else
							sleep 4
			fi
 
 
fi
 
done

in_ioHelper.tcl

in_ioHelper.tcl
## Tcl-Script
## Part of Work: Optimizing the hydraulic efficiency of a Tesla Pump
## Rotor Independent Optimization
## Calculate new speed to match pressure head

# -----------------------------------------------
# Input Parameters
# -----------------------------------------------
	# Procedure argument
	set round [lindex $argv 0]
 
	source in_parameter.tcl
 
	set pi [expr acos(-1)]
	
# -----------------------------------------------
# Helper Procedures
# -----------------------------------------------
proc readSingleValueFile {filename} {
	set fileID [open $filename]
	while {1} {
			set line [gets $fileID]
			if {[eof $fileID]} {
					close $fileID
					break
			}
			return $line
	}
}
 
proc writeSingleValueFile {filename value} {
	set fileID [open $filename "w"]
	puts $fileID $value
	close $fileID
}
 
proc replaceFileString {filename pattern replacement} {
set tmpfile ${filename}.new
set in [open $filename r+]; # to force error if file is not writable
set out [open $tmpfile w]
while {[gets $in line]>=0} {
regsub $pattern $line $replacement line ;# Replace in place
puts $out $line
}
close $in
close $out
file rename -force $tmpfile $filename
} 

	# if { [catch {} fid]} {
		# puts stderr "\n$fid"
		# set error 1
	# }

# -----------------------------------------------
# Main Part
# -----------------------------------------------
 
if {$round == "start"} {
	# Approximate new speed (Ekman Zahl = 1)
	set speed1 [expr $my/$rho*($ekman/($dsksp/2000))*($ekman/($dsksp/2000))]
 
	if {$speed1 > 500.0}	{	
	set speed1 500.0
	}
	
	# Write Files	
	writeSingleValueFile speed1.val $speed1
 
	replaceFileString in_tesrot.java RHOKGKM $rho
	replaceFileString in_tesrot.java KINVISC $my
	replaceFileString in_tesrot.java MASSFLOW $mflow	
	replaceFileString in_tesrot.java RPMVALUE $speed1			
 
}
 
if {$round == "end"} {
	set error 0
	if { [catch {	exec sh in_extractvalues.sh presDrop} fid]} {
		puts stderr "in_extractvalues.sh presDrop\n$fid"
		set error 1		
	}
	if { [catch {	exec sh in_extractvalues.sh torque} fid]} {
		puts stderr "in_extractvalues.sh torque\n$fid"
		set error 1		
	}
	if { [catch {	exec sh in_extractvalues.sh speed} fid]} {
		puts stderr "in_extractvalues.sh speed\n$fid"
		set error 1		
	}	
 
	if { [catch {	set speed [readSingleValueFile speed2.val]} fid]} {
		puts stderr "readSingleValueFile speed2.val\n$fid"
		set error 1	
	}	
	if { [catch {	set presDrop [readSingleValueFile presDrop.val]} fid]} {
		puts stderr "readSingleValueFile presDrop.val\n$fid"
		set error 1		
	}	
	if { [catch {	set torque [readSingleValueFile torque.val]	} fid]} {
		puts stderr "readSingleValueFile torque.val\n$fid"
		set error 1		
	}
 
	if {$error} {
		writeSingleValueFile eff.val 0		
	} else {
	set omega $speed
	set gh [expr ($presDrop)/$rho]
	set power [expr $torque*$omega]
	set efficiency [expr -$mflow*$gh/$power]
 
	if {$efficiency > 1.0} {
		set efficiency -99.0}
 
	writeSingleValueFile eff.val $efficiency		
	}
}

in_icem.tcl

in_icem.tcl
## Tcl-Script for ANSYS ICEM CFD
## Part of Work: Optimizing the hydraulic efficiency of a Tesla Pump
## 
## Creates Geometry and Mesh of a 2D axisymmetrical configuration
## Outputs a ANSYS Fluent Mesh (V6)
## 
## written by Sebastian Engel
## last edit: 18.08.14
## 
##
## todo 
##	some more geometry parameter validity checks
##
##
##
  # Debug settings, should be 0 (0/1)
  set fast 0
	
	# Export name
	set fluentfile "slice"

# -----------------------------------------------
# Input Parameters
# -----------------------------------------------
  # Spacing between disk [mm]
#  set dsksp 0.1 
  # thickness of every disk [mm]
  set dskth 0.25 
  # number of disks
  set dsknm 4
  # Wallspacing [mm]
#  set dskwllsp 0.05
  # Outer disk radius [mm]
#  set dskrdout 60
  # Inner disk radius [mm]
#  set dskrdin 11
  # Inlet pipe radius [mm]
  set inlrd 5.0
  # inlet verrundung [mm]
  set inlvrd 2.0
  # innerer Öffnungswinkel am innerem Radius [°] [0<30]
  set inopang 0.0
 # set outopang 5
  # Diffusor length [mm]
  set difflgh 100
  # Spacing between Outer radius of disk and beginning of Diffusor [mm]
  set diffsp 1.0
  # Inlet length
  set inllgh 30.0
 
source in_parameter.tcl
	
	# BL thickness
	set blth [expr $dsksp/4.]

# -----------------------------------------------
# Mesh Options
# -----------------------------------------------

	# Global mesh factor
	  set nglob 0.7
  # Boundary Layer Nodes
  set nbl 9
  # Disk surface Nodes [n/mm]
  set ndsksrf 10
  # Disk space Nodes [n/mm]
  set ndsksp 25
  # Inlet region Nodes [n/mm]
  set ninlreg 10
  # Diffusor space Nodes [n/mm]
  set ndiffsp 15
  # Diffusor Nodes [n/mm]
  set ndifflgh 7
  # Fillet
  set nfillet 8	
	
  # Inlet length
  set ninllgh 5
	
	# Global factoring
	  # Disk thickness Nodes [n/mm]
		set ndskth [ceil [expr (1+($dskth+2*$blth)/$dsksp)*$ndsksp*$nglob]]
		set ndsksrf [ceil [expr $ndsksrf*$nglob]]
		set ndsksp [ceil [expr $ndsksp*$nglob]]
		set ndiffsp [ceil [expr $ndiffsp*$nglob]]
		set ndifflgh [ceil [expr $ndifflgh*(0.25+0.75*$nglob)]]
		set nfillet [ceil [expr $nfillet*(0.1+0.9*$nglob)]]
		set ninllgh [ceil [expr $ninllgh*(0.25+0.75*$nglob)]]
		set ndsksrf [ceil [expr $ndsksrf*$nglob]]		
		set ninlreg [ceil [expr $ninlreg*(0.2+0.8*$nglob)]]	
		
 
# -----------------------------------------------
# Non-Parameter Geometry Enteties
# -----------------------------------------------

 
 
#  set wdir "C:/CFD/TslASym/MeshScripts"
	set wdir [pwd]
  set staticdir $wdir
 
	ic_chdir $wdir
 
	global env
  set icemenv $env(ICEM_ACN)
  
 
 
 
  # Helper (simplify loop definitions and other operations)
  array set clp {1 2 2 3 3 4 4 1}
  array set dlp {1 5 2 6 3 7 4 8}
  set pi [expr acos(-1)]
  set degtorad [expr $pi/180.]

  # Get coordinates from point
  proc pcrd {point xyz} {
    return [lindex [ic_vset -vec - $point] $xyz]
  }
  # Print output to log, with system time option
 
 
  proc log [list str [list time 1] [list name $fluentfile]] {
      set systemTime [clock seconds]
      if {$time} { ic_mess "[clock format $systemTime -format %H:%M:%S]: $str\n"} else {
		   ic_mess "$str\n"}
 
  }
	
 
 
# Some log describtions	
log "#################################################" 0
log "# Tesla pump Mesh, written by Sebastian Engel   #" 0
log "#################################################" 0
log "Date: [cmd_date]" 0
log "Current User: [cmd_whoami]"
log "System info: [cmd_uname_a]"
log "Free Memory: [cmd_freemem]"
log "ICEM Path: $icemenv"
log "Current Work Directory: $wdir"
if {[ic_batch_mode]} {log "Script executed in batch mode"} else {log "Script executed in Interactive mode"}
log "#################################################" 0	
log "Input Parameter Output"
log "         sym. Disks \[1\]: $dsknm"
log "         Disks gap \[mm\]: $dsksp"
log " Disk outer radius \[mm\]: $dskrdout"
log " Disk inner radius \[mm\]: $dskrdin"
log "    Disk thickness \[mm\]: $dskth"
log "  Disk-Housing gap \[mm\]: $dskwllsp"
log "      Inlet radius \[mm\]: $inlrd"
log "     Fillet radius \[mm\]: $inlvrd"
log "Inner Opening angle \[°\]: $inopang"
log "#################################################" 0	

 
 
 
# -----------------------------------------------
log "Checking some parameters"
# -----------------------------------------------
	if {$dskrdout < $inlrd} { 
		log "Error! Disk smaller than inlet tube."
		break }
	if {$dskrdout < $dskrdin} { 
		log "Error! Inner radius bigger than outer radius"
		break }	
	if {$dskrdout <= $inlvrd} { 
		log "Error! Verrundung groesser als Disk radius."
		break }
	if {$dsknm < 1} { 
		log "Error! Disk number wrong."
		break }
	if {[expr (($dsknm-1)*($dsksp+$dskth))*tan($inopang*$degtorad)] > [expr $dskrdout -$dskrdin]} { 
		log "Error! Opening angle to big. Max Opening angle: [expr atan(($dskrdout -$dskrdin)/(($dsknm-1.)*($dsksp+$dskth)))/$degtorad] deg. Entered angle: $inopang deg."
		break }
	if {$inopang < 0} {
		log "Error! Opening angle smaller than 0"
		break }
	if {$dskrdout <= [expr $inlvrd+$inlrd]} { 
		log "Error! Verrundung und Inlet groesser als Disk radius."
		break }
	if {$inopang >= 30} { 
		log "Warning, possibly bad elements due to high opening angle" }
	if {[expr $inlrd+0.8*$inlvrd] >= $dskrdin} { 
		log "Warning, possibly bad elements due to inlet radius / dskrdin ratio" }	
 
log "Checks completed"

# -----------------------------------------------
log "Start ICEM Session"
# -----------------------------------------------#
	ic_unload_tetin
	ic_hex_unload_blocking 
	ic_unload_mesh 
	ic_empty_tetin
	ic_boco_clear_icons
	ic_csystem_display all 0
	ic_csystem_set_current global
# -----------------------------------------------
log "Geometry creation"
# -----------------------------------------------
	# Points in symmetry plane
	ic_point {} GEOM p.sym.center "0,0,0"
	ic_point {} DISKS p.disk.1.1 "0,[expr $dskrdout],0"
	ic_point {} DISKS p.disk.1.2 "0,[expr $dskrdin],0"
	ic_point {} SYM p.sym.diffsp "0,[expr $dskrdout+$diffsp],0"
	ic_point {} SYM p.sym.diffout "0,[expr $dskrdout+$diffsp+$difflgh],0" 
 
	set dsktipreduc [expr $dskth/3.]
	set dsktiprad [expr (4*$dsktipreduc*$dsktipreduc+$dskth*$dskth)/(8.*$dsktipreduc)]
	
	# Arc centerpoints
	ic_point {} DISKS p.disk.1.a.4	"0,[expr $dskrdout-$dsktiprad],0"
	ic_point {} DISKS p.disk.1.a.2	"0,[expr $dskrdin+$dsktiprad],0"

	# First disk points
	ic_point {} DISKS p.disk.1.4 "[expr $dskth/2.],[expr $dskrdout-$dsktipreduc],0"
	ic_point {} DISKS p.disk.1.3 "[expr $dskth/2.],[expr $dskrdin+$dsktipreduc],0"

	# Connecting symmetry lines
	ic_curve point SYM c.sym.inreg "p.sym.center p.disk.1.2"
	ic_curve point SYM c.sym.diffsp "p.disk.1.1 p.sym.diffsp"
	ic_curve point SYM c.sym.diff "p.sym.diffsp p.sym.diffout"
	
	# Connecting first disk lines
	for {set x 1} {$x <4} {incr x} {
	  ic_curve point DISKS c.disk.1.$x "p.disk.1.$x p.disk.1.$clp($x)"
		incr x
		ic_curve arc_ctr_rad DISKS c.disk.1.$x "p.disk.1.a.$x p.disk.1.$x p.disk.1.$clp($x) $dsktiprad {} {} 0"
	 }
	 
	# Creating remaining disks 
	for {set y 2} {$y<=$dsknm} {incr y} {
		ic_point {} DISKS p.disk.$y.1 "[expr -$dskth/2.+($y-1)*($dsksp+$dskth)] [expr $dskrdout-$dsktipreduc] 0"	
		ic_point {} DISKS p.disk.$y.2 "[expr -$dskth/2.+($y-1)*($dsksp+$dskth)] [expr $dskrdin+tan($inopang*$degtorad)*($dsksp+$dskth)*($y-1)+$dsktipreduc] 0"
		ic_point {} DISKS p.disk.a.$y.2	"[expr ($y-1)*($dsksp+$dskth)],[expr $dskrdin+tan($inopang*$degtorad)*($dsksp+$dskth)*($y-1)+$dsktiprad],0"
		ic_point {} DISKS p.disk.$y.3 "[expr $dskth/2.+($y-1)*($dsksp+$dskth)] [expr $dskrdin+tan($inopang*$degtorad)*($dsksp+$dskth)*($y-1)+$dsktipreduc] 0"
		ic_point {} DISKS p.disk.$y.4 "[expr $dskth/2.+($y-1)*($dsksp+$dskth)] [expr $dskrdout-$dsktipreduc] 0"		
		ic_point {} DISKS p.disk.a.$y.4	"[expr ($y-1)*($dsksp+$dskth)],[expr $dskrdout-$dsktiprad],0"		
		for {set x 1} {$x <4} {incr x} {
			ic_curve point DISKS c.disk.$y.$x "p.disk.$y.$x p.disk.$y.$clp($x)"
			incr x
			ic_curve arc_ctr_rad DISKS c.disk.$y.$x "p.disk.a.$y.$x p.disk.$y.$x p.disk.$y.$clp($x) $dsktiprad {} {} 0"
		}	
	}

 
 
 
	# Creating wall
	ic_point {} WALL p.wll.diffsp "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp] [expr $dskrdout+$diffsp] 0"
	ic_point {} WALL p.wll.out "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp] $dskrdout 0"
	ic_point {} WALL p.wll.outvrd "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp] [expr $inlrd+$inlvrd] 0"
	ic_point {} WALL p.wll.outvrdcenter "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp+$inlvrd] [expr $inlrd+$inlvrd] 0"
	ic_point {} WALL p.wll.inlvrd "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp+$inlvrd] [expr $inlrd] 0"
	ic_point {} WALL p.wll.inl "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp+$inllgh] $inlrd 0"
	ic_point {} INLET p.inl.center "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp+$inllgh] 0 0"

	  # Diffusor
	  for {set x 1} {$x < 11} {incr x} {
	    set temprad [expr [pcrd p.wll.diffsp 1]+$difflgh/10.0*$x]
	    set tempwidth [expr [pcrd p.wll.diffsp 0]*[pcrd p.wll.diffsp 1]/$temprad]
	    if {$x == 0} {
			ic_point {} WALL "p.wll.diffout.$x" "[expr 1.0*$tempwidth] $temprad 0"
			incr x
			set temprad [expr [pcrd p.wll.diffsp 1]+$difflgh/10.0*$x]
			ic_point {} WALL "p.wll.diffout.$x" "[expr 1.0*$tempwidth] $temprad 0"			
			} else {
			ic_point {} WALL "p.wll.diffout.$x" "$tempwidth $temprad 0"
			}
	  }
 
	ic_curve point WALL c.wll.diff "p.wll.diffsp p.wll.diffout.1 p.wll.diffout.2 p.wll.diffout.3 p.wll.diffout.4 p.wll.diffout.5 p.wll.diffout.6 p.wll.diffout.7 p.wll.diffout.8 p.wll.diffout.9 p.wll.diffout.10"
	ic_curve point WALL c.wll.diffsp "p.wll.diffsp p.wll.out"
	ic_curve point WALL c.wll.disk "p.wll.out p.wll.outvrd"
	ic_curve arc_ctr_rad WALL c.wll.vrd "p.wll.outvrdcenter p.wll.outvrd p.wll.inlvrd $inlvrd {} {} 0"
	ic_curve point WALL c.wll.inl "p.wll.inlvrd p.wll.inl"
	ic_curve point OUTLET c.outlet "p.sym.diffout p.wll.diffout.10"
	ic_curve point INLET c.inlet "p.inl.center p.wll.inl"
	ic_curve point AXIS c.axis "p.sym.center p.inl.center"

	# interface
	ic_point {} GEOM p.int.out "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp/2.] $dskrdout 0"
	ic_point {} GEOM p.int.center "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp/2.] 0 0"
# 	ic_curve point INTERIOR c.int.out "p.int.out p.wll.out"
# 	ic_curve point INTERIOR c.int.rd "p.int.center p.int.out"
	ic_point {} GEOM p.int.center.help "[expr ($dsknm-1)*($dsksp+$dskth)+$dskth/2.+$dskwllsp] 0 0"
 
 
 
	ic_curve point POUTLET c.poutlet "p.sym.diffsp p.wll.diffsp"
	ic_point {} PINLET p.inlvrd.center "[pcrd p.wll.inlvrd 0] 0 0"
	ic_curve point PINLET c.pinlet "p.wll.inlvrd p.inlvrd.center"
	
	# Helper Point
	ic_point {} GEOM p.inlvrd.helper.0 "[pcrd p.wll.outvrd 0],[expr [pcrd p.wll.outvrd 1]-$inlvrd*0.6],0"
	ic_point projcurv GEOM p.inlvrd.helper.1 {p.inlvrd.helper.0 c.wll.vrd}
 
 
 
 
 
log "Geometry creation completed"
# -----------------------------------------------
log "Set Blocking"
# -----------------------------------------------
	# Inits (not all needed, untested)
	ic_geo_new_family FLUID
	#ic_boco_set_part_color FLUID
	ic_hex_initialize_mesh 2d new_numbering new_blocking FLUID
	# Creating FLUID blocking
# 	ic_hex_rename_blocking root fluid1
# 	ic_hex_switch_blocking fluid1
	ic_hex_unblank_blocks 
	ic_hex_multi_grid_level 0
	ic_hex_projection_limit 0
	ic_hex_default_bunching_law default 1.2
	ic_hex_floating_grid off
	ic_hex_transfinite_degree 1
	ic_hex_unstruct_face_type several_tris
	#ic_hex_set_unstruct_face_method uniform_quad
	ic_hex_set_n_tetra_smoothing_steps 20
	ic_hex_set_mesh_params GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID -version 110
	ic_hex_error_messages off_minor
# 	ic_hex_switch_blocking fluid1

	# Associate first edges and nodes
	ic_hex_move_node 11 p.sym.center
	ic_hex_move_node 13 p.sym.diffout
	ic_hex_move_node 19 p.int.center.help
	ic_hex_move_node 21 p.wll.diffout.10

	# outlet
	ic_hex_set_edge_projection 13 21 0 1 c.outlet
	ic_hex_split_grid 11 13 p.sym.diffsp m GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID 
	ic_hex_move_node 33 p.sym.diffsp
	ic_hex_move_node 34 p.wll.diffsp
	ic_hex_set_edge_projection 13 33 0 1 c.sym.diff
	ic_hex_set_edge_projection 21 34 0 1 c.wll.diff

 
 
	# Diffusor space split
	ic_hex_split_grid 11 33 p.disk.1.1 m GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID 
	ic_hex_move_node 37 p.disk.1.1
# 	ic_hex_move_node 38 p.wll.out
	ic_hex_set_edge_projection 33 37 0 1 c.sym.diffsp
# 	ic_hex_set_edge_projection 34 38 0 1 c.wll.diffsp

 
 
	# Interface split
	ic_hex_split_grid 37 38 p.int.out m GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID 

 
	# Disks
	for {set x 0} {$x < [expr $dsknm-1]} {incr x} {
		ic_hex_mark_blocks unmark
		ic_hex_mark_blocks superblock 4
		ic_hex_split_grid 37 "[expr 42+$x*12]" "p.disk.[expr $dsknm-$x].4" m GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID marked
		ic_hex_mark_blocks unmark
		ic_hex_mark_blocks superblock 4
		ic_hex_split_grid 37 "[expr 48+$x*12]" "p.disk.[expr $dsknm-$x].1" m GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID marked
		ic_hex_mark_blocks unmark
	}
		ic_hex_mark_blocks unmark
		ic_hex_mark_blocks superblock 4
		ic_hex_split_grid 37 "[expr 30+$dsknm*12]" "p.disk.1.4" m GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID marked
		ic_hex_mark_blocks unmark

 
	# Inner radius split 
	ic_hex_split_grid 41 42 p.disk.1.2 m SYM DISKS WALL FLUID 

 
	# Dynamic Vertex number of discribed points
	set p_disk_x_1_ver {if "$x!=1" "expr 36+($dsknm-$x+1)*12+6" else "expr 37"}
	set p_disk_x_2_ver {expr 41+$dsknm*12+($x-1)*2}
	set p_disk_x_3_ver {expr 42+$dsknm*12+($x-1)*2}
	set p_disk_x_4_ver {expr 36+($dsknm-$x+1)*12}

	# dynamic disks association
	for {set x 1} {$x < [expr $dsknm+1]} {incr x} {
		ic_hex_move_node [eval $p_disk_x_1_ver] "p.disk.$x.1"
		ic_hex_set_edge_projection [eval $p_disk_x_1_ver] [eval $p_disk_x_2_ver] 0 1 "c.disk.$x.1"		;# (p.disk.x.1 p.disk.x.2)
		ic_hex_move_node [eval $p_disk_x_2_ver] "p.disk.$x.2"
		ic_hex_set_edge_projection [eval $p_disk_x_2_ver] [eval $p_disk_x_3_ver] 0 1 "c.disk.$x.2"		;# (p.disk.x.2 p.disk.x.3)
		ic_hex_move_node [eval $p_disk_x_3_ver] "p.disk.$x.3"
		ic_hex_set_edge_projection [eval $p_disk_x_3_ver] [eval $p_disk_x_4_ver] 0 1 "c.disk.$x.3"		;# (p.disk.x.3 p.disk.x.4)
 		ic_hex_move_node [eval $p_disk_x_4_ver] "p.disk.$x.4"
  	ic_hex_set_edge_projection [eval $p_disk_x_4_ver] [eval $p_disk_x_1_ver] 0 1 "c.disk.$x.4"
	}

	# inletregion association
	ic_hex_set_edge_projection 11 [set x 1;eval $p_disk_x_2_ver] 0 1 c.sym.inreg

	# Axis projections 
	set p_disk_x_2_to_axis_ver {expr 35+($dsknm-$x+1)*12+6}
	set p_disk_x_3_to_axis_ver {expr 35+($dsknm-$x+1)*12}
	# Axis Associations
	ic_hex_set_edge_projection 11 [eval $p_disk_x_3_to_axis_ver] 0 1 c.axis
	ic_hex_set_edge_projection [eval $p_disk_x_3_to_axis_ver] [set x 2;eval $p_disk_x_2_to_axis_ver] 0 1 c.axis
	for {set y 2} {$y <= $dsknm} {incr y} {
		ic_hex_set_edge_projection [set x [expr $y-1]; eval $p_disk_x_3_to_axis_ver] [set x $y;eval $p_disk_x_2_to_axis_ver] 0 1 c.axis
		ic_hex_set_edge_projection [eval $p_disk_x_2_to_axis_ver] [eval $p_disk_x_3_to_axis_ver] 0 1 c.axis
	}
	ic_hex_set_edge_projection 41 47 0 1 c.axis
	ic_hex_set_edge_projection 41 19 0 1 c.axis
	
 
 
 
	# O-grids um disks
	for {set x 0} {$x < $dsknm} {incr x} {
		ic_hex_mark_blocks unmark
		ic_hex_mark_blocks superblock "[expr 13+$dsknm*2+$x*2]"
		ic_hex_ogrid 1 not_marked m GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID -version 50
		ic_hex_mark_blocks unmark
	}

	# BL Thickness setting
	set p_disk_x_1_bl_ver {expr 45+$dsknm*14+($x-1)*4}
	set p_disk_x_2_bl_ver {expr 44+$dsknm*14+($x-1)*4}
	set p_disk_x_3_bl_ver {expr 46+$dsknm*14+($x-1)*4}
	set p_disk_x_4_bl_ver {expr 47+$dsknm*14+($x-1)*4}
	set blsqth [expr $blth/sqrt(2.)]
	for {set x 1} {$x <= $dsknm} {incr x} {
		ic_hex_set_node_location x [expr [pcrd p.disk.$x.1 0]-$blth] y [expr [pcrd p.disk.$x.1 1]+$blsqth] -csys global node_numbers "\{[eval $p_disk_x_1_bl_ver]\}"
		ic_hex_set_node_location x [expr [pcrd p.disk.$x.2 0]-$blth] y [expr [pcrd p.disk.$x.2 1]-$blsqth] -csys global node_numbers "\{[eval $p_disk_x_2_bl_ver]\}"
		ic_hex_set_node_location x [expr [pcrd p.disk.$x.3 0]+$blth] y [expr [pcrd p.disk.$x.3 1]-$blsqth] -csys global node_numbers "\{[eval $p_disk_x_3_bl_ver]\}"
		ic_hex_set_node_location x [expr [pcrd p.disk.$x.4 0]+$blth] y [expr [pcrd p.disk.$x.4 1]+$blsqth] -csys global node_numbers "\{[eval $p_disk_x_4_bl_ver]\}"
		ic_hex_set_node_location x [expr [pcrd p.disk.$x.2 0]-$blth] -csys global node_numbers "\{[eval $p_disk_x_2_to_axis_ver]\}"
		ic_hex_set_node_location x [expr [pcrd p.disk.$x.3 0]+$blth] -csys global node_numbers [eval $p_disk_x_3_to_axis_ver]
	}
set x 1
		ic_hex_set_node_location x 0 y [expr [pcrd p.disk.$x.1 1]] -csys global node_numbers "\{[eval $p_disk_x_1_ver]\}"
		ic_hex_set_node_location x 0 y [expr [pcrd p.disk.$x.2 1]] -csys global node_numbers "\{[eval $p_disk_x_2_ver]\}"
		ic_hex_set_node_location x 0 y [expr [pcrd p.disk.$x.1 1]+$blth] -csys global node_numbers "\{[eval $p_disk_x_1_bl_ver]\}"
		ic_hex_set_node_location x 0 y [expr [pcrd p.disk.$x.2 1]-$blth] -csys global node_numbers "\{[eval $p_disk_x_2_bl_ver]\}"

 
 
	# setting mesh between disk BL and interior boundary
	if {$dskwllsp > $dsksp} {
	  set tempx [expr [pcrd p.wll.out 0]-$blth]
	  set tempy [expr [pcrd p.disk.$dsknm.4 1]+$blth]
		set tempx_3 [expr [pcrd p.disk.$dsknm.3 0]+$blth]
		set tempy_3 [expr [pcrd p.disk.$dsknm.3 1]-$blth]
		set tempx_4 [expr [pcrd p.disk.$dsknm.4 0]+$blth]
		set tempy_4 [expr [pcrd p.disk.$dsknm.4 1]+$blth]
	} else {
	  set tempx [expr [pcrd p.wll.out 0]-$dskwllsp/4.]
	  set tempy [expr [pcrd p.disk.$dsknm.4 1]+$dskwllsp/4.]
		set tempx_3 [expr [pcrd p.disk.$dsknm.3 0]+$dskwllsp/4.]
		set tempy_3 [expr [pcrd p.disk.$dsknm.3 1]-$dskwllsp/4.]
		set tempx_4 [expr [pcrd p.disk.$dsknm.4 0]+$dskwllsp/4.]
		set tempy_4 [expr [pcrd p.disk.$dsknm.4 1]+$dskwllsp/4.]		
 
	}	
 
 
 
	ic_hex_set_node_location x $tempx y $tempy -csys global node_numbers "\{42\}"	
	ic_hex_set_node_location x $tempx -csys global node_numbers "\{43\}"
	ic_hex_set_node_location y $tempy -csys global node_numbers "\{38\}"		
	ic_hex_set_node_location x $tempx -csys global node_numbers "\{41\}"	
	set x $dsknm
	ic_hex_set_node_location x $tempx_3 y $tempy_3 -csys global node_numbers "\{[eval $p_disk_x_3_bl_ver]\}"
	ic_hex_set_node_location x $tempx_4 y $tempy_4 -csys global node_numbers "\{[eval $p_disk_x_4_bl_ver]\}"
	ic_hex_set_node_location x $tempx_3 -csys global node_numbers "\{[eval $p_disk_x_3_to_axis_ver]\}"
	ic_hex_set_node_location x [expr [pcrd p.wll.outvrd 0]-$dskwllsp/4.] y $tempy_3 -csys global node_numbers "\{[expr [eval $p_disk_x_3_ver]+1]\}"	
	ic_hex_set_node_location y $tempy_3 -csys global node_numbers "\{[expr [eval $p_disk_x_3_ver]+2]\}"	
	
 
	# Additional Splitting of blocks to reduce skewed mesh between disks
		ic_hex_split_grid [set x 1; eval $p_disk_x_2_ver] [set x 1;eval $p_disk_x_1_ver] 0.9 m SYM DISKS WALL FLUID VORFN
 
		set p_disk_x_1_cw_ver {expr 41 + $dsknm*18+$x*4}
		set p_disk_x_4_cw_ver {expr 43 + $dsknm*18+$x*4}	
		set p_disk_x_1_cbl_ver {expr 42 + $dsknm*18+$x*4}
		set p_disk_x_4_cbl_ver {expr 44 + $dsknm*18+$x*4}	
 
		ic_hex_split_grid [set x 1; eval $p_disk_x_2_ver] [set x 1;eval $p_disk_x_1_cw_ver] 0.1 m SYM DISKS WALL FLUID VORFN	
 
		set p_disk_x_2_cw_ver {expr 45 + $dsknm*22+$x*4}
		set p_disk_x_3_cw_ver {expr 47 + $dsknm*22+$x*4}	
		set p_disk_x_2_cbl_ver {expr 46 + $dsknm*22+$x*4}
		set p_disk_x_3_cbl_ver {expr 48 + $dsknm*22+$x*4}	

 
		# Vertices setting to
		for {set x 1} {$x <= $dsknm} {incr x} {
			set radlgh [expr [pcrd p.disk.$x.1 1]-[pcrd p.disk.$x.2 1]]
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.1 1]-0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_1_cw_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.2 1]+0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_2_cw_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.3 1]+0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_3_cw_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.4 1]-0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_4_cw_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.1 1]-0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_1_cbl_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.2 1]+0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_2_cbl_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.3 1]+0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_3_cbl_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$x.4 1]-0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_4_cbl_ver]\}"	
			}
			
 
		# Vertices at Wall
			set x [expr $dsknm+1]; set radlgh [expr [pcrd p.disk.$dsknm.1 1]-[pcrd p.disk.$dsknm.2 1]]
			ic_hex_set_node_location y [expr [pcrd p.disk.$dsknm.4 1]-0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_1_cw_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$dsknm.4 1]-0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_1_cbl_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$dsknm.3 1]+0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_2_cw_ver]\}"
			ic_hex_set_node_location y [expr [pcrd p.disk.$dsknm.3 1]+0.05*$radlgh] -csys global node_numbers "\{[eval $p_disk_x_2_cbl_ver]\}"
	
 
 
 
	# -----------------------------------------------
	# Making the Inlet Tube
	# -----------------------------------------------
	
 
	#  Split to capture fillet at 80% of its radius
	set x $dsknm
	ic_hex_split_grid 19 [expr [eval $p_disk_x_3_ver]+2] p.inlvrd.helper.1 m SYM DISKS WALL INLET OUTLET AXIS FLUID
 
 
	set p_inlvrd_helper_1 {expr 54+$dsknm*28}
	set p_inlvrd_helper_1_bl {expr 53+$dsknm*28}
	# preparing Inlet tube blocking
	ic_hex_mark_blocks unmark
	ic_hex_mark_blocks superblock 11
	ic_hex_mark_blocks numbers 19 [eval $p_inlvrd_helper_1] edge_neighbors
	ic_hex_mark_blocks numbers 19 41 edge_neighbors
	ic_hex_ogrid 1 m SYM DISKS WALL INLET OUTLET AXIS FLUID -version 50
	ic_hex_mark_blocks unmark
 
	ic_hex_delete_blocks numbers 11
 
	set p_inl_center {expr 57+28*$dsknm}
	set p_inl_bl {expr 58+28*$dsknm}
	set p_wll_inl {expr 59+28*$dsknm}
	set p_wll_rin {expr 42+$dsknm*14}

 
	# Inlet Wall Projections
	ic_hex_create_composite {c.wll.diffsp c.wll.disk c.wll.vrd c.wll.inl}
	ic_hex_set_edge_projection [eval $p_inlvrd_helper_1] [eval $p_wll_rin] 0 1 c.wll.diffsp
	ic_hex_set_edge_projection [eval $p_inlvrd_helper_1] [eval $p_wll_inl] 0 1 c.wll.diffsp
	ic_hex_set_edge_projection [eval $p_wll_rin] [expr 50+$dsknm*26] 0 1 c.wll.diffsp
	ic_hex_set_edge_projection [expr 46+$dsknm*22] [expr 50+$dsknm*26] 0 1 c.wll.diffsp
	ic_hex_set_edge_projection [expr 46+$dsknm*22]  38 0 1 c.wll.diffsp
	ic_hex_set_edge_projection 38 34 0 1 c.wll.diffsp
 
	ic_hex_move_node [eval $p_inlvrd_helper_1] p.inlvrd.helper.1
	
	# Inlet Projections
	ic_hex_set_edge_projection [eval $p_inl_bl] [eval $p_wll_inl] 0 1 c.inlet
	ic_hex_set_edge_projection [eval $p_inl_center] [eval $p_inl_bl] 0 1 c.inlet
 
	ic_hex_move_node [eval $p_inl_center] p.inl.center
	ic_hex_move_node [eval $p_wll_inl] p.wll.inl
	
	# Bl Thickness in Inlet
	ic_hex_place_node [eval $p_inl_bl] curve:c.inlet 0.95
	# Link Inlet BL to Inlet Wall
	ic_hex_link_shape [expr 53+28*$dsknm] [eval $p_inl_bl] [eval $p_inlvrd_helper_1] [eval $p_wll_inl] 1.0

	# Fillet Split
	ic_hex_split_grid [eval $p_inlvrd_helper_1] [eval $p_wll_inl] p.wll.inlvrd m SYM DISKS WALL INLET OUTLET AXIS FLUID
 
	set p_wll_inlvrd {expr 84+32*$dsknm}
	set p_wll_inlvrd_bl {expr 83+32*$dsknm}
	set p_wll_inlvrd_center {expr 81+32*$dsknm}	
 
	ic_hex_set_edge_projection [eval $p_wll_inlvrd] [eval $p_wll_inlvrd_bl] 0 1 c.pinlet
	ic_hex_set_edge_projection [eval $p_wll_inlvrd_bl] [eval $p_wll_inlvrd_center] 0 1 c.pinlet
		
 
	# Fillet Split projections
	ic_hex_move_node [eval $p_wll_inlvrd] p.wll.inlvrd
	set tempx [expr [pcrd p.wll.inlvrd 0]]
	ic_hex_set_node_location x $tempx y [expr [pcrd p.wll.inlvrd 1]-$inlrd*0.05] -csys global node_numbers "\{[eval $p_wll_inlvrd_bl]\}"
	ic_hex_set_node_location x $tempx -csys global node_numbers "\{[eval $p_wll_inlvrd_center]\}"
	ic_hex_set_node_location x [expr [pcrd p.inlvrd.helper.1 0]-$inlrd*0.05] -csys global node_numbers "\{[eval $p_inlvrd_helper_1_bl]\}"
		# Get temp3 from previous calculation (setting mesh between disk BL and interior boundary)
	ic_hex_set_node_location x $tempx_3 -csys global node_numbers "\{[expr [eval $p_inlvrd_helper_1_bl]-1]\}"
	ic_hex_set_node_location x [expr [pcrd p.inlvrd.helper.1 0]] -csys global node_numbers 41	
 
	ic_hex_set_node_location y [pcrd p.inlvrd.helper.1 1] -csys global node_numbers "\{[expr 53+$dsknm*26]\}"	
	
 
	# Remaining Disk blocks deletion
	for {set x 0} {$x<$dsknm} {incr x} {
	  ic_hex_delete_blocks numbers "[expr 13+$dsknm*2+$x*2]"
		if {$x == [expr $dsknm-1]} {
			ic_hex_delete_blocks numbers "[expr 17+$dsknm*8+$x*4]"
		} else {
			ic_hex_delete_blocks numbers "[expr 16+$dsknm*8+$x*4]"
		}
	  ic_hex_delete_blocks numbers "[expr 19+$dsknm*12+$x*4]"		
	}

	# Messflaechen
	ic_hex_set_edge_projection 33 [expr 37+12*$dsknm] 0 1 c.poutlet
	ic_hex_set_edge_projection 43 34 0 1 c.poutlet	
	set x 1
	ic_hex_set_mesh 33 [expr 37+12*$dsknm] linked [eval $p_disk_x_1_bl_ver] [eval $p_disk_x_4_bl_ver] locked
 
 
 
	ic_hex_link_shape 43 44 34 21 1.0	
 
log "Blocking Completed"

# -----------------------------------------------		
log "Edge settings"
# -----------------------------------------------
	for {set x 1} {$x <= $dsknm} {incr x} {
		# Disk BL
		ic_hex_set_mesh [eval $p_disk_x_2_bl_ver] [eval $p_disk_x_2_ver]  n $nbl h1rel 0.00 h2rel 0.04 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked 
		# Disk thickness direction inner radius ( problem with first disk: copy to parallel)
			if {$x == 1} {ic_hex_set_mesh [eval $p_disk_x_3_ver] [eval $p_disk_x_2_ver] n "[ceil [expr $ndskth*([pcrd p.disk.$x.3 0]-[pcrd p.disk.$x.2 0])+5]]" h1rel 0.0 h2rel 0.0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
			} else {	  	ic_hex_set_mesh [eval $p_disk_x_2_ver] [eval $p_disk_x_3_ver] n "[ceil [expr $ndskth*([pcrd p.disk.$x.3 0]-[pcrd p.disk.$x.2 0])+5]]" h1rel 0.0 h2rel 0.0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked}
		# Disk thickness direction outer radius
		ic_hex_set_mesh [eval $p_disk_x_1_ver] [eval $p_disk_x_4_ver] n "[ceil [expr $ndskth*([pcrd p.disk.$x.4 0]-[pcrd p.disk.$x.1 0])+5]]" h1rel 0.0 h2rel 0.0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked	
		# Disk Surface  in three parts
		ic_hex_set_mesh [eval $p_disk_x_4_cw_ver ] [eval $p_disk_x_4_ver] n "[ceil [expr $ndsksp*([pcrd p.disk.$x.4 1]-[pcrd p.disk.$x.3 1])*0.05+5]]" h1rel 0.0 h2rel linked [eval $p_disk_x_1_ver ] [eval $p_disk_x_4_ver ] r1 1.2 r2 1.05 lmax 0 default copy_to_parallel unlocked
		ic_hex_set_mesh [eval $p_disk_x_3_cw_ver] [eval $p_disk_x_4_cw_ver] n "[ceil [expr $ndsksp*([pcrd p.disk.$x.4 1]-[pcrd p.disk.$x.3 1])*0.3+5]]" h1rel linked [eval $p_disk_x_3_cw_ver] [eval $p_disk_x_3_ver] h2rel linked [eval $p_disk_x_4_cw_ver] [eval $p_disk_x_4_ver] r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
		ic_hex_set_mesh [eval $p_disk_x_3_ver] [eval $p_disk_x_3_cw_ver] n "[ceil [expr $ndsksp*([pcrd p.disk.$x.4 1]-[pcrd p.disk.$x.3 1])*0.05+5]]" h1rel linked [eval $p_disk_x_2_ver ] [eval $p_disk_x_3_ver ] h2rel 0.0 r1 1.05 r2 1.2 lmax 0 default copy_to_parallel unlocked
	}
	
 
 
	# Disk space
	for {set y 1} {$y < $dsknm} {incr y} {
		set x $y; set	dskminusone [eval $p_disk_x_4_bl_ver ]
		set x [expr $y+1]; set dskone [eval $p_disk_x_1_bl_ver ]
		ic_hex_set_mesh $dskminusone $dskone n "[ceil [expr $ndsksp*($dsksp-2*$blth)+4]]" h1rel 0 h2rel 0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
 
	}
	
	# Inlet region upper part
	set x 1
	set nn [ceil [expr $ninlreg*([pcrd p.disk.$dsknm.3 1] - [pcrd p.inlvrd.helper.1 1])*2+5]]
	ic_hex_set_mesh [eval $p_inlvrd_helper_1] [eval $p_wll_rin] n $nn h2rel linked [eval $p_disk_x_2_bl_ver] [eval $p_disk_x_2_ver] h1rel 0 r1 1.2 r2 1.05 lmax 0 default copy_to_parallel unlocked
	
	# Inlet region lower part
	set nn [ceil [expr $ninlreg*([pcrd p.inlvrd.helper.1 1])+5]]
	set x 1
	ic_hex_set_mesh 11 [expr 53+$dsknm*26]  n $nn h2rel linked [expr 53+$dsknm*26] [eval $p_disk_x_2_bl_ver] h1rel 0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
		
	# Diffusor space 
	set x 1
	ic_hex_set_mesh [eval $p_disk_x_1_bl_ver] 33 n [ceil [expr $diffsp*$ndiffsp+5]] h1rel linked [eval $p_disk_x_1_bl_ver] [eval $p_disk_x_1_ver] h2rel 0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
	# Diffusor 
	ic_hex_set_mesh 34 21 n [ceil [expr $difflgh*$ndifflgh]] h1rel linked 34 38 h2rel 0.001 r1 1.05 r2 1.05 lmax 0 default copy_to_parallel unlocked
	
 
	# Wall Bl
	set x $dsknm
	ic_hex_set_mesh [expr [eval $p_disk_x_3_ver]+1] [expr [eval $p_disk_x_3_ver]+2] n $nbl h2rel 0.05 h1rel 0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
	# Disk Wall space
	ic_hex_set_mesh [eval $p_disk_x_3_bl_ver] [expr [eval $p_disk_x_3_ver]+1] n [ceil [expr $ndsksp*$dskwllsp/2.+18]] h2rel 0.00 h1rel 0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
	
 
	# Fillet
	ic_hex_set_mesh [eval $p_inlvrd_helper_1] [eval $p_wll_inlvrd] n [ceil [expr $nfillet*$pi*$inlvrd+5]] h1rel linked [eval $p_inlvrd_helper_1] [expr [eval $p_disk_x_3_ver]+2] h2rel 0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
	ic_hex_link_shape [eval $p_inlvrd_helper_1_bl] [eval $p_wll_inlvrd_bl] [eval $p_inlvrd_helper_1] [eval $p_wll_inlvrd] 1.0
	
 
 
	# Inlet length
	ic_hex_set_mesh [eval $p_wll_inlvrd] [eval $p_wll_inl] n [ceil [expr $ninllgh*$inllgh]] h1rel linked [eval $p_wll_inlvrd] [eval $p_inlvrd_helper_1] h2rel 0 r1 1.2 r2 1.2 lmax 0 default copy_to_parallel unlocked
	
 
# -----------------------------------------------
log "Meshing"
# -----------------------------------------------
 if {$fast == 0} {
		# Boundary Conditions
		ic_boco_solver Fluent_V6
		ic_solver_mesh_info Fluent_V6		
		ic_boco_set WALL { { 1  {WALL}  0  } }
		ic_boco_set INLET { { 1  {MASFI}  0  } }
		ic_boco_set FLUID { { 1  {FLUID}  0  } }
		ic_boco_set SYM { { 1  {SYM}  0  } }
		ic_boco_set OUTLET { { 1  {PRESO}  0  } }
		ic_boco_set AXIS { { 1  {AXIS}  0  } }
		ic_boco_set DISKS { { 1  {WALL}  0  } }
#		ic_boco_set PINLET {{1 INTER 0}}
#		ic_boco_set POUTLET {{1 INTER 0}}
 		ic_boco_set PINLET {{1 FAN 0}}
 		ic_boco_set POUTLET {{1 FAN 0}} 
 
    # Convert to Unstruct
  #  cmd_rm "hex.uns"
 #   ic_hex_write_file hex.uns GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID proj 2 dim_to_mesh 2 -family_boco $fluentfile.fbc
    set determinante [split [ic_hex_rh 20 FLUID proj 2 minval 0 -type determinant_27 maxval 1 new_format] " "]
    set skewness [split [ic_hex_rh 20 FLUID proj 2 minval 0 -type eriksson maxval 1 new_format] " "]
    log "Worst Determinant: [lindex $determinante 0]"
    log "Worst Skewness: [lindex $skewness 0]"

    # Try Smoothing
    ic_undo_group_begin 
    ic_hex_smooth 5 GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID elliptic iter_srf 15 iter_vol 5 exp_srf 0.0 exp_vol 0.0 smooth_type 202 nfix_layers -1 rebunch_edges 0 treat_unstruct 0 stabilize_srf 1.0 stabilize_vol 2.0 ortho_distance_srf 0 ortho_distance_vol 0 surface_fitting 1 keep_per_geom 1
    ic_undo_group_end

    # Smoothing succeded?
    set determinante [split [ic_hex_rh 20 FLUID proj 2 minval 0 -type determinant_27 maxval 1 new_format] " "]
    set skewness [split [ic_hex_rh 20 FLUID proj 2 minval 0 -type eriksson maxval 1 new_format] " "]
    if {[lindex $determinante 0]<0.4} {
      ic_undo
      log "No smoothing applied. Determinante too low, [lindex $determinante 0]"
     } else {
      log "After Smoothing Worst Determinant: [lindex $determinante 0]"
    }
    if {[lindex $skewness 0]<0.3} {
      ic_undo
      log "No smoothing applied. Skewness not good enough, [lindex $skewness 0]"
     } else {
      log "After Smooting Worst Skewness: [lindex $skewness 0]"  
    }
	  # Convert Final to Unstruct	
	  cmd_rm "hex.uns"	
    ic_hex_write_file hex.uns GEOM SYM DISKS WALL INLET OUTLET AXIS FLUID proj 2 dim_to_mesh 2 -family_boco $fluentfile.fbc
 
 
log "#################################################" 0
log "Writing mesh" 0
log "#################################################" 0
log "Export Command: $icemenv/icemcfd/output-interfaces/fluent6 -dom hex.uns -b $fluentfile.fbc -dim2d -scale 0.001,0.001,0.001 $fluentfile.msh"
    exec "$icemenv/icemcfd/output-interfaces/fluent6" -dom hex.uns -b $fluentfile.fbc -dim2d -scale 0.001,0.001,0.001 $fluentfile.msh
    log "Mesh exported to $fluentfile.msh"
  } else {
    log "Export skipped."
  }
 
cmd_rm "hex.uns"		
log "Script finshed!"
log "#################################################" 0

in_parameter.tcl

in_parameter.tcl
set dsksp $GAP
set dskrdout $RDOUT
set dskrdin $RDIN
set dskwllsp $WALLSPACE
set rho 1059.0
set my 0.0036
set ekman 2.0	
set mflow 0.00842725
set targetpressure 13332.23

in_extractvalues.sh

in_extractvalues.sh
#!/bin/bash
# Bash-Script
# Part of Work: Optimizing the hydraulic efficiency of a Tesla Pump
# Rotor Independent Optimization
# Read Last line value and write to corresponding file
 
if [ "$1" = "presDrop" ]
then
	echo "$(tail -n 1  < presDrop.csv | cut -d',' -f2)" > presDrop.val
else
	if [ "$1" = "torque" ]
	then
		echo "$(tail -n 1  < torque.csv | cut -d',' -f2)" > torque.val
	else
		if [ "$1" = "speed" ]
		then
			echo "$(tail -n 1  < rotationSpeed.csv | cut -d',' -f2)" > speed2.val
		fi
	fi
fi

in_tesrot.java

in_tesrot.java
// STAR-CCM+ macro: tesopt2.java
// Written by STAR-CCM+ 9.04.009
package macro;
 
import java.util.*;
 
import star.turbulence.*;
import star.kwturb.*;
import star.material.*;
import star.common.*;
import star.base.neo.*;
import star.vis.*;
import star.base.report.*;
import star.flow.*;
import star.segregatedflow.*;
import star.metrics.*;
import star.energy.*;
 
public class in_tesrot extends StarMacro {
 
  public void execute() {
 
    prep();
    solve();
    save();
 
  }
 
  private void prep() {
 
    Simulation sim = getActiveSimulation();  
 
    // Mesh import
    ImportManager importManager_0 = sim.getImportManager();
    importManager_0.importMeshFiles(new StringVector(new String[] {resolvePath("slice.msh")}), NeoProperty.fromString("{\'FileOptions\': [{\'Sequence\': 50}]}"));
 
    FvRepresentation fvRepresentation_0 = ((FvRepresentation) sim.getRepresentationManager().getObject("Volume Mesh"));
 
    Region fluidregion = sim.getRegionManager().getRegion("FLUID");
 
    fvRepresentation_0.generateMeshReport(new NeoObjectVector(new Object[] {fluidregion})); 
 
 
    // Setting physics
    PhysicsContinuum physicsContinuum_fluid = ((PhysicsContinuum) sim.getContinuumManager().getContinuum("Physics 1"));
 
    physicsContinuum_fluid.enable(SteadyModel.class);
    TwoDimensionalModel twoDimensionalModel_0 = physicsContinuum_fluid.getModelManager().getModel(TwoDimensionalModel.class);
    physicsContinuum_fluid.disableModel(twoDimensionalModel_0);        
    physicsContinuum_fluid.enable(AxisymmetricModel.class);
    physicsContinuum_fluid.enable(SingleComponentLiquidModel.class);
    physicsContinuum_fluid.enable(SegregatedFlowModel.class);
    physicsContinuum_fluid.enable(ConstantDensityModel.class);
    physicsContinuum_fluid.enable(TurbulentModel.class);
    physicsContinuum_fluid.enable(RansTurbulenceModel.class);
    physicsContinuum_fluid.enable(KOmegaTurbulence.class);
    physicsContinuum_fluid.enable(SstKwTurbModel.class);
    physicsContinuum_fluid.enable(KwAllYplusWallTreatment.class);
 
    // Low Re Wall function
    KwAllYplusWallTreatment kwAllYplusWallTreatment_0 = physicsContinuum_fluid.getModelManager().getModel(KwAllYplusWallTreatment.class);
    physicsContinuum_fluid.disableModel(kwAllYplusWallTreatment_0);
    physicsContinuum_fluid.enable(KwLowYplusWallTreatment.class);
    physicsContinuum_fluid.enable(AxisymmetricSwirlModel.class);   
 
    // Fluid properties
    SingleComponentLiquidModel fluidModel = physicsContinuum_fluid.getModelManager().getModel(SingleComponentLiquidModel.class);
    Liquid liquid_water = ((Liquid) fluidModel.getMaterial());
 
    ConstantMaterialPropertyMethod fluidMaterialProperty_density = ((ConstantMaterialPropertyMethod) liquid_water.getMaterialProperties().getMaterialProperty(ConstantDensityProperty.class).getMethod());
    fluidMaterialProperty_density.getQuantity().setValue(RHOKGKM);
 
    ConstantMaterialPropertyMethod fluidMaterialProperty_viscosity = ((ConstantMaterialPropertyMethod) liquid_water.getMaterialProperties().getMaterialProperty(DynamicViscosityProperty.class).getMethod());
    fluidMaterialProperty_viscosity.getQuantity().setValue(KINVISC); 
 
 
 
    // Boundary Conditions
    BoundaryInterface bdryInt_inlet = ((BoundaryInterface) sim.getInterfaceManager().getInterface("PINLET"));
    bdryInt_inlet.setInterfaceType(InternalInterface.class);
 
    BoundaryInterface bdryInt_outlet = ((BoundaryInterface) sim.getInterfaceManager().getInterface("POUTLET"));
    bdryInt_outlet.setInterfaceType(InternalInterface.class);
 
    Boundary boundary_disks = fluidregion.getBoundaryManager().getBoundary("DISKS");
    boundary_disks.getConditions().get(WallSlidingOption.class).setSelected(WallSlidingOption.ROTATION_RATE);
    WallRotationProfile wallRotationProfile = boundary_disks.getValues().get(WallRotationProfile.class);
    wallRotationProfile.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(RPMVALUE);
 
    Boundary boundary_inlet = fluidregion.getBoundaryManager().getBoundary("INLET");
    boundary_inlet.getConditions().get(FlowDirectionOption.class).setSelected(FlowDirectionOption.COMPONENTS);
    FlowDirectionProfile flowDirectionProfile_0 = boundary_inlet.getValues().get(FlowDirectionProfile.class);
    flowDirectionProfile_0.getMethod(ConstantVectorProfileMethod.class).getQuantity().setComponents(-1.0, 0.0, 0.0);     
    MassFlowRateProfile massFlowRateProfile_0 = boundary_inlet.getValues().get(MassFlowRateProfile.class);
    massFlowRateProfile_0.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(MASSFLOW);
 
    Boundary boundary_outlet = fluidregion.getBoundaryManager().getBoundary("OUTLET");
    boundary_outlet.getConditions().get(BackflowDirectionOption.class).setSelected(BackflowDirectionOption.EXTRAPOLATED);
 
    // Set Rotation Axis
    ReferenceFrame referenceFrame_0 = fluidregion.getValues().get(ReferenceFrame.class);
    referenceFrame_0.getAxisVector().setComponents(1.0, 0.0, 0.0); 
 
 
 
    // Make Reports, Monitors
 
      // PressureDrop
      PressureDropReport pressureDropReport = sim.getReportManager().createReport(PressureDropReport.class);      
      InterfaceBoundary interfaceBoundary_poutlet = ((InterfaceBoundary) fluidregion.getBoundaryManager().getBoundary("POUTLET [POUTLET]"));
      InterfaceBoundary interfaceBoundary_pinlet = ((InterfaceBoundary) fluidregion.getBoundaryManager().getBoundary("PINLET [PINLET]"));
      pressureDropReport.getParts().setObjects(interfaceBoundary_poutlet);
      pressureDropReport.getLowPressureParts().setObjects(interfaceBoundary_pinlet);
      pressureDropReport.setPresentationName("presDrop");      
 
      sim.getMonitorManager().createMonitorAndPlot(new NeoObjectVector(new Object[] {pressureDropReport}), true, "%1$s Plot");
      ReportMonitor reportMonitor_presDrop = ((ReportMonitor) sim.getMonitorManager().getMonitor("presDrop Monitor"));
      MonitorPlot monitorPlot_presDrop= sim.getPlotManager().createMonitorPlot(new NeoObjectVector(new Object[] {reportMonitor_presDrop}), "presDrop Monitor Plot");
 
      // Pressure 
      PrimitiveFieldFunction primitiveFieldFunction_1 = ((PrimitiveFieldFunction) sim.getFieldFunctionManager().getFunction("Pressure"));
	// Inlet
	AreaAverageReport areaAverageReport_presInlet = sim.getReportManager().createReport(AreaAverageReport.class);
	areaAverageReport_presInlet.setScalar(primitiveFieldFunction_1);
	areaAverageReport_presInlet.setPresentationName("presInlet");	
	areaAverageReport_presInlet.getParts().setObjects(boundary_inlet);	
 
	sim.getMonitorManager().createMonitorAndPlot(new NeoObjectVector(new Object[] {areaAverageReport_presInlet}), true, "%1$s Plot");
	ReportMonitor reportMonitor_presInlet = ((ReportMonitor) sim.getMonitorManager().getMonitor("presInlet Monitor"));
	MonitorPlot monitorPlot_presInlet = sim.getPlotManager().createMonitorPlot(new NeoObjectVector(new Object[] {reportMonitor_presInlet}), "presInlet Monitor Plot");     
 
	// Outlet
	AreaAverageReport areaAverageReport_presOutlet = sim.getReportManager().createReport(AreaAverageReport.class);
	areaAverageReport_presOutlet.setScalar(primitiveFieldFunction_1);
	areaAverageReport_presOutlet.setPresentationName("presOutlet");	
	areaAverageReport_presOutlet.getParts().setObjects(boundary_outlet);  
 
	sim.getMonitorManager().createMonitorAndPlot(new NeoObjectVector(new Object[] {areaAverageReport_presOutlet}), true, "%1$s Plot");
	ReportMonitor reportMonitor_presOutlet = ((ReportMonitor) sim.getMonitorManager().getMonitor("presOutlet Monitor"));
	MonitorPlot monitorPlot_2 = sim.getPlotManager().createMonitorPlot(new NeoObjectVector(new Object[] {reportMonitor_presOutlet}), "presOutlet Monitor Plot");	
 
      // Torque
      MomentReport momentReport_0 = sim.getReportManager().createReport(MomentReport.class);
      momentReport_0.setPresentationName("Torque");
      momentReport_0.getDirection().setComponents(1.0, 0.0, 0.0);
      momentReport_0.getParts().setObjects(boundary_disks);
 
      sim.getMonitorManager().createMonitorAndPlot(new NeoObjectVector(new Object[] {momentReport_0}), true, "%1$s Plot");
      ReportMonitor reportMonitor_Torque = ((ReportMonitor) sim.getMonitorManager().getMonitor("Torque Monitor"));
      MonitorPlot monitorPlot_Torque = sim.getPlotManager().createMonitorPlot(new NeoObjectVector(new Object[] {reportMonitor_Torque}), "Torque Monitor Plot");      
 
 
 
      // Speed
      UserFieldFunction currentSpeedFF = sim.getFieldFunctionManager().createFieldFunction();
      currentSpeedFF.getTypeOption().setSelected(FieldFunctionTypeOption.SCALAR);
      currentSpeedFF.setFunctionName("rotationSpeed");
      currentSpeedFF.setPresentationName("rotationSpeed");
      currentSpeedFF.setDefinition(String.valueOf(RPMVALUE));  
 
 
      ExpressionReport expressionReport_rotationSpeed = sim.getReportManager().createReport(ExpressionReport.class);
      expressionReport_rotationSpeed.setDefinition("$rotationSpeed");
      expressionReport_rotationSpeed.setPresentationName("rotationSpeed");
 
      ReportMonitor reportMonitor_rotationSpeed = expressionReport_rotationSpeed.createMonitor();     
 
      MonitorPlot monitorPlot_rotationSpeed = sim.getPlotManager().createMonitorPlot(new NeoObjectVector(new Object[] {reportMonitor_rotationSpeed}), "rotationSpeed Monitor Plot");
 
 
 
 
 
 
      // Solver Settings
 
	// Reduce Under-Relaxation for pressure
	SegregatedFlowSolver segregatedFlowSolver_0 = ((SegregatedFlowSolver) sim.getSolverManager().getSolver(SegregatedFlowSolver.class));
	PressureSolver pressureSolver_0 = segregatedFlowSolver_0.getPressureSolver();
	pressureSolver_0.setUrf(0.1);      
 
 
    Solution solution_0 = sim.getSolution();
    solution_0.initializeSolution();
 
 
      // Stopping Criteria
 
      ResidualMonitor residualMonitor_xMomentum = ((ResidualMonitor) sim.getMonitorManager().getMonitor("X-momentum"));
      ResidualMonitor residualMonitor_yMomentum = ((ResidualMonitor) sim.getMonitorManager().getMonitor("Y-momentum"));
      ResidualMonitor residualMonitor_zMomentum = ((ResidualMonitor) sim.getMonitorManager().getMonitor("Z-momentum"));
      ResidualMonitor residualMonitor_continuity = ((ResidualMonitor) sim.getMonitorManager().getMonitor("Continuity"));  
      ResidualMonitor residualMonitor_tke = ((ResidualMonitor) sim.getMonitorManager().getMonitor("Tke"));    
      ResidualMonitor residualMonitor_sdr = ((ResidualMonitor) sim.getMonitorManager().getMonitor("Sdr"));        
      MonitorIterationStoppingCriterion stopCrit_xMomentum = residualMonitor_xMomentum.createIterationStoppingCriterion();
      MonitorIterationStoppingCriterion stopCrit_yMomentum = residualMonitor_yMomentum.createIterationStoppingCriterion();
      MonitorIterationStoppingCriterion stopCrit_zMomentum = residualMonitor_zMomentum.createIterationStoppingCriterion();
      MonitorIterationStoppingCriterion stopCrit_continuity = residualMonitor_continuity.createIterationStoppingCriterion();
      MonitorIterationStoppingCriterion stopCrit_tke = residualMonitor_tke.createIterationStoppingCriterion();
      MonitorIterationStoppingCriterion stopCrit_sdr = residualMonitor_sdr.createIterationStoppingCriterion();    
 
 
      MonitorIterationStoppingCriterionMinLimitType stopCritLimit_0 = ((MonitorIterationStoppingCriterionMinLimitType) stopCrit_xMomentum.getCriterionType());
      MonitorIterationStoppingCriterionMinLimitType stopCritLimit_1 = ((MonitorIterationStoppingCriterionMinLimitType) stopCrit_yMomentum.getCriterionType());
      MonitorIterationStoppingCriterionMinLimitType stopCritLimit_2 = ((MonitorIterationStoppingCriterionMinLimitType) stopCrit_zMomentum.getCriterionType());    
      MonitorIterationStoppingCriterionMinLimitType stopCritLimit_3 = ((MonitorIterationStoppingCriterionMinLimitType) stopCrit_continuity.getCriterionType());
      MonitorIterationStoppingCriterionMinLimitType stopCritLimit_4 = ((MonitorIterationStoppingCriterionMinLimitType) stopCrit_tke.getCriterionType());    
      MonitorIterationStoppingCriterionMinLimitType stopCritLimit_5 = ((MonitorIterationStoppingCriterionMinLimitType) stopCrit_sdr.getCriterionType());
 
     StepStoppingCriterion maxStepStoppingCriterion = 
      ((StepStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Maximum Steps"));     
 
      maxStepStoppingCriterion.setMaximumNumberSteps(30000);          
 
      stopCritLimit_0.getLimit().setValue(1.0E-3);
      stopCritLimit_1.getLimit().setValue(1.0E-3);
      stopCritLimit_2.getLimit().setValue(1.0E-3);
      stopCritLimit_3.getLimit().setValue(1.0E-3);
      stopCritLimit_4.getLimit().setValue(1.0E-3);
      stopCritLimit_5.getLimit().setValue(1.0E-3);  
 
  }
 
  private void solve() {
 
 
    // Some Constants
    double targetPressure = 13332.23;
 
    //some Variables
    double rotationRate = RPMVALUE;
    double currentPressure;
    double jump;
    int nextIter = 1000;
 
    boolean presOK = false;
 
 
    // Get Report pressure Drop      
    Simulation sim = getActiveSimulation();      
    Report presDropReport = sim.getReportManager().getReport("presDrop");      
 
 
 
 
    // Get Wall Rotation Object
    Region fluidregion = sim.getRegionManager().getRegion("FLUID");    
    Boundary boundary_disks = fluidregion.getBoundaryManager().getBoundary("DISKS");    
    WallRotationProfile wallRotationProfile = boundary_disks.getValues().get(WallRotationProfile.class);
 
 
    // get Solver Settings
    //get Solver settings
    StepStoppingCriterion maxStepStoppingCriterion = 
      ((StepStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Maximum Steps"));
    MonitorIterationStoppingCriterion continuityStoppingCriterion = 
      ((MonitorIterationStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Continuity Criterion"));
    MonitorIterationStoppingCriterion xMomentumStoppingCriterion = 
      ((MonitorIterationStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("X-momentum Criterion"));  
    MonitorIterationStoppingCriterion yMomentumStoppingCriterion = 
      ((MonitorIterationStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Y-momentum Criterion"));
    MonitorIterationStoppingCriterion zMomentumStoppingCriterion = 
      ((MonitorIterationStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Z-momentum Criterion"));
    MonitorIterationStoppingCriterion tkeStoppingCriterion = 
      ((MonitorIterationStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Tke Criterion"));      
    MonitorIterationStoppingCriterion sdrStoppingCriterion = 
      ((MonitorIterationStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Sdr Criterion"));      
 
 
    UserFieldFunction currentSpeedFF = ((UserFieldFunction) sim.getFieldFunctionManager().getFunction("rotationSpeed"));
 
 
 
 
 
 
 
 
    // Init
    sim.getSimulationIterator().step(nextIter);
 
 
    // Solve
    do {
 
 
      // Get Pressuredrop
      currentPressure = presDropReport.getReportMonitorValue();
 
 
      // Calc new Roation Speed
      if (currentPressure > 0) 
      {
      jump = Math.sqrt(targetPressure/currentPressure)*rotationRate - rotationRate;
      rotationRate = rotationRate + jump*0.1;	
      nextIter = 100;
      } else {
	// If Turbine effect
	rotationRate = rotationRate + 2.0;
	nextIter = 100;
 
      }
 
       // set current speed
      wallRotationProfile.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(rotationRate);
      currentSpeedFF.setDefinition(String.valueOf(rotationRate));            
 
      // Advance solution
     sim.getSimulationIterator().step(nextIter); 
 
      // Get Pressuredrop
      currentPressure = presDropReport.getReportMonitorValue();
 
      if ( currentPressure < (targetPressure+10.0) && currentPressure > (targetPressure-10.0) )
      {
	presOK = true;
      } else
      {
	presOK = false;
      }
 
 
 
     } while(!maxStepStoppingCriterion.getIsSatisfied() && !(presOK && continuityStoppingCriterion.getIsSatisfied() && xMomentumStoppingCriterion.getIsSatisfied() && yMomentumStoppingCriterion.getIsSatisfied() && zMomentumStoppingCriterion.getIsSatisfied() && tkeStoppingCriterion.getIsSatisfied() && sdrStoppingCriterion.getIsSatisfied() ));
 
    // Disable Stopping Criteria
    StepStoppingCriterion stepStoppingCriterion_0 = ((StepStoppingCriterion) sim.getSolverStoppingCriterionManager().getSolverStoppingCriterion("Maximum Steps"));
    stepStoppingCriterion_0.setIsUsed(false);    
 
     sim.getSimulationIterator().step(5000);     
 
 
  }
 
  private void save() {
 
    Simulation sim = getActiveSimulation();  
 
    // Save Simulation data
    sim.saveState(resolvePath("teslaAxisym.sim"));    
 
 
    // Save Screens  
    MonitorPlot monitorPlot_0 = ((MonitorPlot) sim.getPlotManager().getPlot("presDrop Monitor Plot"));
    MonitorPlot monitorPlot_1 = ((MonitorPlot) sim.getPlotManager().getPlot("presInlet Monitor Plot"));    
    MonitorPlot monitorPlot_2 = ((MonitorPlot) sim.getPlotManager().getPlot("presOutlet Monitor Plot"));    
    MonitorPlot monitorPlot_3 = ((MonitorPlot) sim.getPlotManager().getPlot("Torque Monitor Plot"));    
    MonitorPlot monitorPlot_4 = ((MonitorPlot) sim.getPlotManager().getPlot("rotationSpeed Monitor Plot"));    
    ResidualPlot residualPlot_0 = ((ResidualPlot) sim.getPlotManager().getPlot("Residuals"));
 
 
    monitorPlot_0.export(resolvePath("presDrop.csv"), ",");  
    monitorPlot_1.export(resolvePath("presInlet.csv"), ","); 
    monitorPlot_2.export(resolvePath("presOutlet.csv"), ","); 
    monitorPlot_3.export(resolvePath("torque.csv"), ","); 
    monitorPlot_4.export(resolvePath("rotationSpeed.csv"), ","); 
 
    monitorPlot_0.encode(resolvePath("presDrop.png"), "png", 1280, 1024);
    monitorPlot_1.encode(resolvePath("presInlet.png"), "png", 1280, 1024);
    monitorPlot_2.encode(resolvePath("presOutlet.png"), "png", 1280, 1024);
    monitorPlot_3.encode(resolvePath("torque.png"), "png", 1280, 1024);
    monitorPlot_4.encode(resolvePath("rotationSpeed.png"), "png", 1280, 1024);
    residualPlot_0.encode(resolvePath("residuals.png"), "png", 1280, 1024);
 
  }
 
 
 
}

pbs-opal2-script.pbs

pbs-opal2-script.pbs
#!/bin/bash
#!
#! Sample PBS file
#!
#! Name of job
 
#PBS -N TeslaOptAxiSymMain
 
#! Number of nodes (in this case 8 nodes with 4 CPUs each)
#! The total number of nodes passed to mpirun will be nodes*ppn
#! Second entry: Total amount of wall-clock time (true time).
#! 02:00:00 indicates 02 hours
 
#PBS -l nodes=10:ppn=8,walltime=480:00:00
 
#! export all my environment variables to the job (and the loaded modules)
###PBS -V
 
#! Specify the queue
 
#######PBS -q short_eth
#PBS -q long_eth
 
#! Mail to user when job terminates or aborts
#PBS -m ae
#PBS -M sebastian.engel@st.ovgu.de
 
#module add fluent
 
#! Full path to application + application name
application="/home/engel/OPAL2/OPAL2/OPAL2"
 
#! Run options for the application
#options="SendRecv"
options="-opal2.mode=NUMERIC_NORMAL -opal2.file={/home/engel/Tesla/AxiSym/Optimization/Input/Tesla_Optimization.o2script}"
 
#! Work directory
workdir="/home/engel/Tesla/AxiSym/Optimization/Input"
 
 
 
###############################################################
### You should not have to change anything below this line ####
###############################################################
 
#! change the working directory (default is home directory)
 
cd $workdir
 
echo Running on host `hostname`
echo Time is `date`
echo Directory is `pwd`
echo PBS job ID is $PBS_JOBID
echo This jobs runs on the following machines:
echo `cat $PBS_NODEFILE | uniq`
 
#! Create a machine file for MPI
cat $PBS_NODEFILE | uniq > machine.file.$PBS_JOBID
 
numnodes=`wc $PBS_NODEFILE | awk '{ print $1 }'`
 
#! Run the parallel MPI executable (nodes*ppn)
 
#echo "Running mpirun -np $numnodes -machinefile machine.file.$PBS_JOBID $application $options"
mpirun -np 10 -machinefile machine.file.$PBS_JOBID $application $options > output_PBS.txt
echo "Opal started!"
#$application $options
guide/opal-scriptexample-tesla-icem-starccm.txt · Last modified: 2020/02/16 18:52 by seengel
Back to top
CC Attribution-Share Alike 3.0 Unported
chimeric.de = chi`s home Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0