research/incomp-flame-post
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
incomp-flame-post/code/code_gen/post.py
ignis 9ee9721c4c
All checks were successful
CI Test Suite / run-tests (push) Successful in 1m11s
Details
Deploy Documentation / build-and-deploy (push) Successful in 52s
Details
refactor: pipeline redesign & DIP compliance (Phase 4)
2026-06-04 08:47:29 +00:00

1676 lines
63 KiB
Python
Raw Blame History

"""
DNS Post-Processing Code Generator Core (post.py)
==================================================
이 모듈은 난류 및 연소 DNS(Direct Numerical Simulation) 데이터의 후처리를 위한 고성능 Fortran 코드를 생성하는 컴파일러의 코어입니다.
사용자가 정의한 DSL(Domain Specific Language) 입력 식을 읽어 파싱한 후, 다음과 같은 4단계 최적화 컴파일 과정을 거쳐 극도로 최적화된 3차원 루프 Fortran 모듈 코드를 자동 생성합니다.
[컴파일러 파이프라인 4단계 개요]
1. Stage 1 (AST 수집 및 변수 정의):
- Lark 파서가 생성한 AST(Abstract Syntax Tree)를 순회하며 기본 입력 변수(Primary), 계산이 필요한 대입식(Derived),
그리고 통계 물리량 계산을 위한 평균화 변수(Averaged)를 추출하고 필드 메타데이터 객체를 구성합니다.
2. Stage 2 (수치 미분 및 변동량 확장):
- 수식 내에 수치 미분자(ddx, d2dy 등)나 변동량(fluctuation, u')이 존재하면, 이를 물리적으로 차분 연산할 중간
미분 필드(DerivedField) 및 변동량 필드(FluctuationField)로 자동 변환하고 변수 테이블에 등록하여 확장합니다.
3. Stage 3 (의존성 분석 및 위상 정렬):
- 필드 간의 선후 연산 관계를 분석하여 유향 의존성 그래프(Directed Dependency Graph)를 생성합니다.
- 난류 통계의 특성상 평균 연산을 기준으로 "평균치 계산 전의 루프(Pass 1)"와 "평균치를 구한 후 변동량을 계산하는 루프(Pass 2)"로
전체 연산 블록을 논리적으로 분할하고, 각각의 블록 내에서 올바른 순서로 계산되도록 위상 정렬(Topological Sort)을 수행합니다.
4. Stage 4 (수식 기호 최적화, Liveness 분석 및 Buffer Array Pooling):
- SymPy 기호 수학 라이브러리를 이용하여 복잡한 3차원 수식을 대수적으로 간소화하고, 공통 부분 식 제거(CSE)를 적용하여 연산 비용(Flops)을 최적화합니다.
- 격자 데이터가 거대하므로 모든 변수에 개별 3D 배열을 할당하면 메모리가 고갈됩니다. 이를 방지하기 위해 변수들의
생명 주기(Liveness Window)를 수학적으로 추적하고, 동적 메모리 풀을 구축하여 동시에 활성화되지 않는 임시 변수들이
공통의 제한된 버퍼 배열(xyzbuffer0, xyzbuffer1, ...)을 나누어 사용(Array Pooling)하도록 할당하여 메모리 사용량을 최소화합니다.
"""
import sys
from lark import Lark, Visitor, Transformer, v_args, Token
import warnings
from jinja2 import Template
import sympy
from sympy.printing.fortran import FCodePrinter
# --- Registries for SOLID (OCP) Compliance ---
class FunctionRegistry:
def __init__(self):
self._sympy_registry = {}
self._latex_registry = {}
def register_sympy(self, name, sympy_builder):
self._sympy_registry[name] = sympy_builder
def register_latex(self, name, latex_builder):
self._latex_registry[name] = latex_builder
def to_sympy(self, name, *args):
if name in self._sympy_registry:
return self._sympy_registry[name](*args)
return sympy.Function(name)(*args)
def to_latex(self, name, *args):
if name in self._latex_registry:
return self._latex_registry[name](*args)
b = ", ".join(args)
if name.startswith("\\"):
return r"{}{{({})}}".format(name, b)
return r"\mathrm{{{}}}({})".format(name, b)
function_registry = FunctionRegistry()
# Register standard mathematical functions
function_registry.register_sympy("sqrt", lambda *args: sympy.sqrt(args[0]))
function_registry.register_sympy("abs", lambda *args: sympy.Abs(args[0]))
function_registry.register_sympy("log", lambda *args: sympy.log(args[0]))
function_registry.register_sympy("exp", lambda *args: sympy.exp(args[0]))
function_registry.register_sympy("rxn_rate", lambda *args: sympy.Function("rxn_rate")(args[0]))
function_registry.register_latex("sqrt", lambda *args: r"\sqrt{{{}}}".format(", ".join(args)))
function_registry.register_latex("abs", lambda *args: r"\left| {} \right|".format(", ".join(args)))
function_registry.register_latex(r"\log", lambda *args: r"\log{{({})}}".format(", ".join(args)))
function_registry.register_latex(r"\exp", lambda *args: r"\exp{{({})}}".format(", ".join(args)))
function_registry.register_latex(r"\omega", lambda *args: r"\omega{{({})}}".format(", ".join(args)))
class DifferentialOperatorRegistry:
def __init__(self):
self._operators = {}
def register(self, op_name, latex_symbol):
self._operators[op_name] = latex_symbol
def get_latex_symbol(self, op_name):
if op_name in self._operators:
return self._operators[op_name]
# Fallback to dynamic parsing matching original code
fmt = r"\partial_{{{}}}"
coord = op_name[-1] if op_name else ""
return fmt.format(coord + coord if len(op_name) > 3 else coord)
differential_operator_registry = DifferentialOperatorRegistry()
# Register standard derivative operators
differential_operator_registry.register("ddx", r"\partial_{x}")
differential_operator_registry.register("d2dx", r"\partial_{xx}")
differential_operator_registry.register("ddy", r"\partial_{y}")
differential_operator_registry.register("d2dy", r"\partial_{yy}")
differential_operator_registry.register("ddz", r"\partial_{z}")
differential_operator_registry.register("d2dz", r"\partial_{zz}")
class FortranTemplateStore:
REAL_ARRAY_LOOP = """
! {{ comment }}
{% if decls_str -%}
block
{{ decls_str | indent(4, True) }}
{%- endif %}
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
{% if assigns_str -%}
{{ assigns_str | indent(4, True) }}
{{ array }}(i,j,k) = {{ rhs }}
{%- else -%}
{{ array }}(i,j,k) = {{ rhs }}
{%- endif %}
end do
end do
end do
{% if decls_str -%}
end block
{%- endif %}
"""
FLUCTUATION_ARRAY_LOOP = """
! {{ comment }}
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
{{ array }}(i,j,k) = {{ rhs }}
end do
end do
end do
"""
AVG_ARRAY_SUM = """
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
{{ name }}(i) = {{ name }}(i) + {{ arrname }}
end do
end do
end do
"""
AVG_ARRAY_DIVIDE = """
call MPI_ALLREDUCE(MPI_IN_PLACE, {{ name }}, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
{{ name }} = {{ name }} {{ dWeight }} / denum
"""
FMT_DECL_SUBARRAY = """
! - file_handles and mpi_infos
integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_fh
integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_info
integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_filetype
"""
FMT_INIT_SUBARRAY = """
! init subarray datatype for {{ field_name }}
block
integer(4) :: sizes(3), subsizes(3), starts(3)
call MPI_INFO_CREATE({{ field_name }}_info, mpi_err)
call MPI_FILE_OPEN(MPI_COMM_TASK,'export-{{ field_name }}.dat',MPI_MODE_WRONLY+MPI_MODE_CREATE,{{ field_name }}_info,{{ field_name }}_fh,mpi_err)
sizes = (/ nxp, nyp, nzp /)
subsizes = (/ {{ len_xpts }}, {{ ye }} - {{ ys }} + 1, {{ ze }} - {{ zs }} + 1 /)
starts = (/ {{ xs }} - 1, {{ ys }} - 1, {{ zs }} - 1 /)
call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts, MPI_ORDER_FORTRAN, MPI_REAL8, {{ field_name }}_filetype, mpi_err)
call MPI_TYPE_COMMIT({{ field_name }}_filetype, mpi_err)
end block
"""
FMT_FINAL_SUBARRAY = """
! finalize
call MPI_FILE_CLOSE({{ field_name }}_fh, mpi_err)
call MPI_INFO_FREE({{ field_name }}_info, mpi_err)
call MPI_TYPE_FREE({{ field_name }}_filetype, mpi_err)
"""
FMT_CALC_SUBARRAY = """
! write to file via MPI Subarray
count = ({{ len_xpts }}) * ({{ ye }} - {{ ys }} + 1) * ({{ ze }} - {{ zs }} + 1)
offset = export_offset(fidx) * count * 8
call MPI_FILE_WRITE_AT({{ field_name }}_fh, offset, {{ work_array }}, 1, {{ field_name }}_filetype, mpi_status, mpi_err)
"""
FMT_DECL_LEGACY = """
! - file_handles and mpi_infos
integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_fh
integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_info
! - buffer
real(real64), allocatable, dimension(:,:,:) :: {{ field_name }}_export_array
integer, allocatable, dimension(:) :: {{ field_name }}_xpts
"""
FMT_INIT_LEGACY = """
! init
call MPI_INFO_CREATE({{ field_name }}_info, mpi_err)
call MPI_FILE_OPEN(MPI_COMM_TASK,'export-{{ field_name }}.dat',MPI_MODE_WRONLY+MPI_MODE_CREATE,{{ field_name }}_info,{{ field_name }}_fh,mpi_err)
allocate({{ field_name }}_export_array(1:{{ len_xpts }},{{ ys }}:{{ ye }},{{ zs }}:{{ ze }}), stat=ierr)
if (ierr /= 0) then
write(0,*) 'Error: allocation of {{ field_name }}_export_array failed on process', myid
call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
end if
{{ field_name }}_export_array = 0.
allocate({{ field_name }}_xpts(1:{{ len_xpts }}), stat=ierr)
if (ierr /= 0) then
write(0,*) 'Error: allocation of {{ field_name }}_xpts failed on process', myid
call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
end if
{{ xpts_init }}
"""
FMT_FINAL_LEGACY = """
! finalize
call MPI_FILE_CLOSE({{ field_name }}_fh, mpi_err)
call MPI_INFO_FREE({{ field_name }}_info, mpi_err)
deallocate({{ field_name }}_export_array)
deallocate({{ field_name }}_xpts)
"""
FMT_CALC_LEGACY = """
! copy to array for export
do k = {{ zs }}, {{ ze }}
do j = {{ ys }}, {{ ye }}
do i = 1, {{ len_xpts }}
{{ field_name }}_export_array(i,j,k) = {{ work_array }}({{ field_name }}_xpts(i),j,k)
end do
end do
end do
! write to file
count = ({{ len_xpts }}) * ({{ ye }} - {{ ys }} + 1) * ({{ ze }} - {{ zs }} + 1)
offset = export_offset(fidx) * count * 8
call MPI_FILE_WRITE_AT({{ field_name }}_fh, offset, {{ field_name }}_export_array, count, MPI_REAL8, mpi_status, mpi_err)
"""
AVG_ARRAY_WRITE = """
real(real64), dimension(nxp) :: xbuffer
integer :: i
open (200, file="qEdge_X.dat")
write (200,*) output_header
do i=1,nxp
write (200,'({{ num_args }}e20.10)') real(i)*hxp, {{ formatted_avglist }}
end do
close (200)
open (200, file="d1.dat")
{{ deriv1_lines }}
close (200)
open (200, file="d2.dat")
{{ deriv2_lines }}
close (200)
"""
class FortranCodeGenerator(object):
def __init__(self, fdict):
self.fdict = fdict
def generate_code(self, field, alloc=None):
method_name = 'visit_' + field.__class__.__name__ + '_code'
visitor = getattr(self, method_name, self.generic_code)
return visitor(field, alloc)
def generate_decl(self, field):
method_name = 'visit_' + field.__class__.__name__ + '_decl'
visitor = getattr(self, method_name, self.generic_decl)
return visitor(field)
def generate_alloc(self, field):
method_name = 'visit_' + field.__class__.__name__ + '_alloc'
visitor = getattr(self, method_name, self.generic_alloc)
return visitor(field)
def generate_free(self, field):
method_name = 'visit_' + field.__class__.__name__ + '_free'
visitor = getattr(self, method_name, self.generic_free)
return visitor(field)
def generate_avg(self, field):
method_name = 'visit_' + field.__class__.__name__ + '_avg'
visitor = getattr(self, method_name, self.generic_avg)
return visitor(field)
# --- Generic Fallbacks ---
def generic_code(self, field, alloc=None):
return ""
def generic_decl(self, field):
real_array_decl = "real(real64), allocatable, dimension({1}) :: {0}"
return real_array_decl.format(field.name, field.dim)
def generic_alloc(self, field):
return make_allocate(field.name, field.shape)
def generic_free(self, field):
real_array_free = "deallocate({})"
return real_array_free.format(field.name)
def generic_avg(self, field):
return ""
# --- Visit Methods ---
def visit_FieldExporter_code(self, exporter, alloc=None):
exporter.params["work_array"] = exporter.parent.array
if exporter.use_subarray:
return Template(FortranTemplateStore.FMT_CALC_SUBARRAY).render(**exporter.params)
else:
return Template(FortranTemplateStore.FMT_CALC_LEGACY).render(**exporter.params)
def visit_FieldExporter_decl(self, exporter):
if exporter.use_subarray:
return Template(FortranTemplateStore.FMT_DECL_SUBARRAY).render(**exporter.params)
else:
return Template(FortranTemplateStore.FMT_DECL_LEGACY).render(**exporter.params)
def visit_FieldExporter_alloc(self, exporter):
if exporter.use_subarray:
return Template(FortranTemplateStore.FMT_INIT_SUBARRAY).render(**exporter.params)
else:
return Template(FortranTemplateStore.FMT_INIT_LEGACY).render(**exporter.params)
def visit_FieldExporter_free(self, exporter):
if exporter.use_subarray:
return Template(FortranTemplateStore.FMT_FINAL_SUBARRAY).render(**exporter.params)
else:
return Template(FortranTemplateStore.FMT_FINAL_LEGACY).render(**exporter.params)
def visit_Field_code(self, field, alloc=None):
field.array = alloc[field.name] if alloc else field.name
opt = SympyOptimizer.get_instance(self.fdict)
rhs, cse_decls, cse_assigns = opt.optimize_field(field.name, alloc)
decls_str = "\n".join(cse_decls) if cse_decls else ""
assigns_str = "\n".join(cse_assigns) if cse_assigns else ""
calculation_code = Template(FortranTemplateStore.REAL_ARRAY_LOOP).render(
comment=field.comment,
decls_str=decls_str,
assigns_str=assigns_str,
array=field.array,
rhs=rhs
)
export_code = self.generate_code(field.exporter) if field.export_on() else ""
return calculation_code + export_code
def visit_FluctuationField_code(self, field, alloc=None):
field.array = alloc[field.name] if alloc else field.name
rhs = ExpToCode(self.fdict).transform(field.field.exp)
if field.field.is_fluctuation():
rhs = rhs.format(field.w)
return Template(FortranTemplateStore.FLUCTUATION_ARRAY_LOOP).render(
comment=field.comment,
array=field.array,
rhs=rhs
)
def visit_PrimaryField_code(self, field, alloc=None):
return "! {} is read from file".format(field.name)
def visit_PrimaryField_decl(self, field):
return "! {} is read from file".format(field.name)
def visit_PrimaryField_alloc(self, field):
return "! {} is read from file".format(field.name)
def visit_PrimaryField_free(self, field):
return "! {} is read from file".format(field.name)
def visit_DerivedField_code(self, field, alloc=None):
field.array = alloc[field.name] if alloc else field.name
varray = alloc[field.v] if alloc else field.v
return "call {0} ( {2}, {1} )".format(field.op, varray, field.array)
def visit_AveragedField_code(self, field, alloc=None):
arrname = self.fdict[field.tgt].array + "(i,j,k)"
if field.weighted is not None:
arrname = arrname + " * " + field.w.array + "(i,j,k)"
return Template(FortranTemplateStore.AVG_ARRAY_SUM).render(name=field.name, arrname=arrname)
def visit_AveragedField_avg(self, field):
dWeight = (f"/ avg_{field.weighted}" if field.weighted else "")
return Template(FortranTemplateStore.AVG_ARRAY_DIVIDE).render(name=field.name, dWeight=dWeight)
def generate_write_avg(self, avglist):
avgarr = "{}(i)"
deriv1_avgarr = """call ddx1d ( xbuffer, {} ) ; write (200,*) xbuffer"""
deriv2_avgarr = """call d2dx1d ( xbuffer, {} ) ; write (200,*) xbuffer"""
num_args = len(avglist) + 1
formatted_avglist = ", ".join(map(avgarr.format, avglist))
deriv1_lines = "\n".join(map(deriv1_avgarr.format, avglist))
deriv2_lines = "\n".join(map(deriv2_avgarr.format, avglist))
write_avg = Template(FortranTemplateStore.AVG_ARRAY_WRITE).render(
num_args=num_args,
formatted_avglist=formatted_avglist,
deriv1_lines=deriv1_lines,
deriv2_lines=deriv2_lines
)
return write_avg
@v_args(inline=True)
class LarkToSympy(Transformer):
"""Lark AST의 수학적 노드들을 SymPy 기호 수식 객체로 변환하는 Transformer 클래스입니다.
이를 통해 DSL 수식 트리를 SymPy 기호 객체로 매핑하고, SymPy의 강력한 수식 간소화(Simplify)
및 공통 부분식 제거(CSE) 최적화 도구를 파이프라인 하부에서 활용할 수 있게 합니다.
"""
def __init__(self, fdict):
"""Lark-to-SymPy Transformer를 초기화합니다.
Args:
fdict (dict): 변수명에서 물리 필드 정의 객체(FieldBase 자식 클래스들)로 매핑되는 딕셔너리.
"""
self.fdict = fdict
def number(self, numeral):
# 상수를 SymPy Float 기호 객체로 변환
return sympy.Float(float(numeral))
def env(self, name):
# 환경 변수($ 기호로 시작)를 SymPy Symbol로 매핑
return sympy.Symbol(name.value)
def paren(self, val):
# 괄호는 수식 우선순위 보존 후 내부 식 반환
return val
def var(self, name):
# 일반 변수명을 SymPy Symbol로 매핑
return sympy.Symbol(name.value)
def fluc(self, name):
# 변동량 u' 기호를 SymPy가 해석할 수 있도록 u__prime 기호명으로 매핑
return sympy.Symbol(name.value + "__prime")
def dnx(self, partial, b):
# ddx(u)와 같은 수치 미분항을 ddx_u 형태의 하나의 고유 Symbol로 묶어서 취급
signature = f"{partial.data}_{b.value}"
return sympy.Symbol(signature)
def icall(self, op, val):
# sqr, pow3 등의 인라인 전용 함수들을 SymPy 지수 표현식으로 다이렉트 변환
if op.data == "sqr":
return val**2
elif op.data == "pow3":
return val**3
return val
def fcall(self, *args):
# 내장 수학 함수(sqrt, exp, log, abs, rxn_rate) 또는 사용자 정의 UDF들을 SymPy 함수 노드로 매핑
a = args[0]
func_name = a.value if hasattr(a, 'value') else str(a)
if func_name == "udf":
return sympy.Function(a.value if hasattr(a, 'value') else str(a))(*args[1:])
return function_registry.to_sympy(func_name, *args[1:])
def neg(self, val):
return -val
def add(self, a, b):
return a + b
def sub(self, a, b):
return a - b
def mul(self, a, b):
return a * b
def div(self, a, b):
return a / b
def udf(self, a):
return a.value
log = lambda self: "log"
exp = lambda self: "exp"
sqrt = lambda self: "sqrt"
abs = lambda self: "abs"
rxn_rate = lambda self: "rxn_rate"
class ArrayFCodePrinter(FCodePrinter):
"""SymPy 식을 Fortran 코드로 인쇄할 때, 변수들을 적절한 다차원 격자점 인덱스 배열(예: var(i,j,k))
또는 평균 배열(예: var(i)) 형식으로 포맷팅하여 변환 출력해 주는 특수 Printer 클래스입니다.
"""
def __init__(self, settings=None, array_symbols=None, avg_symbols=None):
settings = settings or {}
settings.setdefault('source_format', 'free') # 기본 프리포맷 Fortran 95 준수
settings.setdefault('standard', 95)
super().__init__(settings)
self.array_symbols = array_symbols or {} # 3차원 필드 배열의 실제 물리 버퍼 명칭 매핑
self.avg_symbols = avg_symbols or {} # 평균 필드 배열의 명칭 매핑
def _print_Float(self, expr):
# 배정밀도 실수 상수에 d0 접미사를 명시하여 컴파일러 정밀도 유실 방지
val = str(expr)
if 'e' in val or 'E' in val:
return val.replace('e', 'd').replace('E', 'd')
if '.' not in val:
return val + ".0d0"
return val + "d0"
def _print_Symbol(self, expr):
name = expr.name
# 1. 3D 공간 격자점 버퍼 할당명에 매핑되어 있는 경우 (i,j,k) 차분 격자점 인덱스를 부착
if name in self.array_symbols:
return f"{self.array_symbols[name]}(i,j,k)"
# 2. X방향(i) 1차원 평균화 변수인 경우 (i) 인덱스 부착
if name in self.avg_symbols:
return f"{self.avg_symbols[name]}(i)"
# 3. 변동량 변수(u__prime)인 경우, 원본 격자점 값과 1차원 평균값의 편차식인 '(u(i,j,k) - avg_u(i))'로 자동 인라인 전개
if name.endswith("__prime"):
base = name[:-7]
arr = self.array_symbols.get(base, base)
avg_name = f"avg_{base}"
printed_avg = f"{self.avg_symbols.get(avg_name, avg_name)}(i)"
return f"({arr}(i,j,k) - {printed_avg})"
return name
def _print_Function(self, expr):
# SymPy 내장 함수가 아닌 사용자 정의 함수(rxn_rate, udf 등)의 정상적인 Fortran 출력을 위한 예외 폴백 처리
try:
return super()._print_Function(expr)
except Exception:
args = ", ".join(self.doprint(arg) for arg in expr.args)
return f"{expr.func.__name__}({args})"
class SympyOptimizer:
"""컴파일러 내부에서 SymPy 수학 기호 연산 엔진을 제어하는 최적화 매니저 클래스입니다.
1. 수식 확장 대입(Substitution): 사용자가 정의한 물리 변수 중 디스크로 내보내지(export) 않는
임시 변수들의 수식을 참조하는 하부 수식에 재귀적으로 대입하여 거대한 하나의 식으로 확장합니다.
2. 대수적 간소화(Simplify & Cancel): SymPy를 사용하여 분모 분자 소거 및 삼각/지수 대수 간소화를 적용합니다.
3. 공통 부분식 제거(CSE - Common Subexpression Elimination): 복잡하게 얽힌 다차원 미분 수식 속에서
중복 계산되는 서브 식(예: 공통 항곱)을 탐지하여, 루프 외/내부에서 임시 스칼라 변수(x0, x1 등)에 선언한 후
루프 내 연산에서 한 번만 대입하여 연산하도록 루프 코드를 개조하여 FLOPS와 캐시 로드를 획기적으로 낮춥니다.
"""
_instance = None
@classmethod
def get_instance(cls, fdict):
if cls._instance is None or cls._instance.fdict is not fdict:
cls._instance = cls(fdict)
return cls._instance
def __init__(self, fdict):
self.fdict = fdict
self.sympy_cache = {} # 중복 연산 계산 방지를 위한 SymPy 식 캐시
self.exported_fields = set(
name for name, f in fdict.items()
if hasattr(f, 'attr') and f.attr.get('export')
) # 디스크로 파일 출력이 설정된 최종 목표 변수 집합
self.averaged_targets = set()
self.avg_names = set()
def set_averaged(self, averaged_dict):
# 평균화 대상이 되는 변수 집합을 수집
self.averaged_targets = {a.target for a in averaged_dict.values()}
self.avg_names = set(averaged_dict.keys())
def get_sympy_expr(self, name):
"""특정 변수 명칭에 대해 기호적으로 완전히 확장 대입된 SymPy Expression 객체를 빌드하여 캐싱합니다.
기본 입력 변수(Primary), 미분 변수(Derived), 평균 변수(Averaged), 변동량 변수(Fluctuation)는
치환하지 않고 독립적인 리프 기호(Symbol)로 둡니다. 그 외 중간에서만 활용되는 일반 임시 Derived 변수들은
의존하는 원본 수식들을 탐색하여 재귀적으로 인라인 대체(expand) 시킴으로써, SymPy 기호 엔진이
루프 수식 전체의 논리적 최적화를 전역적으로 수행할 수 있도록 합니다.
"""
if name in self.sympy_cache:
return self.sympy_cache[name]
field = self.fdict[name]
if hasattr(field, 'prime') and field.prime:
expr = sympy.Symbol(name)
self.sympy_cache[name] = expr
return expr
if hasattr(field, 'op'): # DerivedField (ddx 등 미분 결과)
expr = sympy.Symbol(name)
self.sympy_cache[name] = expr
return expr
if hasattr(field, 'weighted'): # AveragedField
expr = sympy.Symbol(name)
self.sympy_cache[name] = expr
return expr
if hasattr(field, 'field') and hasattr(field, 'w'): # FluctuationField
expr = sympy.Symbol(name)
self.sympy_cache[name] = expr
return expr
transformer = LarkToSympy(self.fdict)
expr = transformer.transform(field.exp)
# 재귀적으로 의존 중인 캐싱되지 않은 중간 계산 임시 변수들을 수식 본문으로 확장 치환(Substitution)
expanded_expr = expr
changed = True
while changed:
changed = False
free_syms = list(expanded_expr.free_symbols)
sub_dict = {}
for sym in free_syms:
sym_name = sym.name
if sym_name in self.fdict:
f = self.fdict[sym_name]
is_derived_field = hasattr(f, 'op')
is_averaged_field = hasattr(f, 'weighted')
is_primary_field = hasattr(f, 'prime') and f.prime
is_exported = sym_name in self.exported_fields
is_averaged_target = sym_name in self.averaged_targets
# 최종 출력 대상이거나 물리적 경계 처리가 필요한 핵심 변수가 아니라면,
# 기 기하학적 수식을 최적화하기 위해 현재 식 안으로 대입해 버립니다.
if not (is_derived_field or is_averaged_field or is_primary_field or is_exported or is_averaged_target):
sub_dict[sym] = self.get_sympy_expr(sym_name)
changed = True
if sub_dict:
expanded_expr = expanded_expr.subs(sub_dict)
self.sympy_cache[name] = expanded_expr
return expanded_expr
def calculate_flops_and_heavy(self, expr):
"""수식에서 수행되는 부동 소수점 연산(FLOPS) 및 계산 부하가 큰 연산(Division, Sqrt, Exp 등)의 개수를 평가합니다.
이를 통해 최적화 리포트에 연산 비용 증감률을 리포팅합니다.
"""
flops = 0
heavy = 0
for node in sympy.preorder_traversal(expr):
if isinstance(node, sympy.Add):
flops += len(node.args) - 1
elif isinstance(node, sympy.Mul):
flops += len(node.args) - 1
elif isinstance(node, sympy.Pow):
base, exp = node.args
if exp == 0.5 or exp == -0.5:
flops += 10 # 제곱근 및 역제곱근은 대략 10 flops 가중치 부여
heavy += 1
elif exp == -1:
flops += 4 # 역수 나눗셈은 4 flops 가중치 부여
heavy += 1
elif isinstance(exp, sympy.Integer):
val = abs(int(exp))
if val > 1:
flops += val - 1
else:
flops += 10
heavy += 1
elif isinstance(node, (sympy.Derivative, sympy.Function)):
name = node.func.__name__
if name == 'sqrt':
flops += 10
heavy += 1
elif name in ('exp', 'log', 'sin', 'cos', 'tan', 'rxn_rate'):
flops += 10 # 특수 수학 함수 연산은 10 flops 가중치 및 Heavy 연산자로 평가
heavy += 1
elif name == 'Abs':
flops += 1
else:
flops += 10
heavy += 1
return flops, heavy
def count_3d_loads(self, expr, three_d_arrays):
"""수식 내부에서 참조하는 3차원 격자점 배열의 총 메모리 로드 횟수를 카운트합니다.
메모리 대역폭 한계(Memory-bound)에 부딪히는 HPC 연산의 병목을 평가하는 데 중요한 지표입니다.
"""
count = 0
for node in sympy.preorder_traversal(expr):
if isinstance(node, sympy.Symbol) and node.name in three_d_arrays:
count += 1
return count
def optimize_field(self, name, alloc=None):
"""개별 물리 필드 수식을 SymPy를 사용해 최적화하고 공통 부분식(CSE) 코드를 추출합니다.
Args:
name (str): 최적화할 필드 변수명.
alloc (dict): 변수명에서 버퍼 배열명(xyzbufferN)으로 매핑되는 테이블.
Returns:
rhs (str): 최적화가 완료된 최종 Fortran 우변(Right-Hand Side) 코드 문자열.
cse_decls (list): CSE 추출 결과 생성된 임시 실수 변수 선언 리스트 (예: real(real64) :: x0).
cse_assigns (list): 루프 내부 대입에 사용될 임시 변수의 연산식 리스트.
"""
expr = self.get_sympy_expr(name)
three_d_arrays = {
k for k, v in self.fdict.items()
if hasattr(v, 'dim') and v.dim == ':,:,:'
}
# 최적화 전 비용 수집
before_flops, before_heavy = self.calculate_flops_and_heavy(expr)
before_loads = self.count_3d_loads(expr, three_d_arrays)
# 1차 대수적 간소화 및 약분 소거
simplified_expr = sympy.simplify(expr)
simplified_expr = sympy.cancel(simplified_expr)
# 기호명을 실제 Fortran 격자 배열명(xyzbufferN)으로 매핑하기 위한 준비
array_symbols = {}
for k, v in self.fdict.items():
if hasattr(v, 'array') and v.array:
array_symbols[k] = v.array
elif alloc and k in alloc:
array_symbols[k] = alloc[k]
else:
array_symbols[k] = k
avg_symbols = {k: k for k in getattr(self, 'avg_names', [])}
printer = ArrayFCodePrinter(array_symbols=array_symbols, avg_symbols=avg_symbols)
# 2차: 공통 부분식 제거(CSE) 알고리즘 가동
# 복잡하게 중복 사용되는 수식 구조를 서브 트리로 분리하여 임시 기호(x0, x1 등)와 최종 축소 식으로 분할
replacements, reduced_exprs = sympy.cse(simplified_expr)
reduced_expr = reduced_exprs[0]
# 최적화 후 비용 계산
after_flops = 0
after_heavy = 0
after_loads = 0
for temp_var, temp_expr in replacements:
f_val, h_val = self.calculate_flops_and_heavy(temp_expr)
after_flops += f_val
after_heavy += h_val
after_loads += self.count_3d_loads(temp_expr, three_d_arrays)
f_val, h_val = self.calculate_flops_and_heavy(reduced_expr)
after_flops += f_val
after_heavy += h_val
after_loads += self.count_3d_loads(reduced_expr, three_d_arrays)
def pct_str(before, after):
if before == 0:
return "0.0%" if after == 0 else "+inf%"
diff = after - before
pct = (diff / before) * 100
return f"{pct:+.1f}%"
flops_pct = pct_str(before_flops, after_flops)
heavy_pct = pct_str(before_heavy, after_heavy)
loads_pct = pct_str(before_loads, after_loads)
# 성능 개선 등급 진단
if after_flops < before_flops * 0.5 or after_loads < before_loads * 0.5:
est_speedup = "Highly significant"
elif after_flops < before_flops or after_loads < before_loads:
est_speedup = "Moderate"
else:
est_speedup = "Minimal / Already optimal"
# 컴파일 중 최적화 상세 리포트를 stderr로 콘솔 화면에 로깅
sys.stderr.write(f"\n[SymPy Optimizer Report: {name}]\n")
sys.stderr.write(f"- Floating Point Ops : {before_flops} -> {after_flops} ({flops_pct})\n")
sys.stderr.write(f"- Heavy Ops (Div/Sqrt): {before_heavy} -> {after_heavy} ({heavy_pct})\n")
sys.stderr.write(f"- 3D Array Mem Reads : {before_loads} -> {after_loads} ({loads_pct})\n")
sys.stderr.write(f"=> Estimated Speedup in loop: {est_speedup}\n\n")
cse_decls = []
cse_assigns = []
# CSE 추출 결과를 Fortran 코드로 포매팅
if replacements:
for temp_var, temp_expr in replacements:
# 임시 변수 로컬 스칼라 선언 구문 작성
cse_decls.append(f"real(real64) :: {temp_var}")
# 루프 내에서 먼저 연산 후 할당하는 구문 작성
cse_assigns.append(f"{temp_var} = {printer.doprint(temp_expr)}")
rhs = printer.doprint(reduced_expr)
return rhs, cse_decls, cse_assigns
class CollectDefinitions(Visitor):
"""Lark AST의 전체 노드를 스캔하며 정의문(변수 리스트, 대입식, 평균화 조건)을 수집하여
컴파일러의 빌딩 블록 객체들로 구축하는 Visitor 클래스입니다 (Stage 1 핵심 작동부).
"""
def __init__ (self, primary, derived, averaged):
self.primary = primary # 기본 입력 물리 필드명 집합 (격자 파일 등에서 읽음)
self.derived = derived # 유도/계산식 필드 딕셔너리
self.averaged = averaged # 평균 연산 조건 딕셔너리
def varlist(self, tree):
# 괄호 안 [u, v, w] 같은 형식으로 표기된 기본(Primary) 변수 리스트 파싱
for v in tree.children:
self.primary.add(v.value)
# 의존성 테이블에 기본 변수 객체 등록
self.derived[v.value] = PrimaryField(v.value, self.derived)
def assign_var (self, tree):
# 등호(=) 연산자를 통한 변수 계산식 정의 노드 파싱
# 속성값(예: (export=true, latex="..."))이 존재하면 추출
if len(tree.children) > 2:
lval, lattr, rval = tree.children
else:
lval, rval = tree.children
lattr = None
attr_dict = {}
if lattr is not None:
for t in lattr.children:
k, v = t.children
attr_dict[k.value] = v.value
if lval.value in self.derived:
raise ValueError("duplicate definition of " + lval)
# Field 객체를 생성하여 derived 계산 테이블에 등록
self.derived[lval.value] = Field(lval.value, attr_dict, rval, self.derived)
def assign_avg_var (self, tree):
# avg 또는 avg$w {u, v, ...} 와 같은 평균화 지시어 파싱 및 등록
w = tree.children[0]
targets = tree.children[1:]
if (not w.children) or (w.children[0] is None):
self.averaged[None] = set([x.value for x in targets])
else:
self.averaged[w.children[0].value] = set([x.value for x in targets])
class ExpInspector(Visitor):
"""수식 AST 노드들을 순회하며 해당 식 내부의 변동량(Fluctuation) 포함 여부,
참조 중인 다른 물리 변수 종속성 목록(dep), 그리고 수치 미분항(deriv) 목록을 추출하는 Visitor 클래스입니다.
"""
def __init__(self):
self.fluctuation = False
self.dep = set([])
self.deriv = set([])
@classmethod
def inspect(cls, tree):
self = cls()
return self(tree)
def __call__(self, tree):
self.visit(tree)
return self.fluctuation, self.dep, self.deriv
def fluc(self, tree):
# 식 내부에 변동량 항이 존재함을 기록하고 종속성에 등록
self.fluctuation = True
self.dep.add(tree.children[0].value)
def var(self, tree):
self.dep.add(tree.children[0].value)
def dnx (self, tree):
# ddx(u) 와 같이 수치 미분 노드를 만나면 ddx_u 형태의 중간 미분 변수를 종속성 리스트에 기록
op, v = tree.children
deriv = "{}_{}".format(op.data, v.value)
self.dep.add(deriv)
self.deriv.add((op.data, v.value))
@v_args(inline=True) # 메서드의 서명 형식을 자식 노드 매개변수 나열형으로 변환
class ExpToLatex(Transformer):
"""DSL 수학 표현식 AST를 문서를 검토하기 편하게 LaTeX 문법 수식으로 번역하는 번역기입니다."""
def __init__(self, fdict):
self.fdict = fdict
def arithmatic_rooted(self, name):
try:
exproot = self.fdict[name].exp.data
except AttributeError:
exproot = "something_11fasq2afa3rfzsaerqw23"
return ((exproot == "add") or (exproot == "sub") or
(exproot == "mul") or (exproot == "div"))
def parenthise(self, name):
try:
latex = self.fdict[name].latex
latex_given = self.fdict[name].latex_given
except KeyError:
warnings.warn(name + " is not found")
latex = r"\mathrm{{{}}}".format(name)
latex_given = None
if self.arithmatic_rooted(name) and (latex_given is None):
latex = "(" + latex + ")"
return latex
def number(self, numeral):
return numeral
def env(self, name):
return r"\mathrm{{{}}}".format(name.value)
def paren(self, name):
return "({})".format(str(name))
def var(self, name):
return self.parenthise(name.value)
def fluc(self, name):
# LaTeX 수식 표기용 u'' 프라임 기호 렌더링
return self.parenthise(name.value)+ "''"
def dnx (self, partial, b):
fmt = r"\partial_{{{}}}"
coord = partial.data[-1]
op = fmt.format(coord + coord if len(partial.data) > 3 else coord)
signature = "{}_{}".format(partial.data, b.value)
try:
eq = self.fdict[signature].latex
except KeyError:
eq = op + self.parenthise(b.value)
warnings.warn(signature + " is not found: " + eq)
return eq
def icall (self, a, b):
if a.data == "sqr":
fcode = "({0})^2".format(b)
elif a.data == "pow3":
fcode = "({0})^3".format(b)
else:
fcode = "({0})".format(b)
return fcode
def fcall (self, *args):
a = args[0]
func_name = a.value if hasattr(a, 'value') else str(a)
return function_registry.to_latex(func_name, *args[1:])
def neg(self, b):
fcode = "(-{})".format(b)
return fcode
def add(self, a, b):
fcode = "{} + {}".format(a, b)
return fcode
def sub(self, a, b):
fcode = "{} - {}".format(a, b)
return fcode
def mul(self, a, b):
fcode = "{} {}".format(a, b)
return fcode
def div(self, a, b):
fcode = "{} / {}".format(a, b)
return fcode
log = lambda self : "\log"
exp = lambda self : "\exp"
sqrt = lambda self : "sqrt"
abs = lambda self : "abs"
rxn_rate = lambda self : "\omega"
udf = lambda self, a : a.value
@v_args(inline=True)
class ExpToCode(Transformer):
"""DSL 수학 수식 트리를 Fortran의 3차원 인덱스 코드 형태의 수식 문자열로 다이렉트 매핑하는 번역기입니다.
SymPy 최적화를 사용하지 않는 대체 연산 경로 등에서 보조적으로 사용됩니다.
"""
def __init__(self, fdict):
self.fdict = fdict
def number(self, numeral):
return str(float(numeral))
def env(self, name):
return name.value
def paren(self, name):
return "({})".format(str(name))
def var(self, name):
try:
arrname = self.fdict[name.value].array
except KeyError:
arrname = name.value
return arrname + "(i,j,k)"
def fluc(self, name):
try:
arrname = self.fdict[name.value].array
except KeyError:
arrname = name.value
fmt = "({0}(i,j,k) - {{0}}avg_{1}(i))"
return fmt.format(arrname, name.value)
def dnx (self, partial, b):
signature = "{}_{}".format(partial.data, b.value)
try:
arrname = self.fdict[signature].array
except KeyError:
arrname = signature
return arrname + "(i,j,k)"
def icall (self, a, b):
if a.data == "sqr":
fcode = "(({0})*({0}))".format(b)
elif a.data == "pow3":
fcode = "(({0})*({0})*({0}))".format(b)
else:
fcode = "({0})".format(b)
return fcode
def fcall (self, *args):
a = args[0]
b = ", ".join(args[1:])
fcode = "( {} ( {} ) )".format(a, b)
return fcode
def neg(self, b):
fcode = "( - {} )".format(b)
return fcode
def add(self, a, b):
fcode = "( {} + {} )".format(a, b)
return fcode
def sub(self, a, b):
fcode = "( {} - {} )".format(a, b)
return fcode
def mul(self, a, b):
fcode = "( {} * {} )".format(a, b)
return fcode
def div(self, a, b):
fcode = "( {} / {} )".format(a, b)
return fcode
log = lambda self : "log"
exp = lambda self : "exp"
sqrt = lambda self : "sqrt"
abs = lambda self : "abs"
rxn_rate = lambda self : "rxn_rate"
udf = lambda self, a : a.value
def make_allocate(name, shape, init_zero=True):
"""Fortran 배열을 동적 할당하고 예외 발생 시 프로세스를 안전하게 폭파시키는 할당 코드를 작성해 줍니다."""
alloc_str = f"allocate({name}({shape}), stat=ierr)\n"
alloc_str += f"if (ierr /= 0) then\n"
alloc_str += f" write(0,*) 'Error: allocation of {name} failed on process', myid\n"
alloc_str += f" call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)\n"
alloc_str += f"end if"
if init_zero:
alloc_str += f"\n{name} = 0."
return alloc_str
class DependencyNode(object):
"""의존성 관계 분석을 담당하는 추상 인터페이스 및 노드 클래스입니다."""
def __init__(self, name, fdict):
self.name = name
self.fdict = fdict
self.dep = set([]) # 종속된(의존 중인) 하위 필드 리스트
self.fluc = False # 변동량 편차 계산 대상인지 여부
def depends_on(self, a):
return a in self.dep
def is_fluctuation(self):
return self.fluc
def checkFluctuation(self):
"""본 변수 혹은 의존하고 있는 하위 기호들 중에 변동량(Fluctuation) 관련 계산이 개입되어 있는지
상향식으로 전파 추적하는 재귀 메서드입니다.
"""
fset = set([])
for d in map(self.fdict.get, self.dep):
fset.update(d.checkFluctuation())
if self.is_fluctuation() or len(fset) > 0:
fset.add(self.name)
return fset
def depClosure(self):
"""해당 변수를 계산하기 위해 선행 계산되어야 하는 모든 하위 변수 노드들을
재귀적으로 타고 내려가 총합 폐쇄 집합(Closure Set)으로 묶어 반환합니다.
"""
fset = set(self.dep)
for d in self.dep:
fset.update(self.fdict[d].depClosure())
return fset
class FieldBase (DependencyNode):
"""모든 물리 필드 객체의 최상위 기본 클래스로서 공통 데이터 구조를 정의합니다."""
def __init__ (self, name, fdict):
super(FieldBase, self).__init__(name, fdict)
self.array = name # 실제 메모리 버퍼 상의 지칭 배열 명칭
self.prime = False # 기본 데이터 파일로부터 읽어 들이는 원본 필드인지 여부
self.shape = "nxp,nyp,nzp" # 3차원 그리드 구조
self.dim = ":,:,:"
def export_on (self):
return False
def __repr__ (self):
return self.name
class FieldExporter (object):
"""물리 필드 데이터를 병렬 분산 디스크 시스템으로 직접 추출(Export)하는 고성능 MPI-IO 서브루틴 블록을
정의하는 도메인 데이터 클래스입니다.
"""
mpi_io_decl = """
! field exporter common
integer(kind=MPI_OFFSET_KIND) :: offset
"""
def __init__ (self, name, attr, parent):
self.name = name
self.attr = attr
self.parent = parent
self.params = dict(attr)
self.params.setdefault("xs", 1)
self.params.setdefault("xe", "nxp")
self.params.setdefault("ys", 1)
self.params.setdefault("ye", "nyp")
self.params.setdefault("zs", 1)
self.params.setdefault("ze", "nzp")
self.params.setdefault("field_name", self.name)
self.params.setdefault("len_xpts", f"({self.params['xe']} - {self.params['xs']} + 1)")
fmt_xpts_init = '''
do i = {{ xs }}, {{ xe }}
{{ field_name }}_xpts(i-{{ xs }}+1) = i
end do
'''
self.params.setdefault("xpts_init", Template(fmt_xpts_init).render(**self.params))
try:
# Sampling at listed x coordinates
fmt_xpts_init_list = "{{ field_name }}_xpts = (/ {{ list_xpts }} /)"
raw_xpts = self.params["xpts"]
int_xpts = list(map(int, raw_xpts.split()))
len_xpts = len(int_xpts)
self.params["len_xpts"] = len_xpts
self.params["list_xpts"] = ",".join(map(str, int_xpts))
self.params["xpts_init"] = Template(fmt_xpts_init_list).render(**self.params)
except KeyError:
pass
self.use_subarray = ("xpts" not in self.params)
class Field (FieldBase):
"""일반 대입 계산 변수를 관리하며, 3차원 격자점 루프 코드 생성을 담당하는 핵심 클래스입니다.
SymPy 최적화 및 CSE 적용 코드를 루프 본문에 결합합니다.
"""
def __init__ (self, name, attr, exp, fdict):
super(Field,self).__init__(name, fdict)
self.attr = attr
self.exp = exp
self.fluc, self.dep, self.derivs = ExpInspector.inspect(exp)
self.comment = self.name + " = " + ExpToCode(self.fdict).transform(self.exp)
self.latex_given = self.attr.get("latex")
if self.latex_given is None:
self.latex = ExpToLatex(self.fdict).transform(self.exp)
else:
self.latex = self.latex_given
self.exporter = None
try:
if self.attr["export"]:
self.exporter = FieldExporter(self.name, self.attr, self)
except KeyError:
pass
def export_on (self):
return self.exporter is not None
class FluctuationField (FieldBase):
"""물리 필드의 난류 변동 성분(Fluctuation, u' = u - <u_w>)을 계산하기 위한 변수 클래스입니다.
수식 내의 u' 기호를 평균량과의 차이 수식으로 팽창하여 할당합니다.
"""
def __init__ (self, w, field, fset, fdict):
super(FluctuationField,self).__init__(self.id(w,field), fdict)
if w is not None:
self.w = w + "_"
else:
self.w = ""
self.field = fdict[field]
self.dep = self.field.dep - fset
for df in self.field.dep & fset:
self.dep.add(self.id(w,df))
self.comment = ExpToCode(self.fdict).transform(self.field.exp)
if self.field.is_fluctuation():
self.comment = self.comment.format(self.w)
self.comment = self.name + " = " + self.comment
@classmethod
def id (cls, w, field):
if w:
name = "{}____{}_avg".format(field, w)
else:
name = "{}____avg".format(field)
return name
class PrimaryField (FieldBase):
"""격자 정보나 외부 물리계 수치 데이터(u, v, w, T 등) 파일에서 사전에 로드하여
메모리에 상주하는 기본 원본 입력 필드 클래스입니다.
자체 계산 루프나 동적 할당 코드를 직접 생성하지 않습니다.
"""
def __init__ (self, name, fdict):
super(PrimaryField,self).__init__(name, fdict)
self.derivs = set([])
self.prime = True
self.latex = name
self.latex_given = None
class DerivedField (FieldBase):
"""수치 공간 미분(ddx, d2dy 등)을 수행하여 계산되는 유도 필드 클래스입니다.
Fortran 수치 차분 패키지 서브루틴(Compact.f90 에 구현된 dfnonp, dfp 등)의
동적 호출 코드를 출력합니다.
"""
def __init__ (self, op, v, fdict):
name = "{}_{}".format(op, v)
super(DerivedField,self).__init__(name, fdict)
self.op = op
self.v = v
self.dep = set([v])
partial = differential_operator_registry.get_latex_symbol(op)
self.latex = partial + "(" + fdict[v].latex + ")"
class AveragedField (FieldBase):
"""특정 물리 필드를 격자의 동질 차원(예: X 방향 1D선상 평균)에 대해
공간 통계 평균(Average) 연산을 수행하는 1차원 필드 클래스입니다.
"""
@classmethod
def id (cls, w, tgt):
if w:
return "{}_avg_{}".format(w, tgt)
else:
return "avg_{}".format(tgt)
def __init__ (self, w, tgt, fdict):
name = self.id(w,tgt)
super(AveragedField,self).__init__(name, fdict)
self.shape = "nxp" # Y, Z 차원을 평균화하여 날려버리므로 X방향 크기인 nxp 1차원 배열이 됨
self.dim = ":"
self.target = tgt
tfield = fdict[tgt]
self.fset = tfield.checkFluctuation()
self.latex = r"\left\langle {} \right\rangle".format(tfield.latex)
if not self.fset:
self.tgt = tgt
self.dep.add(tgt)
else:
ftgt = FluctuationField.id(w,tgt)
self.tgt = ftgt
self.dep.add(ftgt)
self.weighted = w
if w:
self.w = fdict[w]
self.dep.add(w)
self.latex += ("_{{{}}}".format(w))
def isWeighted (self):
return self.weighted is not None
def pass1 (self):
return not self.pass2()
def pass2 (self):
# 수식 내에 u' 등 변동량 항이 포함되어 있으면,
# 원본 u의 평균치(<u_w>) 계산이 완료된 '이후'에만 연산이 가능하므로 Pass 2 대상이 됨
return len(self.fset) > 0
class CompilationContext(object):
def __init__(self):
self.primary = set()
self.derived = {}
self.averaged = {}
self.dependency = {}
self.pass1 = []
self.pass2 = []
self.avg1 = set()
self.avg2 = set()
self.alloc1 = {}
self.alloc2 = {}
self.narr = 0
calc_grammar = """
?varlist: "[" [NAME ("," NAME)*] "]"
?start: statement*
?statement: NAME [ attr_list ] "=" sum -> assign_var
| avg "{" [NAME ("," NAME)*] "}" -> assign_avg_var
| varlist
attr_list: "(" [attr_pair ("," attr_pair)*] ")"
?attr_pair: NAME "=" BOOL
| NAME "=" INT
| NAME "=" ESCAPED_STRING
?sum: product
| sum "+" product -> add
| sum "-" product -> sub
?product: atom
| product "*" atom -> mul
| product "/" atom -> div
?atom: NUMBER -> number
| "-" atom -> neg
| NAME -> var
| NAME "'" -> fluc
| "$" NAME -> env
| "(" sum ")" -> paren
| inlinefunc "(" sum ")" -> icall
| mathfunc "(" sum ("," sum)* ")" -> fcall
| derivative "(" NAME ")" -> dnx
avg: "avg" [NAME]
?inlinefunc: "sqr" -> sqr
| "pow3" -> pow3
?mathfunc: "log" -> log
| "exp" -> exp
| "sqrt" -> sqrt
| "abs" -> abs
| "rxn_rate" -> rxn_rate
| "$" NAME -> udf
?derivative: "ddx" -> ddx
| "d2dx" -> d2dx
| "ddy" -> ddy
| "d2dy" -> d2dy
| "ddz" -> ddz
| "d2dz" -> d2dz
%import common.CNAME -> NAME
%import common.NUMBER
%import common.ESCAPED_STRING
%import common.INT
%import common.WS
BOOL: "true" | "false"
COMMENT: /#.*/
%ignore COMMENT
%ignore WS
"""
def tok_to_bool(tok):
"Convert the value of `tok` from string to bool, while maintaining line number & column."
return Token.new_borrow_pos(tok.type, tok.value == "true", tok)
def tok_to_int(tok):
"Convert the value of `tok` from string to int, while maintaining line number & column."
return Token.new_borrow_pos(tok.type, int(tok), tok)
def tok_to_str(tok):
"Convert the value of `tok` from string to string, while maintaining line number & column."
return Token.new_borrow_pos(tok.type, tok.value.strip('"'), tok)
class ParserStage(object):
def __init__(self):
self.parser = Lark(calc_grammar,
parser='lalr',
lexer_callbacks={
'ESCAPED_STRING': tok_to_str,
'INT': tok_to_int,
'BOOL': tok_to_bool
})
def execute(self, terms_raw, ctx):
tree = self.parser.parse(terms_raw)
CollectDefinitions(ctx.primary, ctx.derived, ctx.averaged).visit(tree)
class DerivativeExpansionStage(object):
def execute(self, ctx):
# 1. 계산식 내에 존재하는 고차 차분 미분항(ddx, d2dy 등)을 찾아 중간 DerivedField로 등록
dset = set()
for k, v in ctx.derived.items():
dset.update(v.derivs)
for tup in dset:
a = DerivedField(tup[0], tup[1], ctx.derived)
ctx.derived[a.name] = a
# 2. 통계 물리량 평균화 대상 변수들을 AveragedField 구조체로 구성하고,
# 변동량 계산이 필요한 경우 FluctuationField로 등록
averaged_raw = ctx.averaged
ctx.averaged = {}
for w, tgts in averaged_raw.items():
for t in tgts:
a = AveragedField(w, t, ctx.derived)
ctx.averaged[a.name] = a
# 평균 편차가 동반된 항들에 대해 FluctuationField 생성
for ff in a.fset:
b = FluctuationField(w, ff, a.fset, ctx.derived)
ctx.derived[b.name] = b
# 3. 프로그램 내 모든 필드 간의 1차 의존 관계 그래프(Dependency Graph) 추출
ctx.dependency = {}
for k, v in ctx.derived.items():
ctx.dependency[k] = v.dep
for k, v in ctx.averaged.items():
ctx.dependency[k] = v.dep
class DependencyResolutionStage(object):
def execute(self, ctx):
# 이름 충돌 방지 검증
assert set(ctx.derived.keys()).isdisjoint(ctx.averaged.keys())
# Pass 1과 Pass 2 대상 평균화 변수 논리 분할
ctx.avg1 = set(filter(AveragedField.pass1, ctx.averaged.values()))
ctx.avg2 = set(filter(AveragedField.pass2, ctx.averaged.values()))
# Pass 1 위상 정렬: Pass 1 평균 변수들의 연산에 관여하는 모든 종속 관계를 수집하여 정렬
pass1calc = set(map(repr, ctx.avg1))
for x in ctx.avg1:
pass1calc.update(x.depClosure())
ctx.pass1 = self.sort_vars_new(ctx.dependency, pass1calc - ctx.primary)
# Pass 2 위상 정렬: Pass 2 평균 변수(변동 연산 연계)들의 연산에 관여하는 종속성을 정렬
pass2calc = set(map(repr, ctx.avg2))
for x in ctx.avg2:
pass2calc.update(x.depClosure())
ctx.pass2 = self.sort_vars_new(ctx.dependency, pass2calc - ctx.primary)
def calc_size (self, dependency, ordered, remaining):
count = 0
dep_union = set()
for v in remaining:
dep_union |= set(dependency[v])
for v in ordered:
if v in dep_union:
count += 1
return count
def sort_vars_new (self, dependency, group):
order = []
remain = list(group)
remain.sort()
while len(remain) > 0:
candidate = []
for v in remain:
if set(dependency[v]).isdisjoint(remain):
candidate.append(v)
impact = {}
size0 = self.calc_size(dependency, set(order), set(remain))
for v in candidate:
impact[v] = self.calc_size(dependency, set(order) | set([v]), set(remain) - set([v])) - size0
candidate.sort(key=impact.get)
order.append(candidate[0])
remain.remove(candidate[0])
return order
class SympyOptimizationStage(object):
def execute(self, ctx):
# 1. SymPy 수식 최적화 엔진 초기화
opt = SympyOptimizer.get_instance(ctx.derived)
opt.set_averaged(ctx.averaged)
# 2. SymPy가 대입식 치환(Substitution) 과정에서 제거한 불필요한 의존성 관계를
# 의존성 그래프에 즉각 반영하여 실제 계산을 위한 의존성 체인을 슬림하게 정리합니다.
updated_dependency = {}
for name, dep_set in ctx.dependency.items():
if name in ctx.derived and isinstance(ctx.derived[name], Field):
expr = opt.get_sympy_expr(name)
# 실제 정리된 SymPy 식에 잔존하는 자유 기호 명칭들만 추출
free_sym_names = {sym.name for sym in expr.free_symbols}
valid_deps = {dep for dep in free_sym_names if dep in ctx.derived or dep in ctx.primary}
updated_dependency[name] = valid_deps
else:
updated_dependency[name] = dep_set
ctx.dependency = updated_dependency
self.array_name = "xyzbuffer{}" # Pooling용 공용 3차원 버퍼 이름 포맷
# 3. Pass 1 및 Pass 2 연산 순서 배열들에 대해 각각 버퍼 공유 매핑(Pooling) 수행
narr1, alloc1 = (self.allocate_arr(ctx, ctx.pass1))
narr2, alloc2 = (self.allocate_arr(ctx, ctx.pass2))
# 전체 프로그램에서 필요한 동적 공유 3D 버퍼 배열의 최대 크기 설정
ctx.narr = max(narr1, narr2)
ctx.alloc1 = alloc1 # Pass 1 변수명 -> 공용 xyzbuffer 명칭 매핑 정보
ctx.alloc2 = alloc2 # Pass 2 변수명 -> 공용 xyzbuffer 명칭 매핑 정보
def liveness (self, ctx, l1, g):
import numpy as np
img = np.zeros((len(l1), len(l1)))
for i, v in enumerate(l1):
for j in range(i, len(l1)):
img[i,i:j+1] = img[i,i:j+1] + (1 if v in g[l1[j]] else 0)
return img > 0
def allocate_arr (self, ctx, l):
import numpy as np
dg = ctx.dependency
mask = self.liveness(ctx, l, dg)
try:
narr = mask.astype(int).sum(axis=0).max()
except ValueError:
narr = 0
array_pool = set([self.array_name.format(i) for i in range(narr)])
livesets = [set([])] + [set(np.asarray(l)[row]) for row in mask.T]
var2arr = { p : p for p in ctx.primary }
for i, (s0, s1) in enumerate(zip(livesets[:-1], livesets[1:])):
array_pool.update(map(var2arr.get, s0 - s1))
for new in s1 - s0:
var2arr[new] = array_pool.pop()
return narr, var2arr
class FortranProgramWriter(object):
def write(self, ctx):
from resources.m_template import mod_form
allvar = dict(ctx.derived)
allvar.update(ctx.averaged)
generator = FortranCodeGenerator(allvar)
# 평균 변수들을 순서대로 분배
set1 = sorted([a.name for a in filter(AveragedField.pass1, ctx.averaged.values())])
set2 = sorted([a.name for a in filter(AveragedField.pass2, ctx.averaged.values())])
# 외부 디스크 파일 익스포트 활성화 여부 확인
set_export_on = list(filter(lambda x: x.export_on(), ctx.derived.values()))
ffmt = 'logical, parameter :: pass2_required={}'
declf = ffmt.format('.true.' if len(set2) > 0 else '.false.')
hfmt = 'character (len = *), parameter :: output_header="{}"'
declh = hfmt.format(" ".join(["x"] + set1 + set2))
# 공용 Pooling 3D 버퍼 배열의 선언/할당 코드 생성
declarr, allocarr, freearr = self.work_array_codes(ctx.narr)
# 1차원 평균 물리량 배열들의 선언/할당 코드 생성
declavg = "\n".join(generator.generate_decl(ctx.averaged[v]) for v in sorted(ctx.averaged))
allocavg = "\n".join(generator.generate_alloc(ctx.averaged[v]) for v in sorted(ctx.averaged))
freeavg = "\n".join(generator.generate_free(ctx.averaged[v]) for v in sorted(ctx.averaged))
# 병렬 파일 쓰기(MPI Subarray)를 위한 MPI 리소스 선언/할당 코드 생성
decl_export = "\n".join(generator.generate_decl(v.exporter) for v in set_export_on)
alloc_export = "\n".join(generator.generate_alloc(v.exporter) for v in set_export_on)
free_export = "\n".join(generator.generate_free(v.exporter) for v in set_export_on)
# Pass 1과 Pass 2 루프 내부 본문 연산 코드들을 생성 (각자 버퍼 맵 alloc1, alloc2 적용)
sub_calc1 = "\n".join(generator.generate_code(allvar[v], ctx.alloc1) for v in ctx.pass1 if v in ctx.averaged or v in ctx.alloc1)
sub_calc2 = "\n".join(generator.generate_code(allvar[v], ctx.alloc2) for v in ctx.pass2 if v in ctx.averaged or v in ctx.alloc2)
# 평균 누적 연산 코드 생성
sub_avg1 = "\n".join(generator.generate_avg(allvar[v]) for v in set1)
sub_avg2 = "\n".join(generator.generate_avg(allvar[v]) for v in set2)
# 파일 쓰기 루틴 코드 생성
sub_write_avg = generator.generate_write_avg(set1+set2)
md = {}
md["module_name"] = "terms"
md["module_data"] = "\n".join((declf, declh, declavg, FieldExporter.mpi_io_decl, decl_export, declarr))
md["module_init"] = "\n".join((allocavg, alloc_export, allocarr))
md["module_finalize"] = "\n".join((freeavg, free_export, freearr))
md["module_pass1"] = sub_calc1
md["module_pass1_avg"] = sub_avg1
md["module_pass2"] = sub_calc2
md["module_pass2_avg"] = sub_avg2
md["module_write_result"] = sub_write_avg
print(Template(mod_form).render(**md))
def work_array_codes(self, narr):
array_name = "xyzbuffer{}"
array_names = [array_name.format(i) for i in range(narr)]
real_array_decl = "real(real64), allocatable, dimension(:,:,:) :: {0}"
decl = "\n".join([real_array_decl.format(v) for v in array_names])
alloc = "\n".join([make_allocate(v, "nxp,nyp,nzp") for v in array_names])
free = "\n".join(["deallocate({})".format(v) for v in array_names])
return decl, alloc, free
class ReportWriter(object):
def write(self, ctx):
import json
dg = {k:list(v) for k,v in ctx.dependency.items()}
with open("ir2.py", "w") as irf:
print("g = ", json.dumps(dg, indent=4), file=irf)
print("l1 = ", json.dumps(ctx.pass1, indent=4), file=irf)
print("l2 = ", json.dumps(ctx.pass2, indent=4), file=irf)
print("avg1 = ", json.dumps(list(map(repr, ctx.avg1)), indent=4), file=irf)
print("avg2 = ", json.dumps(list(map(repr, ctx.avg2)), indent=4), file=irf)
class LatexWriter(object):
def write(self, ctx):
latex_lines = ["{"]
for avg in ctx.averaged.values():
latex_lines.append(' "{}" : r"${}$",'.format(avg.name, avg.latex))
latex_lines.append("}")
print("\n".join(latex_lines))
Reference in a new issue View git blame Copy permalink