diff --git a/liskov/src/icon4py/liskov/logger.py b/common/src/icon4py/common/logger.py similarity index 100% rename from liskov/src/icon4py/liskov/logger.py rename to common/src/icon4py/common/logger.py diff --git a/liskov/README.md b/liskov/README.md index ffa4a5580f..ca8236fc8f 100644 --- a/liskov/README.md +++ b/liskov/README.md @@ -15,10 +15,18 @@ The icon4py-liskov package includes the `icon_liskov` CLI tool which takes a for To use the `icon_liskov` tool, run the following command: ```bash -icon_liskov [--profile] [--metadatagen] +icon_liskov [--profile] [--metadatagen] [--ppser] ``` -Where `input_filepath` is the path to the input file to be processed, and `output_filepath` is the path to the output file. The optional `--profile` flag adds nvtx profile statements to the stencils. +The following are descriptions of the arguments and options: + +- input_filepath: path to the input file to be processed. +- output_filepath: path to the output file. +- profile flag: adds nvtx profile statements to the stencils (optional). +- metadatagen flag: generates a metadata header at the top of the file which includes information on icon_liskov such as the version used. +- ppser flag: activates serialisation mode and will trigger the generation of ppser serialisation statements serialising all variables at the start and end of each stencil directive. + +**Note**: By default the data will be saved at the default folder location of the currently run experiment and will have a prefix of `liskov-serialisation`. ### Preprocessor directives diff --git a/liskov/src/icon4py/liskov/cli.py b/liskov/src/icon4py/liskov/cli.py index 8e0fe27436..b53eed9d8f 100644 --- a/liskov/src/icon4py/liskov/cli.py +++ b/liskov/src/icon4py/liskov/cli.py @@ -15,8 +15,8 @@ import click -from icon4py.liskov.logger import setup_logger -from icon4py.liskov.pipeline import ( +from icon4py.common.logger import setup_logger +from icon4py.liskov.pipeline.collection import ( load_gt4py_stencils, parse_fortran_file, run_code_generation, @@ -26,51 +26,82 @@ logger = setup_logger(__name__) -@click.command("icon_liskov") +@click.group(invoke_without_command=True) +@click.pass_context @click.argument( - "input_filepath", + "input_path", type=click.Path( exists=True, dir_okay=False, resolve_path=True, path_type=pathlib.Path ), ) @click.argument( - "output_filepath", + "output_path", type=click.Path(dir_okay=False, resolve_path=True, path_type=pathlib.Path), ) +def main(ctx, input_path, output_path): + """Command line interface for interacting with the ICON-Liskov DSL Preprocessor. + + Arguments: + INPUT_PATH: Path to input file containing Liskov directives. + OUTPUT_PATH: Path to new file to be generated. + """ + if ctx.invoked_subcommand is None: + click.echo( + "Need to choose one of the following commands:\nintegrate\nserialise" + ) + else: + ctx.ensure_object(dict) + ctx.obj["INPUT"] = input_path + ctx.obj["OUTPUT"] = output_path + + +@main.command() @click.option( - "--profile", "-p", is_flag=True, help="Add nvtx profile statements to stencils." + "--profile", + "-p", + is_flag=True, + help="Add nvtx profile statements to integration code.", ) @click.option( "--metadatagen", "-m", is_flag=True, - help="Add metadata header with information about program (requires git).", + help="Add metadata header with information about program.", ) -def main( - input_filepath: pathlib.Path, - output_filepath: pathlib.Path, - profile: bool, - metadatagen: bool, -) -> None: - """Command line interface for interacting with the ICON-Liskov DSL Preprocessor. - - Usage: - icon_liskov [-p] [-m] +@click.pass_context +def integrate(ctx, profile, metadatagen): + mode = "integration" + inp = ctx.obj["INPUT"] + out = ctx.obj["OUTPUT"] - Options: - -p --profile Add nvtx profile statements to stencils. - -m --metadatagen Add metadata header with information about program (requires git). - - Arguments: - input_filepath Path to the input file to process. - output_filepath Path to the output file to generate. - """ - parsed = parse_fortran_file(input_filepath, output_filepath) - parsed_checked = load_gt4py_stencils(parsed) + iface = parse_fortran_file(inp, out, mode) + iface_gt4py = load_gt4py_stencils(iface) run_code_generation( - parsed_checked, input_filepath, output_filepath, profile, metadatagen + inp, + out, + mode, + iface_gt4py, + profile=profile, + metadatagen=metadatagen, ) +@main.command() +@click.option( + "--multinode", + is_flag=True, + type=bool, + help="Specify whether it is a multinode run. Will generate mpi rank information.", + default=False, +) +@click.pass_context +def serialise(ctx, multinode): + mode = "serialisation" + inp = ctx.obj["INPUT"] + out = ctx.obj["OUTPUT"] + iface = parse_fortran_file(inp, out, mode) + run_code_generation(inp, out, mode, iface, multinode=multinode) + + if __name__ == "__main__": main() diff --git a/liskov/src/icon4py/liskov/py.typed.py b/liskov/src/icon4py/liskov/codegen/integration/__init__.py similarity index 100% rename from liskov/src/icon4py/liskov/py.typed.py rename to liskov/src/icon4py/liskov/codegen/integration/__init__.py diff --git a/liskov/src/icon4py/liskov/parsing/deserialise.py b/liskov/src/icon4py/liskov/codegen/integration/deserialise.py similarity index 81% rename from liskov/src/icon4py/liskov/parsing/deserialise.py rename to liskov/src/icon4py/liskov/codegen/integration/deserialise.py index fa88ae048c..46d09eaca8 100644 --- a/liskov/src/icon4py/liskov/parsing/deserialise.py +++ b/liskov/src/icon4py/liskov/codegen/integration/deserialise.py @@ -11,15 +11,14 @@ # # SPDX-License-Identifier: GPL-3.0-or-later -from dataclasses import dataclass -from typing import Any, Callable, Optional, Protocol, Type +from typing import Any, Optional, Protocol, Type +import icon4py.liskov.parsing.parse import icon4py.liskov.parsing.types as ts -from icon4py.liskov.codegen.interface import ( +from icon4py.common.logger import setup_logger +from icon4py.liskov.codegen.integration.interface import ( BoundsData, - CodeGenInput, DeclareData, - DeserialisedDirectives, EndCreateData, EndIfData, EndProfileData, @@ -27,13 +26,14 @@ FieldAssociationData, ImportsData, InsertData, + IntegrationCodeInterface, StartCreateData, StartProfileData, StartStencilData, UnusedDirective, ) -from icon4py.liskov.common import Step -from icon4py.liskov.logger import setup_logger +from icon4py.liskov.codegen.shared.deserialise import Deserialiser +from icon4py.liskov.codegen.shared.types import CodeGenInput from icon4py.liskov.parsing.exceptions import ( DirectiveSyntaxError, MissingBoundsError, @@ -89,13 +89,11 @@ def __call__( ... -@dataclass class DataFactoryBase: directive_cls: Type[ts.ParsedDirective] dtype: Type[CodeGenInput] -@dataclass class OptionalMultiUseDataFactory(DataFactoryBase): def __call__( self, parsed: ts.ParsedDict, **kwargs: Any @@ -106,48 +104,38 @@ def __call__( else: deserialised = [] for directive in extracted: - deserialised.append( - self.dtype( - startln=directive.startln, endln=directive.endln, **kwargs - ) - ) + deserialised.append(self.dtype(startln=directive.startln, **kwargs)) return deserialised -@dataclass class RequiredSingleUseDataFactory(DataFactoryBase): def __call__(self, parsed: ts.ParsedDict) -> CodeGenInput: extracted = extract_directive(parsed["directives"], self.directive_cls)[0] - return self.dtype(startln=extracted.startln, endln=extracted.endln) + return self.dtype(startln=extracted.startln) -@dataclass class EndCreateDataFactory(RequiredSingleUseDataFactory): - directive_cls: Type[ts.ParsedDirective] = ts.EndCreate + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.EndCreate dtype: Type[EndCreateData] = EndCreateData -@dataclass class ImportsDataFactory(RequiredSingleUseDataFactory): - directive_cls: Type[ts.ParsedDirective] = ts.Imports + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.Imports dtype: Type[ImportsData] = ImportsData -@dataclass class EndIfDataFactory(OptionalMultiUseDataFactory): - directive_cls: Type[ts.ParsedDirective] = ts.EndIf + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.EndIf dtype: Type[EndIfData] = EndIfData -@dataclass class EndProfileDataFactory(OptionalMultiUseDataFactory): - directive_cls: Type[ts.ParsedDirective] = ts.EndProfile + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.EndProfile dtype: Type[EndProfileData] = EndProfileData -@dataclass class StartCreateDataFactory(DataFactoryBase): - directive_cls: Type[ts.ParsedDirective] = ts.StartCreate + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.StartCreate dtype: Type[StartCreateData] = StartCreateData def __call__(self, parsed: ts.ParsedDict) -> StartCreateData: @@ -159,14 +147,11 @@ def __call__(self, parsed: ts.ParsedDict) -> StartCreateData: if named_args: extra_fields = named_args["extra_fields"].split(",") - return self.dtype( - startln=directive.startln, endln=directive.endln, extra_fields=extra_fields - ) + return self.dtype(startln=directive.startln, extra_fields=extra_fields) -@dataclass class DeclareDataFactory(DataFactoryBase): - directive_cls: Type[ts.ParsedDirective] = ts.Declare + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.Declare dtype: Type[DeclareData] = DeclareData @staticmethod @@ -185,7 +170,6 @@ def __call__(self, parsed: ts.ParsedDict) -> list[DeclareData]: deserialised.append( self.dtype( startln=directive.startln, - endln=directive.endln, declarations=named_args, ident_type=ident_type, suffix=suffix, @@ -194,9 +178,8 @@ def __call__(self, parsed: ts.ParsedDict) -> list[DeclareData]: return deserialised -@dataclass class StartProfileDataFactory(DataFactoryBase): - directive_cls: Type[ts.ParsedDirective] = ts.StartProfile + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.StartProfile dtype: Type[StartProfileData] = StartProfileData def __call__(self, parsed: ts.ParsedDict) -> list[StartProfileData]: @@ -206,18 +189,13 @@ def __call__(self, parsed: ts.ParsedDict) -> list[StartProfileData]: named_args = parsed["content"]["StartProfile"][i] stencil_name = _extract_stencil_name(named_args, directive) deserialised.append( - self.dtype( - name=stencil_name, - startln=directive.startln, - endln=directive.endln, - ) + self.dtype(name=stencil_name, startln=directive.startln) ) return deserialised -@dataclass class EndStencilDataFactory(DataFactoryBase): - directive_cls: Type[ts.ParsedDirective] = ts.EndStencil + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.EndStencil dtype: Type[EndStencilData] = EndStencilData def __call__(self, parsed: ts.ParsedDict) -> list[EndStencilData]: @@ -232,7 +210,6 @@ def __call__(self, parsed: ts.ParsedDict) -> list[EndStencilData]: self.dtype( name=stencil_name, startln=directive.startln, - endln=directive.endln, noendif=noendif, noprofile=noprofile, ) @@ -240,9 +217,8 @@ def __call__(self, parsed: ts.ParsedDict) -> list[EndStencilData]: return deserialised -@dataclass class StartStencilDataFactory(DataFactoryBase): - directive_cls: Type[ts.ParsedDirective] = ts.StartStencil + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.StartStencil dtype: Type[StartStencilData] = StartStencilData def __call__(self, parsed: ts.ParsedDict) -> list[StartStencilData]: @@ -282,7 +258,6 @@ def __call__(self, parsed: ts.ParsedDict) -> list[StartStencilData]: fields=fields_w_tolerance, bounds=bounds, startln=directive.startln, - endln=directive.endln, acc_present=acc_present, mergecopy=mergecopy, copies=copies, @@ -377,9 +352,8 @@ def _update_tolerances( return fields -@dataclass class InsertDataFactory(DataFactoryBase): - directive_cls: Type[ts.ParsedDirective] = ts.Insert + directive_cls: Type[ts.ParsedDirective] = icon4py.liskov.parsing.parse.Insert dtype: Type[InsertData] = InsertData def __call__(self, parsed: ts.ParsedDict) -> list[InsertData]: @@ -388,15 +362,13 @@ def __call__(self, parsed: ts.ParsedDict) -> list[InsertData]: for i, directive in enumerate(extracted): content = parsed["content"]["Insert"][i] deserialised.append( - self.dtype( - startln=directive.startln, endln=directive.endln, content=content # type: ignore - ) + self.dtype(startln=directive.startln, content=content) # type: ignore ) return deserialised -class DirectiveDeserialiser(Step): - _FACTORIES: dict[str, Callable] = { +class IntegrationCodeDeserialiser(Deserialiser): + _FACTORIES = { "StartCreate": StartCreateDataFactory(), "EndCreate": EndCreateDataFactory(), "Imports": ImportsDataFactory(), @@ -408,27 +380,4 @@ class DirectiveDeserialiser(Step): "EndProfile": EndProfileDataFactory(), "Insert": InsertDataFactory(), } - - def __call__(self, directives: ts.ParsedDict) -> DeserialisedDirectives: - """Deserialise the provided parsed directives to a DeserialisedDirectives object. - - Args: - directives: The parsed directives to deserialise. - - Returns: - A DeserialisedDirectives object containing the deserialised directives. - - Note: - The method uses the `_FACTORIES` class attribute to create the appropriate - factory object for each directive type, and uses these objects to deserialise - the parsed directives. The DeserialisedDirectives class is a dataclass - containing the deserialised versions of the different directives. - """ - logger.info("Deserialising directives ...") - deserialised = dict() - - for key, func in self._FACTORIES.items(): - ser = func(directives) - deserialised[key] = ser - - return DeserialisedDirectives(**deserialised) + _INTERFACE_TYPE = IntegrationCodeInterface diff --git a/liskov/src/icon4py/liskov/codegen/exceptions.py b/liskov/src/icon4py/liskov/codegen/integration/exceptions.py similarity index 100% rename from liskov/src/icon4py/liskov/codegen/exceptions.py rename to liskov/src/icon4py/liskov/codegen/integration/exceptions.py diff --git a/liskov/src/icon4py/liskov/codegen/generate.py b/liskov/src/icon4py/liskov/codegen/integration/generate.py similarity index 60% rename from liskov/src/icon4py/liskov/codegen/generate.py rename to liskov/src/icon4py/liskov/codegen/integration/generate.py index 26db2e626f..cc8af30138 100644 --- a/liskov/src/icon4py/liskov/codegen/generate.py +++ b/liskov/src/icon4py/liskov/codegen/integration/generate.py @@ -11,13 +11,15 @@ # # SPDX-License-Identifier: GPL-3.0-or-later -from dataclasses import dataclass -from typing import Any, Optional, Sequence, Type +from typing import Any -import gt4py.eve as eve -from gt4py.eve.codegen import TemplatedGenerator - -from icon4py.liskov.codegen.f90 import ( +from icon4py.common.logger import setup_logger +from icon4py.liskov.codegen.integration.interface import ( + IntegrationCodeInterface, + StartStencilData, + UnusedDirective, +) +from icon4py.liskov.codegen.integration.template import ( DeclareStatement, DeclareStatementGenerator, EndCreateStatement, @@ -40,42 +42,26 @@ StartProfileStatementGenerator, StartStencilStatement, StartStencilStatementGenerator, - generate_fortran_code, ) -from icon4py.liskov.codegen.interface import ( - CodeGenInput, - DeserialisedDirectives, - StartStencilData, - UnusedDirective, -) -from icon4py.liskov.common import Step +from icon4py.liskov.codegen.shared.generate import CodeGenerator +from icon4py.liskov.codegen.shared.types import GeneratedCode from icon4py.liskov.external.metadata import CodeMetadata -from icon4py.liskov.logger import setup_logger logger = setup_logger(__name__) -@dataclass -class GeneratedCode: - """A class for storing generated f90 code and its line number information.""" - - source: str - startln: int - endln: int - - -class IntegrationGenerator(Step): +class IntegrationCodeGenerator(CodeGenerator): def __init__( self, - directives: DeserialisedDirectives, - profile: bool, - metadata_gen: bool, + interface: IntegrationCodeInterface, + profile: bool = False, + metadatagen: bool = False, ): + super().__init__() self.profile = profile - self.directives = directives - self.generated: list[GeneratedCode] = [] - self.metadata_gen = metadata_gen + self.interface = interface + self.metadatagen = metadatagen def __call__(self, data: Any = None) -> list[GeneratedCode]: """Generate all f90 code for integration. @@ -94,48 +80,25 @@ def __call__(self, data: Any = None) -> list[GeneratedCode]: self._generate_insert() return self.generated - def _generate( - self, - parent_node: Type[eve.Node], - code_generator: Type[TemplatedGenerator], - startln: int, - endln: int, - **kwargs: CodeGenInput | Sequence[CodeGenInput] | Optional[bool] | Any, - ) -> None: - """Add a GeneratedCode object to the `generated` attribute with the given source code and line number information. - - Args: - parent_node: The parent node of the code to be generated. - code_generator: The code generator to use for generating the code. - startln: The start line number of the generated code. - endln: The end line number of the generated code. - **kwargs: Additional keyword arguments to be passed to the code generator. - """ - source = generate_fortran_code(parent_node, code_generator, **kwargs) - code = GeneratedCode(source=source, startln=startln, endln=endln) - self.generated.append(code) - def _generate_metadata(self) -> None: """Generate metadata about the current liskov execution.""" - if self.metadata_gen: + if self.metadatagen: logger.info("Generating icon-liskov metadata.") self._generate( MetadataStatement, MetadataStatementGenerator, startln=0, - endln=0, metadata=CodeMetadata(), ) def _generate_declare(self) -> None: """Generate f90 code for declaration statements.""" - for i, declare in enumerate(self.directives.Declare): + for i, declare in enumerate(self.interface.Declare): logger.info("Generating DECLARE statement.") self._generate( DeclareStatement, DeclareStatementGenerator, - self.directives.Declare[i].startln, - self.directives.Declare[i].endln, + self.interface.Declare[i].startln, declare_data=declare, ) @@ -147,19 +110,18 @@ def _generate_start_stencil(self) -> None: """ i = 0 - while i < len(self.directives.StartStencil): - stencil = self.directives.StartStencil[i] + while i < len(self.interface.StartStencil): + stencil = self.interface.StartStencil[i] logger.info(f"Generating START statement for {stencil.name}") try: - next_stencil = self.directives.StartStencil[i + 1] + next_stencil = self.interface.StartStencil[i + 1] except IndexError: pass if stencil.mergecopy and next_stencil.mergecopy: stencil = StartStencilData( startln=stencil.startln, - endln=next_stencil.endln, name=stencil.name + "_" + next_stencil.name, fields=stencil.fields + next_stencil.fields, bounds=stencil.bounds, @@ -173,7 +135,6 @@ def _generate_start_stencil(self) -> None: StartStencilStatement, StartStencilStatementGenerator, stencil.startln, - next_stencil.endln, stencil_data=stencil, profile=self.profile, ) @@ -181,8 +142,7 @@ def _generate_start_stencil(self) -> None: self._generate( StartStencilStatement, StartStencilStatementGenerator, - self.directives.StartStencil[i].startln, - self.directives.StartStencil[i].endln, + self.interface.StartStencil[i].startln, stencil_data=stencil, profile=self.profile, ) @@ -194,17 +154,16 @@ def _generate_end_stencil(self) -> None: Args: profile: A boolean indicating whether to include profiling calls in the generated code. """ - for i, stencil in enumerate(self.directives.StartStencil): + for i, stencil in enumerate(self.interface.StartStencil): logger.info(f"Generating END statement for {stencil.name}") self._generate( EndStencilStatement, EndStencilStatementGenerator, - self.directives.EndStencil[i].startln, - self.directives.EndStencil[i].endln, + self.interface.EndStencil[i].startln, stencil_data=stencil, profile=self.profile, - noendif=self.directives.EndStencil[i].noendif, - noprofile=self.directives.EndStencil[i].noprofile, + noendif=self.interface.EndStencil[i].noendif, + noprofile=self.interface.EndStencil[i].noprofile, ) def _generate_imports(self) -> None: @@ -213,9 +172,8 @@ def _generate_imports(self) -> None: self._generate( ImportsStatement, ImportsStatementGenerator, - self.directives.Imports.startln, - self.directives.Imports.endln, - stencils=self.directives.StartStencil, + self.interface.Imports.startln, + stencils=self.interface.StartStencil, ) def _generate_create(self) -> None: @@ -224,64 +182,52 @@ def _generate_create(self) -> None: self._generate( StartCreateStatement, StartCreateStatementGenerator, - self.directives.StartCreate.startln, - self.directives.StartCreate.endln, - stencils=self.directives.StartStencil, - extra_fields=self.directives.StartCreate.extra_fields, + self.interface.StartCreate.startln, + stencils=self.interface.StartStencil, + extra_fields=self.interface.StartCreate.extra_fields, ) self._generate( EndCreateStatement, EndCreateStatementGenerator, - self.directives.EndCreate.startln, - self.directives.EndCreate.endln, + self.interface.EndCreate.startln, ) def _generate_endif(self) -> None: """Generate f90 code for endif statements.""" - if self.directives.EndIf != UnusedDirective: - for endif in self.directives.EndIf: # type: ignore + if self.interface.EndIf != UnusedDirective: + for endif in self.interface.EndIf: # type: ignore logger.info("Generating ENDIF statement.") - self._generate( - EndIfStatement, - EndIfStatementGenerator, - endif.startln, - endif.endln, - ) + self._generate(EndIfStatement, EndIfStatementGenerator, endif.startln) def _generate_profile(self) -> None: """Generate additional nvtx profiling statements.""" if self.profile: - if self.directives.StartProfile != UnusedDirective: - for start in self.directives.StartProfile: # type: ignore + if self.interface.StartProfile != UnusedDirective: + for start in self.interface.StartProfile: # type: ignore logger.info("Generating nvtx start statement.") self._generate( StartProfileStatement, StartProfileStatementGenerator, start.startln, - start.endln, name=start.name, ) - if self.directives.EndProfile != UnusedDirective: - for end in self.directives.EndProfile: # type: ignore + if self.interface.EndProfile != UnusedDirective: + for end in self.interface.EndProfile: # type: ignore logger.info("Generating nvtx end statement.") self._generate( - EndProfileStatement, - EndProfileStatementGenerator, - end.startln, - end.endln, + EndProfileStatement, EndProfileStatementGenerator, end.startln ) def _generate_insert(self) -> None: """Generate free form statement from insert directive.""" - if self.directives.Insert != UnusedDirective: - for insert in self.directives.Insert: # type: ignore + if self.interface.Insert != UnusedDirective: + for insert in self.interface.Insert: # type: ignore logger.info("Generating free form statement.") self._generate( InsertStatement, InsertStatementGenerator, insert.startln, - insert.endln, content=insert.content, ) diff --git a/liskov/src/icon4py/liskov/codegen/interface.py b/liskov/src/icon4py/liskov/codegen/integration/interface.py similarity index 96% rename from liskov/src/icon4py/liskov/codegen/interface.py rename to liskov/src/icon4py/liskov/codegen/integration/interface.py index dfa879f0cd..e7fa5385e1 100644 --- a/liskov/src/icon4py/liskov/codegen/interface.py +++ b/liskov/src/icon4py/liskov/codegen/integration/interface.py @@ -14,17 +14,13 @@ from dataclasses import dataclass from typing import Optional, Sequence +from icon4py.liskov.codegen.shared.types import CodeGenInput + class UnusedDirective: ... -@dataclass -class CodeGenInput: - startln: int - endln: int - - @dataclass class BoundsData: hlower: str @@ -104,7 +100,7 @@ class InsertData(CodeGenInput): @dataclass -class DeserialisedDirectives: +class IntegrationCodeInterface: StartStencil: Sequence[StartStencilData] EndStencil: Sequence[EndStencilData] Declare: Sequence[DeclareData] diff --git a/liskov/src/icon4py/liskov/codegen/f90.py b/liskov/src/icon4py/liskov/codegen/integration/template.py similarity index 89% rename from liskov/src/icon4py/liskov/codegen/f90.py rename to liskov/src/icon4py/liskov/codegen/integration/template.py index 6ecc3e99ff..65c0785a0a 100644 --- a/liskov/src/icon4py/liskov/codegen/f90.py +++ b/liskov/src/icon4py/liskov/codegen/integration/template.py @@ -13,16 +13,14 @@ import re from dataclasses import asdict -from typing import Optional, Sequence, Type +from typing import Optional import gt4py.eve as eve from gt4py.eve.codegen import JinjaTemplate as as_jinja from gt4py.eve.codegen import TemplatedGenerator -from icon4py.icon4pygen.bindings.utils import format_fortran_code -from icon4py.liskov.codegen.exceptions import UndeclaredFieldError -from icon4py.liskov.codegen.interface import ( - CodeGenInput, +from icon4py.liskov.codegen.integration.exceptions import UndeclaredFieldError +from icon4py.liskov.codegen.integration.interface import ( DeclareData, StartStencilData, ) @@ -33,32 +31,6 @@ def enclose_in_parentheses(string: str) -> str: return f"({string})" -def generate_fortran_code( - parent_node: Type[eve.Node], - code_generator: Type[TemplatedGenerator], - **kwargs: CodeGenInput | Sequence[CodeGenInput] | Optional[bool], -) -> str: - """ - Generate Fortran code for the given parent node and code generator. - - Args: - parent_node: A subclass of eve.Node that represents the parent node. - code_generator: A subclass of TemplatedGenerator that will be used - to generate the code. - **kwargs: Arguments to be passed to the parent node constructor. - This can be a single CodeGenInput value, a sequence of CodeGenInput - values, or a boolean value, which is required by some parent nodes which - require a profile argument. - - Returns: - A string containing the formatted Fortran code. - """ - parent = parent_node(**kwargs) - source = code_generator.apply(parent) - formatted_source = format_fortran_code(source) - return formatted_source - - class BoundsFields(eve.Node): vlower: str vupper: str @@ -120,7 +92,7 @@ class MetadataStatementGenerator(TemplatedGenerator): ! GENERATED WITH ICON-LISKOV !+-+-+-+-+-+-+-+-+-+ +-+-+-+-+ +-+-+-+-+-+-+-+-+-+-+-+ ! Generated on: {{ _this_node.metadata.generated_on }} - ! Input filepath: {{ _this_node.metadata.cli_params['input_filepath'] }} + ! Input filepath: {{ _this_node.metadata.cli_params['input_path'] }} ! Profiling active: {{ _this_node.metadata.cli_params['profile'] }} ! Git version tag: {{ _this_node.metadata.version }} !+-+-+-+-+-+-+-+-+-+ +-+-+-+-+ +-+-+-+-+-+-+-+-+-+-+-+ diff --git a/liskov/src/icon4py/liskov/codegen/serialisation/__init__.py b/liskov/src/icon4py/liskov/codegen/serialisation/__init__.py new file mode 100644 index 0000000000..15dfdb0098 --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/serialisation/__init__.py @@ -0,0 +1,12 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later diff --git a/liskov/src/icon4py/liskov/codegen/serialisation/deserialise.py b/liskov/src/icon4py/liskov/codegen/serialisation/deserialise.py new file mode 100644 index 0000000000..f477747260 --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/serialisation/deserialise.py @@ -0,0 +1,225 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later +import uuid + +import icon4py.liskov.parsing.parse +import icon4py.liskov.parsing.types as ts +from icon4py.common.logger import setup_logger +from icon4py.liskov.codegen.integration.deserialise import ( + TOLERANCE_ARGS, + DeclareDataFactory, + _extract_stencil_name, + pop_item_from_dict, +) +from icon4py.liskov.codegen.serialisation.interface import ( + FieldSerialisationData, + ImportData, + InitData, + Metadata, + SavepointData, + SerialisationCodeInterface, +) +from icon4py.liskov.codegen.shared.deserialise import Deserialiser +from icon4py.liskov.parsing.utils import extract_directive + + +logger = setup_logger(__name__) + +KEYS_TO_REMOVE = [ + "accpresent", + "mergecopy", + "copies", + "horizontal_lower", + "horizontal_upper", + "vertical_lower", + "vertical_upper", + "name", +] + +SKIP_VARS = ["ikoffset", "pos_on_tplane_e_1", "pos_on_tplane_e_2"] + + +class InitDataFactory: + dtype = InitData + + def __call__(self, parsed: ts.ParsedDict) -> InitData: + return self.dtype(startln=0, directory=".", prefix="liskov-serialisation") + + +class SavepointDataFactory: + dtype = SavepointData + + def __call__(self, parsed: ts.ParsedDict) -> list[SavepointData]: + """Create a list of Start and End Savepoints for each Start and End Stencil directive.""" + start_stencil = extract_directive( + parsed["directives"], icon4py.liskov.parsing.parse.StartStencil + ) + end_stencil = extract_directive( + parsed["directives"], icon4py.liskov.parsing.parse.EndStencil + ) + gpu_fields = self.get_gpu_fields(parsed) + + repeated = self._find_repeated_stencils(parsed["content"]) + + deserialised = [] + + for i, (start, end) in enumerate(zip(start_stencil, end_stencil)): + named_args = parsed["content"]["StartStencil"][i] + stencil_name = _extract_stencil_name(named_args, start) + + if stencil_name in repeated: + stencil_name = f"{stencil_name}_{str(uuid.uuid4())}" + + field_names = self._remove_unnecessary_keys(named_args) + + metadata = [ + Metadata(key=k, value=v) + for k, v in self._get_timestep_variables(stencil_name).items() + ] + + fields = self._make_fields(field_names, gpu_fields) + + for intent, ln in [("start", start.startln), ("end", end.startln)]: + savepoint = self.dtype( + subroutine=f"{stencil_name}", + intent=intent, + startln=ln, + fields=fields, + metadata=metadata, + ) + deserialised.append(savepoint) + + return deserialised + + def get_gpu_fields(self, parsed: ts.ParsedDict) -> set[str]: + """Get declared fields which will be loaded on GPU and thus need to be serialised using accdata.""" + declare = DeclareDataFactory()(parsed) + fields = [] + for d in declare: + for f in d.declarations: + fields.append(f) + return set(fields) + + @staticmethod + def _remove_unnecessary_keys(named_args: dict) -> dict: + """Remove unnecessary keys from named_args, and only return field names.""" + copy = named_args.copy() + [pop_item_from_dict(copy, k, None) for k in KEYS_TO_REMOVE] + for tol in TOLERANCE_ARGS: + for k in copy.copy().keys(): + if k.endswith(tol): + pop_item_from_dict(copy, k, None) + return copy + + @staticmethod + def _make_fields(named_args: dict, gpu_fields: set) -> list[FieldSerialisationData]: + """Create a list of FieldSerialisationData objects based on named arguments.""" + fields = [ + FieldSerialisationData( + variable=variable, + association="z_hydro_corr(:,:,1)" # special case + if association == "z_hydro_corr(:,nlev,1)" + else association, + decomposed=False, + dimension=None, + typespec=None, + typename=None, + ptr_var=None, + device="gpu" if variable in gpu_fields else "cpu", + ) + for variable, association in named_args.items() + if variable not in SKIP_VARS + ] + return fields + + @staticmethod + def _get_timestep_variables(stencil_name: str) -> dict: + """Get the corresponding timestep metadata variables for the stencil.""" + timestep_variables = {} + + diffusion_stencil_names = [ + "apply", + "calculate", + "enhance", + "update", + "temporary", + "diffusion", + ] + if any(name in stencil_name for name in diffusion_stencil_names): + timestep_variables["jstep"] = "jstep_ptr" + timestep_variables["diffctr"] = "diffctr" + + if "mo_velocity_advection" in stencil_name: + timestep_variables["jstep"] = "jstep_ptr" + timestep_variables["nstep"] = "nstep_ptr" + timestep_variables["istep"] = "istep" + + if "mo_intp_rbf" in stencil_name: + timestep_variables["jstep"] = "jstep_ptr" + timestep_variables["mo_intp_rbf_ctr"] = "mo_intp_rbf_ctr" + + if "mo_math_divrot" in stencil_name: + timestep_variables["jstep"] = "jstep_ptr" + timestep_variables["mo_math_divrot_ctr"] = "mo_math_divrot_ctr" + + if "grad_green_gauss" in stencil_name: + timestep_variables["jstep"] = "jstep_ptr" + timestep_variables["grad_green_gauss_ctr"] = "grad_green_gauss_ctr" + + if "mo_icon_interpolation_scalar" in stencil_name: + timestep_variables["jstep"] = "jstep_ptr" + timestep_variables[ + "mo_icon_interpolation_ctr" + ] = "mo_icon_interpolation_ctr" + + if "mo_advection_traj" in stencil_name: + timestep_variables["jstep"] = "jstep_ptr" + timestep_variables["mo_advection_traj_ctr"] = "mo_advection_traj_ctr" + + if "mo_solve_nonhydro" in stencil_name: + timestep_variables["istep"] = "istep" + timestep_variables["mo_solve_nonhydro_ctr"] = "mo_solve_nonhydro_ctr" + + return timestep_variables + + def _find_repeated_stencils(self, content): + stencil_names = {} + repeated_names = [] + for stencil in content["StartStencil"]: + name = stencil["name"] + if name in stencil_names: + if stencil_names[name] not in repeated_names: + repeated_names.append(stencil_names[name]) + repeated_names.append(name) + else: + stencil_names[name] = name + return set(repeated_names) + + +class ImportDataFactory: + dtype = ImportData + + def __call__(self, parsed: ts.ParsedDict) -> ImportData: + imports = extract_directive( + parsed["directives"], icon4py.liskov.parsing.parse.Imports + )[0] + return self.dtype(startln=imports.startln) + + +class SerialisationCodeDeserialiser(Deserialiser): + _FACTORIES = { + "Init": InitDataFactory(), + "Savepoint": SavepointDataFactory(), + "Import": ImportDataFactory(), + } + _INTERFACE_TYPE = SerialisationCodeInterface diff --git a/liskov/src/icon4py/liskov/codegen/serialisation/generate.py b/liskov/src/icon4py/liskov/codegen/serialisation/generate.py new file mode 100644 index 0000000000..7fda920d2c --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/serialisation/generate.py @@ -0,0 +1,73 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later +from typing import Any + +from icon4py.common.logger import setup_logger +from icon4py.liskov.codegen.serialisation.interface import ( + SerialisationCodeInterface, +) +from icon4py.liskov.codegen.serialisation.template import ( + ImportStatement, + ImportStatementGenerator, + SavepointStatement, + SavepointStatementGenerator, +) +from icon4py.liskov.codegen.shared.generate import CodeGenerator +from icon4py.liskov.codegen.shared.types import GeneratedCode + + +logger = setup_logger(__name__) + + +class SerialisationCodeGenerator(CodeGenerator): + def __init__(self, interface: SerialisationCodeInterface, multinode: bool = False): + super().__init__() + self.interface = interface + self.multinode = multinode + + def __call__(self, data: Any = None) -> list[GeneratedCode]: + """Generate all f90 code for integration.""" + self._generate_import() + self._generate_savepoints() + return self.generated + + def _generate_import(self) -> None: + if self.multinode: + self._generate( + ImportStatement, + ImportStatementGenerator, + self.interface.Import.startln, + ) + + def _generate_savepoints(self) -> None: + init_complete = False + for i, savepoint in enumerate(self.interface.Savepoint): + logger.info("Generating pp_ser savepoint statement.") + if init_complete: + self._generate( + SavepointStatement, + SavepointStatementGenerator, + self.interface.Savepoint[i].startln, + savepoint=savepoint, + multinode=self.multinode, + ) + else: + self._generate( + SavepointStatement, + SavepointStatementGenerator, + self.interface.Savepoint[i].startln, + savepoint=savepoint, + init=self.interface.Init, + multinode=self.multinode, + ) + init_complete = True diff --git a/liskov/src/icon4py/liskov/codegen/serialisation/interface.py b/liskov/src/icon4py/liskov/codegen/serialisation/interface.py new file mode 100644 index 0000000000..301a4d2366 --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/serialisation/interface.py @@ -0,0 +1,60 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from dataclasses import dataclass +from typing import Optional + +from icon4py.liskov.codegen.shared.types import CodeGenInput + + +@dataclass +class InitData(CodeGenInput): + directory: str + prefix: str + + +@dataclass +class Metadata: + key: str + value: str + + +@dataclass +class FieldSerialisationData: + variable: str + association: str + decomposed: bool = False + device: Optional[str] = "cpu" + dimension: Optional[list[str]] = None + typespec: Optional[str] = None + typename: Optional[str] = None + ptr_var: Optional[str] = None + + +@dataclass +class SavepointData(CodeGenInput): + subroutine: str # todo: change to name + intent: str + fields: list[FieldSerialisationData] + metadata: Optional[list[Metadata]] + + +class ImportData(CodeGenInput): + ... + + +@dataclass +class SerialisationCodeInterface: + Import: ImportData + Init: InitData + Savepoint: list[SavepointData] diff --git a/liskov/src/icon4py/liskov/codegen/serialisation/template.py b/liskov/src/icon4py/liskov/codegen/serialisation/template.py new file mode 100644 index 0000000000..120e73bed2 --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/serialisation/template.py @@ -0,0 +1,142 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later +from dataclasses import asdict +from typing import Optional + +import gt4py.eve as eve +from gt4py.eve.codegen import JinjaTemplate as as_jinja +from gt4py.eve.codegen import TemplatedGenerator + +from icon4py.liskov.codegen.serialisation.interface import InitData, SavepointData + + +class Field(eve.Node): + variable: str + association: str + decomposed: bool + dimension: Optional[list[str]] + typespec: Optional[str] + typename: Optional[str] + ptr_var: Optional[str] + device: str + + +class StandardFields(eve.Node): + fields: list[Field] + + +class DecomposedFields(StandardFields): + ... + + +class DecomposedFieldDeclarations(DecomposedFields): + ... + + +class SavepointStatement(eve.Node): + savepoint: SavepointData + init: Optional[InitData] = eve.datamodels.field(default=None) + multinode: bool + standard_fields: StandardFields = eve.datamodels.field(init=False) + decomposed_fields: DecomposedFields = eve.datamodels.field(init=False) + decomposed_field_declarations: DecomposedFieldDeclarations = eve.datamodels.field( + init=False + ) + + def __post_init__(self): + self.standard_fields = StandardFields( + fields=[ + Field(**asdict(f)) for f in self.savepoint.fields if not f.decomposed + ] + ) + self.decomposed_fields = DecomposedFields( + fields=[Field(**asdict(f)) for f in self.savepoint.fields if f.decomposed] + ) + self.decomposed_field_declarations = DecomposedFieldDeclarations( + fields=[Field(**asdict(f)) for f in self.savepoint.fields if f.decomposed] + ) + + +class SavepointStatementGenerator(TemplatedGenerator): + SavepointStatement = as_jinja( + """ + {{ decomposed_field_declarations }} + + {% if _this_node.init %} + !$ser init directory="{{_this_node.init.directory}}" prefix="{{ _this_node.init.prefix }}" {% if _this_node.multinode %}mpi_rank=get_my_mpi_work_id(){% endif %} + {% endif %} + + !$ser savepoint {{ _this_node.savepoint.subroutine }}_{{ _this_node.savepoint.intent }} {% if _this_node.savepoint.metadata %} {%- for m in _this_node.savepoint.metadata -%} {{ m.key }}={{ m.value }} {% endfor -%} {% endif %} + + {{ decomposed_fields }} + + {{ standard_fields }} + """ + ) + + StandardFields = as_jinja( + """ + {% for f in _this_node.fields %} + PRINT *, 'Serializing {{ f.variable }}={{ f.association }}' + {% if f.dimension %} + IF (SIZE({{ f.variable }}) > 0) THEN + !$ser {% if f.device == 'gpu' %}accdata {% else %}data {% endif %}{{ f.variable }}={{ f.association }} + ELSE + PRINT *, 'Warning: Array {{ f.variable }} has size 0. Not serializing array.' + ENDIF + {% else %} + !$ser {% if f.device == 'gpu' %}accdata {% else %}data {% endif %}{{ f.variable }}={{ f.association }} + {% endif %} + {% endfor %} + """ + ) + + DecomposedFieldDeclarations = as_jinja( + """ + {% for f in _this_node.fields %} + !$ser verbatim {{ f.typespec }}, dimension({{ ",".join(f.dimension) }}), allocatable :: {{ f.variable }}_{{ f.ptr_var}}({{ ",".join(f.dimension) }}) + {% endfor %} + """ + ) + + DecomposedFields = as_jinja( + """ + {% for f in _this_node.fields %} + !$ser verbatim allocate({{ f.variable }}_{{ f.ptr_var}}({{ f.alloc_dims }})) + !$ser verbatim {{ f.variable }}_{{ f.ptr_var}} = {{ f.variable }}%{{ f.ptr_var}} + !$ser {% if f.device == 'gpu' %}accdata {% else %}data {% endif %}{{ f.variable }}_{{ f.ptr_var}}={{ f.variable }}_{{ f.ptr_var}} + !$ser verbatim deallocate({{ f.variable }}_{{ f.ptr_var}}) + {% endfor %} + """ + ) + + def visit_DecomposedFields(self, node: DecomposedFields): + def generate_size_strings(colon_list, var_name): + size_strings = [] + for i in range(len(colon_list)): + size_strings.append(f"size({var_name}, {i + 1})") + return size_strings + + for f in node.fields: + f.variable = f.variable.replace(f"_{f.ptr_var}", "") + f.alloc_dims = ",".join(generate_size_strings(f.dimension, f.variable)) + + return self.generic_visit(node) + + +class ImportStatement(eve.Node): + ... + + +class ImportStatementGenerator(TemplatedGenerator): + ImportStatement = as_jinja(" USE mo_mpi, ONLY: get_my_mpi_work_id") diff --git a/liskov/src/icon4py/liskov/codegen/shared/__init__.py b/liskov/src/icon4py/liskov/codegen/shared/__init__.py new file mode 100644 index 0000000000..15dfdb0098 --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/shared/__init__.py @@ -0,0 +1,12 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later diff --git a/liskov/src/icon4py/liskov/codegen/shared/deserialise.py b/liskov/src/icon4py/liskov/codegen/shared/deserialise.py new file mode 100644 index 0000000000..6cb02811ab --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/shared/deserialise.py @@ -0,0 +1,52 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from typing import Callable + +import icon4py.liskov.parsing.types as ts +from icon4py.common.logger import setup_logger +from icon4py.liskov.codegen.integration.interface import IntegrationCodeInterface +from icon4py.liskov.codegen.serialisation.interface import ( + SerialisationCodeInterface, +) +from icon4py.liskov.pipeline.definition import Step + + +logger = setup_logger(__name__) + + +class Deserialiser(Step): + _FACTORIES: dict[str, Callable] = {} + _INTERFACE_TYPE: SerialisationCodeInterface | IntegrationCodeInterface + + def __call__(self, directives: ts.ParsedDict): + """Deserialises parsed directives into an Interface object. + + Args: + directives: A dictionary of parsed directives. + + Returns: + An Interface object containing the deserialised directives. + + This method is responsible for deserialising parsed directives into an Interface object of the given _INTERFACE_TYPE. + It uses the `_FACTORIES` dictionary of factory functions to create the appropriate factory object for each directive type. + The resulting deserialised objects are then used to create an Interface object which can be used for code generation. + """ + logger.info("Deserialising directives ...") + deserialised = dict() + + for key, factory in self._FACTORIES.items(): + obj = factory(directives) + deserialised[key] = obj + + return self._INTERFACE_TYPE(**deserialised) diff --git a/liskov/src/icon4py/liskov/codegen/shared/generate.py b/liskov/src/icon4py/liskov/codegen/shared/generate.py new file mode 100644 index 0000000000..47931e3af4 --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/shared/generate.py @@ -0,0 +1,76 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import abc +from typing import Any, Optional, Sequence, Type + +import gt4py.eve as eve +from gt4py.eve.codegen import TemplatedGenerator + +from icon4py.icon4pygen.bindings.utils import format_fortran_code +from icon4py.liskov.codegen.shared.types import CodeGenInput, GeneratedCode +from icon4py.liskov.pipeline.definition import Step + + +class CodeGenerator(Step): + def __init__(self) -> None: + self.generated: list[GeneratedCode] = [] + + @abc.abstractmethod + def __call__(self, data: Any) -> list[GeneratedCode]: + ... + + @staticmethod + def _generate_fortran_code( + parent_node: Type[eve.Node], + code_generator: Type[TemplatedGenerator], + **kwargs: CodeGenInput | Sequence[CodeGenInput] | Optional[bool], + ) -> str: + """ + Generate Fortran code for the given parent node and code generator. + + Args: + parent_node: A subclass of eve.Node that represents the parent node. + code_generator: A subclass of TemplatedGenerator that will be used + to generate the code. + **kwargs: Arguments to be passed to the parent node constructor. + This can be a single CodeGenInput value, a sequence of CodeGenInput + values, or a boolean value, which is required by some parent nodes which + require a profile argument. + + Returns: + A string containing the formatted Fortran code. + """ + parent = parent_node(**kwargs) + source = code_generator.apply(parent) + formatted_source = format_fortran_code(source) + return formatted_source + + def _generate( + self, + parent_node: Type[eve.Node], + code_generator: Type[TemplatedGenerator], + startln: int, + **kwargs: CodeGenInput | Sequence[CodeGenInput] | Optional[bool] | Any, + ) -> None: + """Add a GeneratedCode object to the `generated` attribute with the given source code and line number information. + + Args: + parent_node: The parent node of the code to be generated. + code_generator: The code generator to use for generating the code. + startln: The start line number of the generated code. + **kwargs: Additional keyword arguments to be passed to the code generator. + """ + source = self._generate_fortran_code(parent_node, code_generator, **kwargs) + code = GeneratedCode(startln=startln, source=source) + self.generated.append(code) diff --git a/liskov/src/icon4py/liskov/codegen/shared/types.py b/liskov/src/icon4py/liskov/codegen/shared/types.py new file mode 100644 index 0000000000..7420b6f532 --- /dev/null +++ b/liskov/src/icon4py/liskov/codegen/shared/types.py @@ -0,0 +1,24 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from dataclasses import dataclass + + +@dataclass +class CodeGenInput: + startln: int + + +@dataclass +class GeneratedCode(CodeGenInput): + source: str diff --git a/liskov/src/icon4py/liskov/codegen/write.py b/liskov/src/icon4py/liskov/codegen/shared/write.py similarity index 88% rename from liskov/src/icon4py/liskov/codegen/write.py rename to liskov/src/icon4py/liskov/codegen/shared/write.py index 6ce3567510..588d90e6d3 100644 --- a/liskov/src/icon4py/liskov/codegen/write.py +++ b/liskov/src/icon4py/liskov/codegen/shared/write.py @@ -13,18 +13,18 @@ from pathlib import Path from typing import List -from icon4py.liskov.codegen.generate import GeneratedCode -from icon4py.liskov.common import Step -from icon4py.liskov.logger import setup_logger +from icon4py.common.logger import setup_logger +from icon4py.liskov.codegen.shared.types import GeneratedCode from icon4py.liskov.parsing.types import DIRECTIVE_IDENT +from icon4py.liskov.pipeline.definition import Step logger = setup_logger(__name__) -class IntegrationWriter(Step): +class CodegenWriter(Step): def __init__(self, input_filepath: Path, output_filepath: Path) -> None: - """Initialize an IntegrationWriter instance with a list of generated code. + """Initialize an CodegenWriter instance. Args: input_filepath: Path to file containing directives. @@ -34,10 +34,12 @@ def __init__(self, input_filepath: Path, output_filepath: Path) -> None: self.output_filepath = output_filepath def __call__(self, generated: List[GeneratedCode]) -> None: - """Write a file containing generated code, with the DSL directives removed in the same directory as filepath using a new suffix. + """Write a new file containing the generated code. + + Any !$DSL directives are removed from the file. Args: - generated: A list of GeneratedCode instances representing the generated code that will be written to a file. + generated: A list of GeneratedCode instances representing the generated code that will be written to the file. """ current_file = self._read_file() with_generated_code = self._insert_generated_code(current_file, generated) @@ -45,7 +47,7 @@ def __call__(self, generated: List[GeneratedCode]) -> None: self._write_file(without_directives) def _read_file(self) -> List[str]: - """Read the lines of a file into a list. + """Read the lines of the input file into a list. Returns: A list of strings representing the lines of the file. diff --git a/liskov/src/icon4py/liskov/external/gt4py.py b/liskov/src/icon4py/liskov/external/gt4py.py index c0096a92af..d02397b1db 100644 --- a/liskov/src/icon4py/liskov/external/gt4py.py +++ b/liskov/src/icon4py/liskov/external/gt4py.py @@ -17,14 +17,14 @@ from gt4py.next.ffront.decorator import Program +from icon4py.common.logger import setup_logger from icon4py.icon4pygen.metadata import get_stencil_info -from icon4py.liskov.codegen.interface import DeserialisedDirectives -from icon4py.liskov.common import Step +from icon4py.liskov.codegen.integration.interface import IntegrationCodeInterface from icon4py.liskov.external.exceptions import ( IncompatibleFieldError, UnknownStencilError, ) -from icon4py.liskov.logger import setup_logger +from icon4py.liskov.pipeline.definition import Step logger = setup_logger(__name__) @@ -33,10 +33,10 @@ class UpdateFieldsWithGt4PyStencils(Step): _STENCIL_PACKAGES = ["atm_dyn_iconam", "advection"] - def __init__(self, parsed: DeserialisedDirectives): + def __init__(self, parsed: IntegrationCodeInterface): self.parsed = parsed - def __call__(self, data: Any = None) -> DeserialisedDirectives: + def __call__(self, data: Any = None) -> IntegrationCodeInterface: logger.info("Updating parsed fields with data from icon4py stencils...") for s in self.parsed.StartStencil: diff --git a/liskov/src/icon4py/liskov/external/metadata.py b/liskov/src/icon4py/liskov/external/metadata.py index e1b5ff2c47..756d92dcea 100644 --- a/liskov/src/icon4py/liskov/external/metadata.py +++ b/liskov/src/icon4py/liskov/external/metadata.py @@ -30,7 +30,9 @@ def generated_on(self) -> str: def cli_params(self) -> dict[str, Any]: try: ctx = click.get_current_context() - return ctx.params + params = ctx.params.copy() + params.update(ctx.parent.params) + return params except Exception as e: raise MissingClickContextError( f"Cannot fetch click context in this thread as no click command has been executed.\n {e}" diff --git a/liskov/src/icon4py/liskov/parsing/parse.py b/liskov/src/icon4py/liskov/parsing/parse.py index 93b42d0cb0..3543fb7c37 100644 --- a/liskov/src/icon4py/liskov/parsing/parse.py +++ b/liskov/src/icon4py/liskov/parsing/parse.py @@ -13,14 +13,16 @@ import collections import shutil import sys +from dataclasses import dataclass, field from pathlib import Path -from typing import Sequence +from typing import Optional, Sequence, Type import icon4py.liskov.parsing.types as ts -from icon4py.liskov.common import Step -from icon4py.liskov.logger import setup_logger +from icon4py.common.logger import setup_logger from icon4py.liskov.parsing.exceptions import UnsupportedDirectiveError +from icon4py.liskov.parsing.types import ParsedDirective, RawDirective from icon4py.liskov.parsing.validation import VALIDATORS +from icon4py.liskov.pipeline.definition import Step REPLACE_CHARS = [ts.DIRECTIVE_IDENT, "&", "\n"] @@ -68,7 +70,7 @@ def _determine_type( typed = [] for raw in raw_directives: found = False - for directive in ts.SUPPORTED_DIRECTIVES: + for directive in SUPPORTED_DIRECTIVES: if directive.pattern in raw.string: typed.append(directive(raw.string, raw.startln, raw.endln)) # type: ignore found = True @@ -108,3 +110,128 @@ def _parse(directives: Sequence[ts.ParsedDirective]) -> ts.ParsedContent: content = d.get_content() parsed_content[d.type_name].append(content) return parsed_content + + +class TypedDirective(RawDirective): + pattern: str + + @property + def type_name(self) -> str: + return self.__class__.__name__ + + def __eq__(self, other: object) -> bool: + if not isinstance(other, TypedDirective): + raise NotImplementedError + return self.string == other.string + + +@dataclass(eq=False) +class WithArguments(TypedDirective): + regex: str = field(default=r"(.+?)=(.+?)", init=False) + + def get_content(self) -> dict[str, str]: + args = self.string.replace(f"{self.pattern}", "") + delimited = args[1:-1].split(";") + content = { + strip_whitespace(a.split("=")[0].strip()): strip_whitespace(a.split("=")[1]) + for a in delimited + } + return content + + +@dataclass(eq=False) +class WithOptionalArguments(TypedDirective): + regex: str = field(default=r"(?:.+?=.+?|)", init=False) + + def get_content(self) -> Optional[dict[str, str]]: + args = self.string.replace(f"{self.pattern}", "")[1:-1] + if len(args) > 0: + content = dict([args.split("=")]) + return content + return None + + +@dataclass(eq=False) +class WithoutArguments(TypedDirective): + # matches an empty string at the beginning of a line + regex: str = field(default=r"^(?![\s\S])", init=False) + + def get_content(self) -> dict: + return {} + + +@dataclass(eq=False) +class FreeForm(TypedDirective): + # matches any string inside brackets + regex: str = field(default=r"(.+?)", init=False) + + def get_content(self) -> str: + args = self.string.replace(f"{self.pattern}", "") + return args[1:-1] + + +def strip_whitespace(string: str) -> str: + """ + Remove all whitespace characters from the given string. + + Args: + string: The string to remove whitespace from. + + Returns: + The input string with all whitespace removed. + """ + return "".join(string.split()) + + +class StartStencil(WithArguments): + pattern = "START STENCIL" + + +class EndStencil(WithArguments): + pattern = "END STENCIL" + + +class Declare(WithArguments): + pattern = "DECLARE" + + +class Imports(WithoutArguments): + pattern = "IMPORTS" + + +class StartCreate(WithOptionalArguments): + pattern = "START CREATE" + + +class EndCreate(WithoutArguments): + pattern = "END CREATE" + + +class EndIf(WithoutArguments): + pattern = "ENDIF" + + +class StartProfile(WithArguments): + pattern = "START PROFILE" + + +class EndProfile(WithoutArguments): + pattern = "END PROFILE" + + +class Insert(FreeForm): + pattern = "INSERT" + + +SUPPORTED_DIRECTIVES: Sequence[Type[ParsedDirective]] = [ + StartStencil, + EndStencil, + Imports, + Declare, + StartCreate, + EndCreate, + EndIf, + StartProfile, + EndProfile, + Insert, +] diff --git a/liskov/src/icon4py/liskov/parsing/scan.py b/liskov/src/icon4py/liskov/parsing/scan.py index 84c44955eb..f13c1611aa 100644 --- a/liskov/src/icon4py/liskov/parsing/scan.py +++ b/liskov/src/icon4py/liskov/parsing/scan.py @@ -15,9 +15,9 @@ from typing import Any import icon4py.liskov.parsing.types as ts -from icon4py.liskov.common import Step -from icon4py.liskov.logger import setup_logger +from icon4py.common.logger import setup_logger from icon4py.liskov.parsing.exceptions import DirectiveSyntaxError +from icon4py.liskov.pipeline.definition import Step logger = setup_logger(__name__) diff --git a/liskov/src/icon4py/liskov/parsing/types.py b/liskov/src/icon4py/liskov/parsing/types.py index 9fb5ee29da..f24028ec01 100644 --- a/liskov/src/icon4py/liskov/parsing/types.py +++ b/liskov/src/icon4py/liskov/parsing/types.py @@ -10,13 +10,11 @@ # distribution for a copy of the license or check . # # SPDX-License-Identifier: GPL-3.0-or-later -from dataclasses import dataclass, field +from dataclasses import dataclass from typing import ( Any, - Optional, Protocol, Sequence, - Type, TypeAlias, TypedDict, runtime_checkable, @@ -55,113 +53,3 @@ class RawDirective: string: str startln: int endln: int - - -class TypedDirective(RawDirective): - pattern: str - - @property - def type_name(self) -> str: - return self.__class__.__name__ - - def __eq__(self, other: object) -> bool: - if not isinstance(other, TypedDirective): - raise NotImplementedError - return self.string == other.string - - -@dataclass(eq=False) -class WithArguments(TypedDirective): - regex: str = field(default=r"(.+?)=(.+?)", init=False) - - def get_content(self) -> dict[str, str]: - args = self.string.replace(f"{self.pattern}", "") - delimited = args[1:-1].split(";") - content = {a.split("=")[0].strip(): a.split("=")[1] for a in delimited} - return content - - -@dataclass(eq=False) -class WithOptionalArguments(TypedDirective): - regex: str = field(default=r"(?:.+?=.+?|)", init=False) - - def get_content(self) -> Optional[dict[str, str]]: - args = self.string.replace(f"{self.pattern}", "")[1:-1] - if len(args) > 0: - content = dict([args.split("=")]) - return content - return None - - -@dataclass(eq=False) -class WithoutArguments(TypedDirective): - # matches an empty string at the beginning of a line - regex: str = field(default=r"^(?![\s\S])", init=False) - - def get_content(self) -> dict: - return {} - - -@dataclass(eq=False) -class FreeForm(TypedDirective): - # matches any string inside brackets - regex: str = field(default=r"(.+?)", init=False) - - def get_content(self) -> str: - args = self.string.replace(f"{self.pattern}", "") - return args[1:-1] - - -class StartStencil(WithArguments): - pattern = "START STENCIL" - - -class EndStencil(WithArguments): - pattern = "END STENCIL" - - -class Declare(WithArguments): - pattern = "DECLARE" - - -class Imports(WithoutArguments): - pattern = "IMPORTS" - - -class StartCreate(WithOptionalArguments): - pattern = "START CREATE" - - -class EndCreate(WithoutArguments): - pattern = "END CREATE" - - -class EndIf(WithoutArguments): - pattern = "ENDIF" - - -class StartProfile(WithArguments): - pattern = "START PROFILE" - - -class EndProfile(WithoutArguments): - pattern = "END PROFILE" - - -class Insert(FreeForm): - pattern = "INSERT" - - -# When adding a new directive this list must be updated. -SUPPORTED_DIRECTIVES: Sequence[Type[ParsedDirective]] = [ - StartStencil, - EndStencil, - Imports, - Declare, - StartCreate, - EndCreate, - EndIf, - StartProfile, - EndProfile, - Insert, -] diff --git a/liskov/src/icon4py/liskov/parsing/utils.py b/liskov/src/icon4py/liskov/parsing/utils.py index 8a6f1d743b..35d4a1a4ed 100644 --- a/liskov/src/icon4py/liskov/parsing/utils.py +++ b/liskov/src/icon4py/liskov/parsing/utils.py @@ -12,7 +12,7 @@ # SPDX-License-Identifier: GPL-3.0-or-later from typing import Sequence, Type -import icon4py.liskov.parsing.types as ts +from icon4py.liskov.parsing import types as ts def flatten_list_of_dicts(list_of_dicts: list[dict]) -> dict: @@ -36,11 +36,6 @@ def string_to_bool(string: str) -> bool: raise ValueError(f"Cannot convert '{string}' to a boolean.") -def print_parsed_directive(directive: ts.ParsedDirective) -> str: - """Print a parsed directive, including its contents, and start and end line numbers.""" - return f"Directive: {directive.string}, start line: {directive.startln}, end line: {directive.endln}\n" - - def extract_directive( directives: Sequence[ts.ParsedDirective], required_type: Type[ts.ParsedDirective], @@ -50,6 +45,11 @@ def extract_directive( return directives +def print_parsed_directive(directive: ts.ParsedDirective) -> str: + """Print a parsed directive, including its contents, and start and end line numbers.""" + return f"Directive: {directive.string}, start line: {directive.startln}, end line: {directive.endln}\n" + + def remove_directive_types( directives: Sequence[ts.ParsedDirective], exclude_types: Sequence[Type[ts.ParsedDirective]], diff --git a/liskov/src/icon4py/liskov/parsing/validation.py b/liskov/src/icon4py/liskov/parsing/validation.py index 616803e709..70e4e4cde8 100644 --- a/liskov/src/icon4py/liskov/parsing/validation.py +++ b/liskov/src/icon4py/liskov/parsing/validation.py @@ -16,8 +16,9 @@ from pathlib import Path from typing import Match, Optional, Protocol +import icon4py.liskov.parsing.parse import icon4py.liskov.parsing.types as ts -from icon4py.liskov.logger import setup_logger +from icon4py.common.logger import setup_logger from icon4py.liskov.parsing.exceptions import ( DirectiveSyntaxError, RepeatedDirectiveError, @@ -117,12 +118,12 @@ def _validate_directive_uniqueness( repeated = remove_directive_types( [d for d in directives if directives.count(d) > 1], [ - ts.StartStencil, - ts.EndStencil, - ts.EndIf, - ts.EndProfile, - ts.StartProfile, - ts.Insert, + icon4py.liskov.parsing.parse.StartStencil, + icon4py.liskov.parsing.parse.EndStencil, + icon4py.liskov.parsing.parse.EndIf, + icon4py.liskov.parsing.parse.EndProfile, + icon4py.liskov.parsing.parse.StartProfile, + icon4py.liskov.parsing.parse.Insert, ], ) if repeated: @@ -136,12 +137,12 @@ def _validate_required_directives( ) -> None: """Check that all required directives are used at least once.""" expected = [ - ts.Declare, - ts.Imports, - ts.StartCreate, - ts.EndCreate, - ts.StartStencil, - ts.EndStencil, + icon4py.liskov.parsing.parse.Declare, + icon4py.liskov.parsing.parse.Imports, + icon4py.liskov.parsing.parse.StartCreate, + icon4py.liskov.parsing.parse.EndCreate, + icon4py.liskov.parsing.parse.StartStencil, + icon4py.liskov.parsing.parse.EndStencil, ] for expected_type in expected: if not any([isinstance(d, expected_type) for d in directives]): @@ -171,13 +172,23 @@ def _validate_stencil_directives( directives (list[ts.ParsedDirective]): List of stencil directives to validate. """ stencil_directives = [ - d for d in directives if isinstance(d, (ts.StartStencil, ts.EndStencil)) + d + for d in directives + if isinstance( + d, + ( + icon4py.liskov.parsing.parse.StartStencil, + icon4py.liskov.parsing.parse.EndStencil, + ), + ) ] stencil_counts: dict = {} for directive in stencil_directives: stencil_name = self.extract_arg_from_directive(directive.string, "name") stencil_counts[stencil_name] = stencil_counts.get(stencil_name, 0) + ( - 1 if isinstance(directive, ts.StartStencil) else -1 + 1 + if isinstance(directive, icon4py.liskov.parsing.parse.StartStencil) + else -1 ) unbalanced_stencils = [ diff --git a/liskov/src/icon4py/liskov/pipeline/__init__.py b/liskov/src/icon4py/liskov/pipeline/__init__.py new file mode 100644 index 0000000000..15dfdb0098 --- /dev/null +++ b/liskov/src/icon4py/liskov/pipeline/__init__.py @@ -0,0 +1,12 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later diff --git a/liskov/src/icon4py/liskov/pipeline.py b/liskov/src/icon4py/liskov/pipeline/collection.py similarity index 51% rename from liskov/src/icon4py/liskov/pipeline.py rename to liskov/src/icon4py/liskov/pipeline/collection.py index d452ce6f73..5bc730c36c 100644 --- a/liskov/src/icon4py/liskov/pipeline.py +++ b/liskov/src/icon4py/liskov/pipeline/collection.py @@ -10,49 +10,75 @@ # distribution for a copy of the license or check . # # SPDX-License-Identifier: GPL-3.0-or-later - from pathlib import Path -from icon4py.liskov.codegen.generate import IntegrationGenerator -from icon4py.liskov.codegen.interface import DeserialisedDirectives -from icon4py.liskov.codegen.write import IntegrationWriter -from icon4py.liskov.common import Step, linear_pipeline +from icon4py.liskov.codegen.integration.deserialise import ( + IntegrationCodeDeserialiser, +) +from icon4py.liskov.codegen.integration.generate import IntegrationCodeGenerator +from icon4py.liskov.codegen.integration.interface import IntegrationCodeInterface +from icon4py.liskov.codegen.serialisation.deserialise import ( + SerialisationCodeDeserialiser, +) +from icon4py.liskov.codegen.serialisation.generate import ( + SerialisationCodeGenerator, +) +from icon4py.liskov.codegen.shared.write import CodegenWriter from icon4py.liskov.external.gt4py import UpdateFieldsWithGt4PyStencils -from icon4py.liskov.parsing.deserialise import DirectiveDeserialiser from icon4py.liskov.parsing.parse import DirectivesParser from icon4py.liskov.parsing.scan import DirectivesScanner +from icon4py.liskov.pipeline.definition import Step, linear_pipeline + + +DESERIALISERS = { + "integration": IntegrationCodeDeserialiser, + "serialisation": SerialisationCodeDeserialiser, +} + +CODEGENS = { + "integration": IntegrationCodeGenerator, + "serialisation": SerialisationCodeGenerator, +} @linear_pipeline -def parse_fortran_file(input_filepath: Path, output_filepath: Path) -> list[Step]: +def parse_fortran_file( + input_filepath: Path, + output_filepath: Path, + deserialiser_type: str, + **kwargs, +) -> list[Step]: """Execute a pipeline to parse and deserialize directives from a file. The pipeline consists of three steps: DirectivesScanner, DirectivesParser, and DirectiveDeserialiser. The DirectivesScanner scans the file for directives, the DirectivesParser parses the directives into a dictionary, and the DirectiveDeserialiser deserializes the dictionary into a - DeserialisedDirectives object. + its corresponding Interface object. Args: input_filepath: Path to the input file to process. output_filepath: Path to the output file to generate. + deserialiser_type: What deserialiser to use. Returns: - DeserialisedDirectives: The deserialized directives object. + IntegrationCodeInterface | SerialisationCodeInterface: The interface object. """ + deserialiser = DESERIALISERS[deserialiser_type] + return [ DirectivesScanner(input_filepath), DirectivesParser(input_filepath, output_filepath), - DirectiveDeserialiser(), + deserialiser(**kwargs), ] @linear_pipeline -def load_gt4py_stencils(parsed: DeserialisedDirectives) -> list[Step]: - """Execute a pipeline to update fields of a DeserialisedDirectives object with GT4Py stencils. +def load_gt4py_stencils(parsed: IntegrationCodeInterface) -> list[Step]: + """Execute a pipeline to update fields of a IntegrationCodeInterface object with GT4Py stencils. Args: - parsed: The input DeserialisedDirectives object. + parsed: The input IntegrationCodeInterface object. Returns: The updated object with fields containing information from GT4Py stencils. @@ -62,26 +88,25 @@ def load_gt4py_stencils(parsed: DeserialisedDirectives) -> list[Step]: @linear_pipeline def run_code_generation( - parsed: DeserialisedDirectives, input_filepath: Path, output_filepath: Path, - profile: bool, - metadatagen: bool, + codegen_type: str, + *args, + **kwargs, ) -> list[Step]: - """Execute a pipeline to generate and write code for a set of directives. - - The pipeline consists of two steps: IntegrationGenerator and IntegrationWriter. The IntegrationGenerator generates - code based on the parsed directives and profile flag. The IntegrationWriter writes the generated code to the - specified filepath. + """Execute a pipeline to generate and write code. Args: - parsed: The deserialized directives object. input_filepath: The original file containing the DSL preprocessor directives. output_filepath: The file path to write the generated code to. - profile: A flag to indicate if profiling information should be included in the generated code. - metadatagen: A flag to indicate if a metadata header should be included in the generated code. + codegen_type: Which type of code generator to use. + + Note: + Additional positional and keyword arguments are passed to the code generator. """ + code_generator = CODEGENS[codegen_type] + return [ - IntegrationGenerator(parsed, profile, metadatagen), - IntegrationWriter(input_filepath, output_filepath), + code_generator(*args, **kwargs), + CodegenWriter(input_filepath, output_filepath), ] diff --git a/liskov/src/icon4py/liskov/common.py b/liskov/src/icon4py/liskov/pipeline/definition.py similarity index 100% rename from liskov/src/icon4py/liskov/common.py rename to liskov/src/icon4py/liskov/pipeline/definition.py diff --git a/liskov/src/icon4py/liskov/py.typed b/liskov/src/icon4py/liskov/py.typed new file mode 100644 index 0000000000..e69de29bb2 diff --git a/liskov/tests/conftest.py b/liskov/tests/conftest.py index fd870ff584..81a94b4faf 100644 --- a/liskov/tests/conftest.py +++ b/liskov/tests/conftest.py @@ -16,9 +16,6 @@ import pytest from click.testing import CliRunner -import icon4py.liskov.parsing.types as ts -from icon4py.liskov.parsing.scan import DirectivesScanner - @pytest.fixture def make_f90_tmpfile(tmp_path) -> Path: @@ -40,15 +37,3 @@ def _make_f90_tmpfile(content: str): @pytest.fixture def cli(): return CliRunner() - - -def scan_for_directives(fpath: Path) -> list[ts.RawDirective]: - collector = DirectivesScanner(fpath) - return collector() - - -def insert_new_lines(fname: Path, lines: list[str]) -> None: - """Append new lines into file.""" - with open(fname, "a") as f: - for ln in lines: - f.write(f"{ln}\n") diff --git a/liskov/tests/test_cli.py b/liskov/tests/test_cli.py index 75a6114da5..6fdffde422 100644 --- a/liskov/tests/test_cli.py +++ b/liskov/tests/test_cli.py @@ -11,48 +11,59 @@ # # SPDX-License-Identifier: GPL-3.0-or-later +import itertools + import pytest -from samples.fortran_samples import ( + +from icon4py.liskov.cli import main +from icon4py.testutils.liskov_fortran_samples import ( CONSECUTIVE_STENCIL, FREE_FORM_STENCIL, MULTIPLE_STENCILS, NO_DIRECTIVES_STENCIL, + REPEATED_STENCILS, SINGLE_STENCIL, ) -from icon4py.liskov.cli import main - @pytest.fixture def outfile(tmp_path): return str(tmp_path / "gen.f90") -@pytest.mark.parametrize("file", [NO_DIRECTIVES_STENCIL]) -def test_cli_no_directives(make_f90_tmpfile, cli, file, outfile): - fpath = str(make_f90_tmpfile(content=file)) - result = cli.invoke(main, [fpath, outfile]) - assert result.exit_code == 0 +test_cases = [] + +files = [ + ("NO_DIRECTIVES", NO_DIRECTIVES_STENCIL), + ("SINGLE", SINGLE_STENCIL), + ("CONSECUTIVE", CONSECUTIVE_STENCIL), + ("FREE_FORM", FREE_FORM_STENCIL), + ("MULTIPLE", MULTIPLE_STENCILS), + ("REPEATED", REPEATED_STENCILS), +] + +flags = {"serialise": ["--multinode"], "integrate": ["-p", "-m"]} + +for file_name, file_content in files: + for cmd in flags.keys(): + flag_combinations = [] + for r in range(1, len(flags[cmd]) + 1): + flag_combinations.extend(itertools.combinations(flags[cmd], r)) + for flags_selected in flag_combinations: + args = (file_name, file_content, cmd, list(flags_selected)) + test_cases.append(args) @pytest.mark.parametrize( - "file, profile", - [ - (NO_DIRECTIVES_STENCIL, False), - (SINGLE_STENCIL, False), - (CONSECUTIVE_STENCIL, False), - (FREE_FORM_STENCIL, False), - (MULTIPLE_STENCILS, False), - (SINGLE_STENCIL, True), - (CONSECUTIVE_STENCIL, True), - (FREE_FORM_STENCIL, True), - (MULTIPLE_STENCILS, True), + "file_name, file_content, cmd, cmd_flags", + test_cases, + ids=[ + "file={}, command={}, flags={}".format(file_name, cmd, ",".join(cmd_flags)) + for file_name, file_content, cmd, cmd_flags in test_cases ], ) -def test_cli(make_f90_tmpfile, cli, file, outfile, profile): - fpath = str(make_f90_tmpfile(content=file)) - args = [fpath, outfile] - if profile: - args.append("--profile") +def test_cli(make_f90_tmpfile, cli, outfile, file_name, file_content, cmd, cmd_flags): + fpath = str(make_f90_tmpfile(content=file_content)) + args = [fpath, outfile, cmd, *cmd_flags] result = cli.invoke(main, args) assert result.exit_code == 0 diff --git a/liskov/tests/test_deserialiser.py b/liskov/tests/test_directives_deserialiser.py similarity index 76% rename from liskov/tests/test_deserialiser.py rename to liskov/tests/test_directives_deserialiser.py index f32c4e2c04..bb10c24a1b 100644 --- a/liskov/tests/test_deserialiser.py +++ b/liskov/tests/test_directives_deserialiser.py @@ -15,8 +15,20 @@ import pytest -import icon4py.liskov.parsing.types as ts -from icon4py.liskov.codegen.interface import ( +import icon4py.liskov.parsing.parse +from icon4py.liskov.codegen.integration.deserialise import ( + DeclareDataFactory, + EndCreateDataFactory, + EndIfDataFactory, + EndProfileDataFactory, + EndStencilDataFactory, + ImportsDataFactory, + InsertDataFactory, + StartCreateDataFactory, + StartProfileDataFactory, + StartStencilDataFactory, +) +from icon4py.liskov.codegen.integration.interface import ( BoundsData, DeclareData, EndCreateData, @@ -29,18 +41,6 @@ StartCreateData, StartProfileData, ) -from icon4py.liskov.parsing.deserialise import ( - DeclareDataFactory, - EndCreateDataFactory, - EndIfDataFactory, - EndProfileDataFactory, - EndStencilDataFactory, - ImportsDataFactory, - InsertDataFactory, - StartCreateDataFactory, - StartProfileDataFactory, - StartStencilDataFactory, -) from icon4py.liskov.parsing.exceptions import ( DirectiveSyntaxError, MissingBoundsError, @@ -49,12 +49,40 @@ @pytest.mark.parametrize( - "factory_class, directive_type, startln, endln, string, expected", + "factory_class, directive_type, string, startln, endln, expected", [ - (EndCreateDataFactory, ts.EndCreate, "END CREATE", 2, 2, EndCreateData), - (ImportsDataFactory, ts.Imports, "IMPORTS", 3, 3, ImportsData), - (EndIfDataFactory, ts.EndIf, "ENDIF", 4, 4, EndIfData), - (EndProfileDataFactory, ts.EndProfile, "END PROFILE", 5, 5, EndProfileData), + ( + EndCreateDataFactory, + icon4py.liskov.parsing.parse.EndCreate, + "END CREATE", + 2, + 2, + EndCreateData, + ), + ( + ImportsDataFactory, + icon4py.liskov.parsing.parse.Imports, + "IMPORTS", + 3, + 3, + ImportsData, + ), + ( + EndIfDataFactory, + icon4py.liskov.parsing.parse.EndIf, + "ENDIF", + 4, + 4, + EndIfData, + ), + ( + EndProfileDataFactory, + icon4py.liskov.parsing.parse.EndProfile, + "END PROFILE", + 5, + 5, + EndProfileData, + ), ], ) def test_data_factories_no_args( @@ -72,7 +100,6 @@ def test_data_factories_no_args( assert isinstance(result, expected) assert result.startln == startln - assert result.endln == endln @pytest.mark.parametrize( @@ -83,8 +110,10 @@ def test_data_factories_no_args( EndStencilData, { "directives": [ - ts.EndStencil("END STENCIL(name=foo)", 5, 5), - ts.EndStencil( + icon4py.liskov.parsing.parse.EndStencil( + "END STENCIL(name=foo)", 5, 5 + ), + icon4py.liskov.parsing.parse.EndStencil( "END STENCIL(name=bar; noendif=true; noprofile=true)", 20, 20 ), ], @@ -101,7 +130,9 @@ def test_data_factories_no_args( EndStencilData, { "directives": [ - ts.EndStencil("END STENCIL(name=foo; noprofile=true)", 5, 5) + icon4py.liskov.parsing.parse.EndStencil( + "END STENCIL(name=foo; noprofile=true)", 5, 5 + ) ], "content": {"EndStencil": [{"name": "foo"}]}, }, @@ -111,8 +142,12 @@ def test_data_factories_no_args( StartProfileData, { "directives": [ - ts.StartProfile("START PROFILE(name=foo)", 5, 5), - ts.StartProfile("START PROFILE(name=bar)", 20, 20), + icon4py.liskov.parsing.parse.StartProfile( + "START PROFILE(name=foo)", 5, 5 + ), + icon4py.liskov.parsing.parse.StartProfile( + "START PROFILE(name=bar)", 20, 20 + ), ], "content": {"StartProfile": [{"name": "foo"}, {"name": "bar"}]}, }, @@ -121,7 +156,11 @@ def test_data_factories_no_args( StartProfileDataFactory, StartProfileData, { - "directives": [ts.StartProfile("START PROFILE(name=foo)", 5, 5)], + "directives": [ + icon4py.liskov.parsing.parse.StartProfile( + "START PROFILE(name=foo)", 5, 5 + ) + ], "content": {"StartProfile": [{"name": "foo"}]}, }, ), @@ -130,7 +169,7 @@ def test_data_factories_no_args( DeclareData, { "directives": [ - ts.Declare( + icon4py.liskov.parsing.parse.Declare( "DECLARE(vn=nlev,nblks_c; w=nlevp1,nblks_e; suffix=dsl; type=LOGICAL)", 5, 5, @@ -152,7 +191,9 @@ def test_data_factories_no_args( InsertDataFactory, InsertData, { - "directives": [ts.Insert("INSERT(content=foo)", 5, 5)], + "directives": [ + icon4py.liskov.parsing.parse.Insert("INSERT(content=foo)", 5, 5) + ], "content": {"Insert": ["foo"]}, }, ), @@ -169,7 +210,11 @@ def test_data_factories_with_args(factory, target, mock_data): [ ( { - "directives": [ts.StartCreate("START CREATE(extra_fields=foo)", 5, 5)], + "directives": [ + icon4py.liskov.parsing.parse.StartCreate( + "START CREATE(extra_fields=foo)", 5, 5 + ) + ], "content": {"StartCreate": [{"extra_fields": "foo"}]}, }, ["foo"], @@ -177,7 +222,9 @@ def test_data_factories_with_args(factory, target, mock_data): ( { "directives": [ - ts.StartCreate("START CREATE(extra_fields=foo,xyz)", 5, 5) + icon4py.liskov.parsing.parse.StartCreate( + "START CREATE(extra_fields=foo,xyz)", 5, 5 + ) ], "content": {"StartCreate": [{"extra_fields": "foo,xyz"}]}, }, @@ -185,7 +232,9 @@ def test_data_factories_with_args(factory, target, mock_data): ), ( { - "directives": [ts.StartCreate("START CREATE()", 5, 5)], + "directives": [ + icon4py.liskov.parsing.parse.StartCreate("START CREATE()", 5, 5) + ], "content": {"StartCreate": [None]}, }, None, @@ -207,8 +256,12 @@ def test_start_create_factory(mock_data, extra_fields): EndStencilData, { "directives": [ - ts.EndStencil("END STENCIL(name=foo)", 5, 5), - ts.EndStencil("END STENCIL(name=bar; noendif=foo)", 20, 20), + icon4py.liskov.parsing.parse.EndStencil( + "END STENCIL(name=foo)", 5, 5 + ), + icon4py.liskov.parsing.parse.EndStencil( + "END STENCIL(name=bar; noendif=foo)", 20, 20 + ), ], "content": { "EndStencil": [{"name": "foo"}, {"name": "bar", "noendif": "foo"}] diff --git a/liskov/tests/test_external.py b/liskov/tests/test_external.py index c97afed2d5..84a9acf2f3 100644 --- a/liskov/tests/test_external.py +++ b/liskov/tests/test_external.py @@ -17,9 +17,9 @@ import pytest from gt4py.next.ffront.decorator import Program -from icon4py.liskov.codegen.interface import ( - DeserialisedDirectives, +from icon4py.liskov.codegen.integration.interface import ( FieldAssociationData, + IntegrationCodeInterface, StartStencilData, ) from icon4py.liskov.external.exceptions import ( @@ -60,7 +60,7 @@ def test_stencil_collector_invalid_member(): os.remove(path) -mock_deserialised_directives = DeserialisedDirectives( +mock_deserialised_directives = IntegrationCodeInterface( StartStencil=[ StartStencilData( name="apply_nabla2_to_w", @@ -77,7 +77,6 @@ def test_stencil_collector_invalid_member(): ], bounds=None, startln=None, - endln=None, acc_present=False, mergecopy=False, copies=True, diff --git a/liskov/tests/test_generation.py b/liskov/tests/test_generation.py index a83f66b8e5..fb61064055 100644 --- a/liskov/tests/test_generation.py +++ b/liskov/tests/test_generation.py @@ -13,11 +13,10 @@ import pytest -from icon4py.liskov.codegen.generate import IntegrationGenerator -from icon4py.liskov.codegen.interface import ( +from icon4py.liskov.codegen.integration.generate import IntegrationCodeGenerator +from icon4py.liskov.codegen.integration.interface import ( BoundsData, DeclareData, - DeserialisedDirectives, EndCreateData, EndIfData, EndProfileData, @@ -25,15 +24,28 @@ FieldAssociationData, ImportsData, InsertData, + IntegrationCodeInterface, StartCreateData, StartProfileData, StartStencilData, ) - # TODO: fix tests to adapt to new custom output fields +from icon4py.liskov.codegen.serialisation.generate import ( + SerialisationCodeGenerator, +) +from icon4py.liskov.codegen.serialisation.interface import ( + FieldSerialisationData, + ImportData, + InitData, + Metadata, + SavepointData, + SerialisationCodeInterface, +) + + @pytest.fixture -def serialised_directives(): +def integration_code_interface(): start_stencil_data = StartStencilData( name="stencil1", fields=[ @@ -65,32 +77,28 @@ def serialised_directives(): ], bounds=BoundsData("1", "10", "-1", "-10"), startln=1, - endln=2, acc_present=False, mergecopy=False, copies=True, ) end_stencil_data = EndStencilData( - name="stencil1", startln=3, endln=4, noendif=False, noprofile=False + name="stencil1", startln=3, noendif=False, noprofile=False ) declare_data = DeclareData( startln=5, - endln=6, declarations={"field2": "(nproma, p_patch%nlev, p_patch%nblks_e)"}, ident_type="REAL(wp)", suffix="before", ) - imports_data = ImportsData(startln=7, endln=8) - start_create_data = StartCreateData( - extra_fields=["foo", "bar"], startln=9, endln=10 - ) - end_create_data = EndCreateData(startln=11, endln=11) - endif_data = EndIfData(startln=12, endln=12) - start_profile_data = StartProfileData(startln=13, endln=13, name="test_stencil") - end_profile_data = EndProfileData(startln=14, endln=14) - insert_data = InsertData(startln=15, endln=15, content="print *, 'Hello, World!'") + imports_data = ImportsData(startln=7) + start_create_data = StartCreateData(extra_fields=["foo", "bar"], startln=9) + end_create_data = EndCreateData(startln=11) + endif_data = EndIfData(startln=12) + start_profile_data = StartProfileData(startln=13, name="test_stencil") + end_profile_data = EndProfileData(startln=14) + insert_data = InsertData(startln=15, content="print *, 'Hello, World!'") - return DeserialisedDirectives( + return IntegrationCodeInterface( StartStencil=[start_stencil_data], EndStencil=[end_stencil_data], Declare=[declare_data], @@ -201,12 +209,14 @@ def expected_insert_source(): @pytest.fixture -def generator(serialised_directives): - return IntegrationGenerator(serialised_directives, profile=True, metadata_gen=False) +def integration_code_generator(integration_code_interface): + return IntegrationCodeGenerator( + integration_code_interface, profile=True, metadatagen=False + ) -def test_generate( - generator, +def test_integration_code_generation( + integration_code_generator, expected_start_create_source, expected_end_create_source, expected_imports_source, @@ -219,7 +229,7 @@ def test_generate( expected_insert_source, ): # Check that the generated code snippets are as expected - generated = generator() + generated = integration_code_generator() assert len(generated) == 10 assert generated[0].source == expected_start_create_source assert generated[1].source == expected_end_create_source @@ -231,3 +241,107 @@ def test_generate( assert generated[7].source == expected_start_profile_source assert generated[8].source == expected_end_profile_source assert generated[9].source == expected_insert_source + + +@pytest.fixture +def serialisation_code_interface(): + interface = { + "Import": ImportData(startln=0), + "Init": InitData(startln=1, directory=".", prefix="liskov-serialisation"), + "Savepoint": [ + SavepointData( + startln=9, + subroutine="apply_nabla2_to_vn_in_lateral_boundary", + intent="start", + fields=[ + FieldSerialisationData( + variable="z_nabla2_e", + association="z_nabla2_e(:,:,1)", + decomposed=False, + dimension=None, + typespec=None, + typename=None, + ptr_var=None, + ), + ], + metadata=[ + Metadata(key="jstep", value="jstep_ptr"), + Metadata(key="diffctr", value="diffctr"), + ], + ), + SavepointData( + startln=38, + subroutine="apply_nabla2_to_vn_in_lateral_boundary", + intent="end", + fields=[ + FieldSerialisationData( + variable="z_nabla2_e", + association="z_nabla2_e(:,:,1)", + decomposed=False, + dimension=None, + typespec=None, + typename=None, + ptr_var=None, + ), + FieldSerialisationData( + variable="vn", + association="p_nh_prog%vn(:,:,1)", + decomposed=False, + dimension=None, + typespec=None, + typename=None, + ptr_var=None, + ), + ], + metadata=[ + Metadata(key="jstep", value="jstep_ptr"), + Metadata(key="diffctr", value="diffctr"), + ], + ), + ], + } + + return SerialisationCodeInterface(**interface) + + +@pytest.fixture +def expected_savepoints(): + return [ + " USE mo_mpi, ONLY: get_my_mpi_work_id", + """ + !$ser init directory="." prefix="liskov-serialisation" mpi_rank=get_my_mpi_work_id() + + !$ser savepoint apply_nabla2_to_vn_in_lateral_boundary_start jstep=jstep_ptr diffctr=diffctr + + PRINT *, 'Serializing z_nabla2_e=z_nabla2_e(:,:,1)' + + !$ser data z_nabla2_e=z_nabla2_e(:,:,1)""", + """ + !$ser savepoint apply_nabla2_to_vn_in_lateral_boundary_end jstep=jstep_ptr diffctr=diffctr + + PRINT *, 'Serializing z_nabla2_e=z_nabla2_e(:,:,1)' + + !$ser data z_nabla2_e=z_nabla2_e(:,:,1) + + PRINT *, 'Serializing vn=p_nh_prog%vn(:,:,1)' + + !$ser data vn=p_nh_prog%vn(:,:,1)""", + ] + + +@pytest.mark.parametrize("multinode", [False, True]) +def test_serialisation_code_generation( + serialisation_code_interface, expected_savepoints, multinode +): + generated = SerialisationCodeGenerator( + serialisation_code_interface, multinode=multinode + )() + + if multinode: + assert len(generated) == 3 + assert generated[0].source == expected_savepoints[0] + assert generated[1].source == expected_savepoints[1] + assert generated[2].source == expected_savepoints[2] + else: + assert len(generated) == 2 + assert generated[1].source == expected_savepoints[2] diff --git a/liskov/tests/test_parser.py b/liskov/tests/test_parser.py index 8101988c51..6fa5095de9 100644 --- a/liskov/tests/test_parser.py +++ b/liskov/tests/test_parser.py @@ -15,17 +15,21 @@ from collections import defaultdict import pytest -from conftest import insert_new_lines, scan_for_directives from pytest import mark -from samples.fortran_samples import ( - MULTIPLE_STENCILS, - NO_DIRECTIVES_STENCIL, - SINGLE_STENCIL, -) +import icon4py.liskov.parsing.parse import icon4py.liskov.parsing.types as ts from icon4py.liskov.parsing.exceptions import UnsupportedDirectiveError from icon4py.liskov.parsing.parse import DirectivesParser +from icon4py.testutils.liskov_fortran_samples import ( + MULTIPLE_STENCILS, + NO_DIRECTIVES_STENCIL, + SINGLE_STENCIL, +) +from icon4py.testutils.liskov_test_utils import ( + insert_new_lines, + scan_for_directives, +) def test_parse_no_input(): @@ -37,21 +41,23 @@ def test_parse_no_input(): "directive, string, startln, endln, expected_content", [ ( - ts.Imports("IMPORTS()", 1, 1), + icon4py.liskov.parsing.parse.Imports("IMPORTS()", 1, 1), "IMPORTS()", 1, 1, defaultdict(list, {"Imports": [{}]}), ), ( - ts.StartCreate("START CREATE(extra_fields=foo)", 2, 2), + icon4py.liskov.parsing.parse.StartCreate( + "START CREATE(extra_fields=foo)", 2, 2 + ), "START CREATE()", 2, 2, defaultdict(list, {"StartCreate": [{"extra_fields": "foo"}]}), ), ( - ts.StartStencil( + icon4py.liskov.parsing.parse.StartStencil( "START STENCIL(name=mo_nh_diffusion_06; vn=p_patch%p%vn; foo=abc)", 3, 4 ), "START STENCIL(name=mo_nh_diffusion_06; vn=p_patch%p%vn; foo=abc)", diff --git a/liskov/tests/test_scanner.py b/liskov/tests/test_scanner.py index 64dcc4aedf..6e7c5c2252 100644 --- a/liskov/tests/test_scanner.py +++ b/liskov/tests/test_scanner.py @@ -16,11 +16,14 @@ import pytest from pytest import mark -from samples.fortran_samples import DIRECTIVES_SAMPLE, NO_DIRECTIVES_STENCIL from icon4py.liskov.parsing.exceptions import DirectiveSyntaxError from icon4py.liskov.parsing.scan import DirectivesScanner from icon4py.liskov.parsing.types import RawDirective +from icon4py.testutils.liskov_fortran_samples import ( + DIRECTIVES_SAMPLE, + NO_DIRECTIVES_STENCIL, +) ALLOWED_EOL_CHARS = [")", "&"] diff --git a/liskov/tests/test_serialisation_deserialiser.py b/liskov/tests/test_serialisation_deserialiser.py new file mode 100644 index 0000000000..5fe597de65 --- /dev/null +++ b/liskov/tests/test_serialisation_deserialiser.py @@ -0,0 +1,133 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import pytest + +from icon4py.liskov.codegen.serialisation.deserialise import ( + InitDataFactory, + SavepointDataFactory, +) +from icon4py.liskov.codegen.serialisation.interface import ( + FieldSerialisationData, + InitData, + Metadata, + SavepointData, +) +from icon4py.liskov.parsing.parse import ( + Declare, + EndCreate, + EndProfile, + EndStencil, + Imports, + StartCreate, + StartProfile, + StartStencil, +) + + +@pytest.fixture +def parsed_dict(): + parsed = { + "directives": [ + Imports(string="IMPORTS()", startln=0, endln=0), + StartCreate( + string="START CREATE()", + startln=2, + endln=2, + ), + Declare( + string="DECLARE(vn=nproma,p_patch%nlev,p_patch%nblks_e; suffix=dsl)", + startln=4, + endln=4, + ), + Declare( + string="DECLARE(vn= nproma,p_patch%nlev,p_patch%nblks_e; a=nproma,p_patch%nlev,p_patch%nblks_e; b=nproma,p_patch%nlev,p_patch%nblks_e; type=REAL(vp))", + startln=6, + endln=7, + ), + StartStencil( + string="START STENCIL(name=apply_nabla2_to_vn_in_lateral_boundary; z_nabla2_e=z_nabla2_e(:, :, 1); area_edge=p_patch%edges%area_edge(:,1); fac_bdydiff_v=fac_bdydiff_v; vn=p_nh_prog%vn(:,:,1); vertical_lower=1; vertical_upper=nlev; horizontal_lower=i_startidx; horizontal_upper=i_endidx; accpresent=True)", + startln=9, + endln=14, + ), + StartProfile( + string="START PROFILE(name=apply_nabla2_to_vn_in_lateral_boundary)", + startln=35, + endln=35, + ), + EndProfile( + string="END PROFILE()", + startln=37, + endln=37, + ), + EndStencil( + string="END STENCIL(name=apply_nabla2_to_vn_in_lateral_boundary; noprofile=True)", + startln=38, + endln=38, + ), + EndCreate( + string="END CREATE()", + startln=39, + endln=39, + ), + ], + "content": { + "Imports": [{}], + "StartCreate": [None], + "Declare": [ + {"vn": "nproma,p_patch%nlev,p_patch%nblks_e", "suffix": "dsl"}, + { + "vn": "nproma,p_patch%nlev,p_patch%nblks_e", + "a": "nproma,p_patch%nlev,p_patch%nblks_e", + "b": "nproma,p_patch%nlev,p_patch%nblks_e", + "type": "REAL(vp)", + }, + ], + "StartStencil": [ + { + "name": "apply_nabla2_to_vn_in_lateral_boundary", + "z_nabla2_e": "z_nabla2_e(:,:,1)", + "area_edge": "p_patch%edges%area_edge(:,1)", + "fac_bdydiff_v": "fac_bdydiff_v", + "vn": "p_nh_prog%vn(:,:,1)", + "vertical_lower": "1", + "vertical_upper": "nlev", + "horizontal_lower": "i_startidx", + "horizontal_upper": "i_endidx", + "accpresent": "True", + } + ], + "StartProfile": [{"name": "apply_nabla2_to_vn_in_lateral_boundary"}], + "EndProfile": [{}], + "EndStencil": [ + {"name": "apply_nabla2_to_vn_in_lateral_boundary", "noprofile": "True"} + ], + "EndCreate": [{}], + }, + } + return parsed + + +def test_init_data_factory(parsed_dict): + init = InitDataFactory()(parsed_dict) + assert isinstance(init, InitData) + assert init.directory == "." + assert init.prefix == "liskov-serialisation" + + +def test_savepoint_data_factory(parsed_dict): + savepoints = SavepointDataFactory()(parsed_dict) + assert len(savepoints) == 2 + assert any([isinstance(sp, SavepointData) for sp in savepoints]) + assert any([isinstance(f, FieldSerialisationData) for f in savepoints[0].fields]) + assert any([isinstance(m, Metadata) for m in savepoints[0].metadata]) diff --git a/liskov/tests/test_utils.py b/liskov/tests/test_utils.py index 1a327af09d..bd3d3313ae 100644 --- a/liskov/tests/test_utils.py +++ b/liskov/tests/test_utils.py @@ -14,8 +14,8 @@ import pytest -import icon4py.liskov.parsing.types as ts -from icon4py.liskov.parsing.types import Imports, StartCreate +import icon4py.liskov.parsing.parse +from icon4py.liskov.parsing.parse import Imports, StartCreate from icon4py.liskov.parsing.utils import ( extract_directive, print_parsed_directive, @@ -64,6 +64,6 @@ def test_string_to_bool(string, expected): def test_print_parsed_directive(): - directive = ts.Imports("IMPORTS()", 1, 1) + directive = icon4py.liskov.parsing.parse.Imports("IMPORTS()", 1, 1) expected_output = "Directive: IMPORTS(), start line: 1, end line: 1\n" assert print_parsed_directive(directive) == expected_output diff --git a/liskov/tests/test_validation.py b/liskov/tests/test_validation.py index 03bd6e5b53..ba6488cfa4 100644 --- a/liskov/tests/test_validation.py +++ b/liskov/tests/test_validation.py @@ -12,9 +12,7 @@ # SPDX-License-Identifier: GPL-3.0-or-later import pytest -from conftest import insert_new_lines, scan_for_directives from pytest import mark -from samples.fortran_samples import MULTIPLE_STENCILS, SINGLE_STENCIL from icon4py.liskov.parsing.exceptions import ( DirectiveSyntaxError, @@ -22,14 +20,22 @@ RequiredDirectivesError, UnbalancedStencilDirectiveError, ) -from icon4py.liskov.parsing.parse import DirectivesParser -from icon4py.liskov.parsing.types import ( +from icon4py.liskov.parsing.parse import ( Declare, + DirectivesParser, Imports, StartCreate, StartStencil, ) from icon4py.liskov.parsing.validation import DirectiveSyntaxValidator +from icon4py.testutils.liskov_fortran_samples import ( + MULTIPLE_STENCILS, + SINGLE_STENCIL, +) +from icon4py.testutils.liskov_test_utils import ( + insert_new_lines, + scan_for_directives, +) @mark.parametrize( diff --git a/liskov/tests/test_writer.py b/liskov/tests/test_writer.py index 2a91077910..41db096973 100644 --- a/liskov/tests/test_writer.py +++ b/liskov/tests/test_writer.py @@ -14,8 +14,8 @@ from pathlib import Path from tempfile import TemporaryDirectory -from icon4py.liskov.codegen.generate import GeneratedCode -from icon4py.liskov.codegen.write import DIRECTIVE_IDENT, IntegrationWriter +from icon4py.liskov.codegen.shared.types import GeneratedCode +from icon4py.liskov.codegen.shared.write import DIRECTIVE_IDENT, CodegenWriter def test_write_from(): @@ -28,8 +28,8 @@ def test_write_from(): f.write("!$DSL\n some code\n another line") # create an instance of IntegrationWriter and write generated code - generated = [GeneratedCode("generated code", 1, 3)] - integration_writer = IntegrationWriter(input_filepath, output_filepath) + generated = [GeneratedCode(1, "generated code")] + integration_writer = CodegenWriter(input_filepath, output_filepath) integration_writer(generated) # check that the generated code was inserted into the file @@ -49,14 +49,14 @@ def test_remove_directives(): "!$DSL another directive", ] expected_output = ["some code", "another line"] - assert IntegrationWriter._remove_directives(current_file) == expected_output + assert CodegenWriter._remove_directives(current_file) == expected_output def test_insert_generated_code(): current_file = ["some code", "another line"] generated = [ - GeneratedCode("generated code2", 5, 6), - GeneratedCode("generated code1", 1, 3), + GeneratedCode(5, "generated code2"), + GeneratedCode(1, "generated code1"), ] expected_output = [ "some code", @@ -65,8 +65,7 @@ def test_insert_generated_code(): "generated code2\n", ] assert ( - IntegrationWriter._insert_generated_code(current_file, generated) - == expected_output + CodegenWriter._insert_generated_code(current_file, generated) == expected_output ) @@ -77,7 +76,7 @@ def test_write_file(): output_filepath = input_filepath.with_suffix(".gen") generated_code = ["some code", "another line"] - writer = IntegrationWriter(input_filepath, output_filepath) + writer = CodegenWriter(input_filepath, output_filepath) writer._write_file(generated_code) # check that the generated code was written to the file diff --git a/pyutils/README.md b/pyutils/README.md index 931007a2b4..d57eed8d3f 100644 --- a/pyutils/README.md +++ b/pyutils/README.md @@ -1,21 +1,34 @@ # icon4py-pyutils -## Description +## icon4pygen -Python utilities for ICON4Py. +The `icon4pygen` tool generates GridTools C++ code as well as Fortran and C++ bindings for an icon4py fencil, so that it can be executed from within ICON. -### icongen +### Usage: -Generates from GT4Py Field View programs +`icon4pygen [OPTIONS] FENCIL BLOCK_SIZE LEVELS_PER_THREAD OUTPATH` -- GridTools C++ `gridtools::fn` code, -- field metadata to be used in dawn bindings generator for GT4Py. +#### Arguments: -## Installation instructions +``` +FENCIL: The fencil to generate code for. It can be specified as :, where is the dotted name of the containing module and is the name of the fencil. +BLOCK_SIZE: The number of threads per block to use in a cuda kernel. +LEVELS_PER_THREAD: How many k-levels to process per thread. +OUTPATH: The folder in which to write all generated code. +``` + +#### Options + +``` +--is_global: Flag to indicate if the stencil is global. +--imperative: Flag to indicate if the generated code should be written imperatively. +``` + +#### Example: -Check `README.md` file in the root of the repository. +`icon4pygen icon4py.atm_dyn_iconam.mo_velocity_advection_stencil_1:mo_velocity_advection_stencil_1 128 4 /path/to/output/folder` -## Autocomplete +#### Autocomplete In order to turn on autocomplete in your shell for `icon4pygen` you need to execute the following in your shell: @@ -24,3 +37,28 @@ eval "$(_ICON4PYGEN_COMPLETE=bash_source icon4pygen)" ``` To permanently enable autocomplete on your system add the above statement to your `~/.bashrc` file. + +## f2ser + +This tool is designed to parse a well-defined Fortran granule interface and generate ppser statements for each variable in the interface. It uses the `f2py` library to perform the parsing and `liskov` for the generation tasks. + +### Usage + +`f2ser [OPTIONS] GRANULE_PATH OUTPUT_FILEPATH` + +### Arguments + +``` +GRANULE_PATH A path to the Fortran source file to be parsed. +OUTPUT_FILEPATH A path to the output Fortran source file to be generated. +``` + +### Options + +``` +--dependencies PATH Optional list of dependency paths. +--directory TEXT The directory to serialise the variables to. +--prefix TEXT The prefix to use for each serialised variable. +``` + +**Note:** The output of f2ser still has to be preprocessed using `pp_ser.py`, which then yields a compilable unit. The serialised files will have `f2ser` as their prefix in the default folder location of the experiment. diff --git a/pyutils/setup.cfg b/pyutils/setup.cfg index de7bddd913..99f1fb95e6 100644 --- a/pyutils/setup.cfg +++ b/pyutils/setup.cfg @@ -53,4 +53,5 @@ exclude = [options.entry_points] console_scripts = + f2ser = icon4py.f2ser.cli:main icon4pygen = icon4py.icon4pygen.cli:main diff --git a/pyutils/src/icon4py/f2ser/__init__.py b/pyutils/src/icon4py/f2ser/__init__.py new file mode 100644 index 0000000000..15dfdb0098 --- /dev/null +++ b/pyutils/src/icon4py/f2ser/__init__.py @@ -0,0 +1,12 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later diff --git a/pyutils/src/icon4py/f2ser/cli.py b/pyutils/src/icon4py/f2ser/cli.py new file mode 100644 index 0000000000..8a65333249 --- /dev/null +++ b/pyutils/src/icon4py/f2ser/cli.py @@ -0,0 +1,79 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import pathlib +from typing import Optional + +import click + +from icon4py.f2ser.deserialise import ParsedGranuleDeserialiser +from icon4py.f2ser.parse import GranuleParser +from icon4py.liskov.codegen.serialisation.generate import ( + SerialisationCodeGenerator, +) +from icon4py.liskov.codegen.shared.write import CodegenWriter + + +@click.command("icon_f2ser") +@click.argument( + "granule_path", + type=click.Path( + exists=True, dir_okay=False, resolve_path=True, path_type=pathlib.Path + ), +) +@click.argument( + "output_filepath", + type=click.Path(dir_okay=False, resolve_path=True, path_type=pathlib.Path), +) +@click.option( + "--dependencies", + "-d", + multiple=True, + type=click.Path(exists=True), + help="Optional list of dependency paths.", +) +@click.option( + "--directory", type=str, help="Directory to serialise variables to.", default="." +) +@click.option( + "--prefix", type=str, help="Prefix to use for serialised files.", default="f2ser" +) +@click.option( + "--multinode", + is_flag=True, + type=bool, + help="Specify whether it is a multinode run.", + default=False, +) +def main( + granule_path: pathlib.Path, + dependencies: Optional[list[pathlib.Path]], + output_filepath: pathlib.Path, + directory: str, + prefix: str, + multinode: bool, +) -> None: + """Command line interface for interacting with the ICON-f2ser serialization Preprocessor. + + Arguments: + GRANULE_PATH: A path to the Fortran source file to be parsed. + OUTPUT_FILEPATH: A path to the output Fortran source file to be generated. + """ + parsed = GranuleParser(granule_path, dependencies)() + interface = ParsedGranuleDeserialiser(parsed, directory=directory, prefix=prefix)() + generated = SerialisationCodeGenerator(interface, multinode=multinode)() + CodegenWriter(granule_path, output_filepath)(generated) + + +if __name__ == "__main__": + main() diff --git a/pyutils/src/icon4py/f2ser/deserialise.py b/pyutils/src/icon4py/f2ser/deserialise.py new file mode 100644 index 0000000000..daeee81950 --- /dev/null +++ b/pyutils/src/icon4py/f2ser/deserialise.py @@ -0,0 +1,175 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from icon4py.f2ser.parse import CodegenContext, ParsedGranule +from icon4py.liskov.codegen.serialisation.interface import ( + FieldSerialisationData, + ImportData, + InitData, + SavepointData, + SerialisationCodeInterface, +) + + +class ParsedGranuleDeserialiser: + def __init__( + self, + parsed: ParsedGranule, + directory: str = ".", + prefix: str = "f2ser", + ): + self.parsed = parsed + self.directory = directory + self.prefix = prefix + self.data = {"Savepoint": [], "Init": ..., "Import": ...} + + def __call__(self) -> SerialisationCodeInterface: + """Deserialise the parsed granule and returns a serialisation interface. + + Returns: + A `SerialisationInterface` object representing the deserialised data. + """ + self._merge_out_inout_fields() + self._make_savepoints() + self._make_init_data() + self._make_imports() + return SerialisationCodeInterface(**self.data) + + def _make_savepoints(self) -> None: + """Create savepoints for each subroutine and intent in the parsed granule. + + Returns: + None. + """ + for subroutine_name, intent_dict in self.parsed.subroutines.items(): + for intent, var_dict in intent_dict.items(): + self._create_savepoint(subroutine_name, intent, var_dict) + + def _create_savepoint( + self, subroutine_name: str, intent: str, var_dict: dict + ) -> None: + """Create a savepoint for the given variables. + + Args: + subroutine_name: The name of the subroutine. + intent: The intent of the fields to be serialised. + var_dict: A dictionary representing the variables to be saved. + + Returns: + None. + """ + field_vals = {k: v for k, v in var_dict.items() if isinstance(v, dict)} + fields = [ + FieldSerialisationData( + variable=var_name, + association=self._create_association(var_data, var_name), + decomposed=var_data["decomposed"] + if var_data.get("decomposed") + else False, + dimension=var_data.get("dimension"), + typespec=var_data.get("typespec"), + typename=var_data.get("typename"), + ptr_var=var_data.get("ptr_var"), + ) + for var_name, var_data in field_vals.items() + ] + + self.data["Savepoint"].append( + SavepointData( + subroutine=subroutine_name, + intent=intent, + startln=self._get_codegen_line(var_dict["codegen_ctx"], intent), + fields=fields, + metadata=None, + ) + ) + + @staticmethod + def get_slice_expression(var_name: str, dimension: str) -> str: + """Return a string representing a slice expression for a given variable name and dimension. + + Args: + var_name (str): The name of the variable. + dimension (str): The dimension of the variable. + + Returns: + str: A string representing a slice expression. + """ + idx = dimension.split()[-1].lstrip("-+") + return f"{var_name}({idx}:)" + + def _create_association(self, var_data: dict, var_name: str) -> str: + """Create an association between a variable and its data. + + Args: + var_data (dict): A dictionary containing information about the variable. + var_name (str): The name of the variable. + + Returns: + str: A string representing the association between the variable and its data. + """ + dimension = var_data.get("dimension") + if dimension is not None: + return ( + self.get_slice_expression(var_name, dimension[0]) + if ":" not in dimension + else f"{var_name}({','.join(dimension)})" + ) + return var_name + + def _make_init_data(self) -> None: + """Create an `InitData` object and sets it to the `Init` key in the `data` dictionary. + + Returns: + None. + """ + first_intent_in_subroutine = [ + var_dict + for intent_dict in self.parsed.subroutines.values() + for intent, var_dict in intent_dict.items() + if intent == "in" + ][0] + startln = self._get_codegen_line( + first_intent_in_subroutine["codegen_ctx"], "init" + ) + self.data["Init"] = InitData( + startln=startln, + directory=self.directory, + prefix=self.prefix, + ) + + def _merge_out_inout_fields(self): + """Merge the `inout` fields into the `in` and `out` fields in the `parsed` dictionary. + + Returns: + None. + """ + for _, intent_dict in self.parsed.subroutines.items(): + if "inout" in intent_dict: + intent_dict["in"].update(intent_dict["inout"]) + intent_dict["out"].update(intent_dict["inout"]) + del intent_dict["inout"] + + @staticmethod + def _get_codegen_line(ctx: CodegenContext, intent: str): + if intent == "in": + return ctx.last_declaration_ln + elif intent == "out": + return ctx.end_subroutine_ln + elif intent == "init": + return ctx.first_declaration_ln + else: + raise ValueError(f"Unrecognized intent: {intent}") + + def _make_imports(self): + self.data["Import"] = ImportData(startln=self.parsed.last_import_ln) diff --git a/pyutils/src/icon4py/f2ser/exceptions.py b/pyutils/src/icon4py/f2ser/exceptions.py new file mode 100644 index 0000000000..318edf9dde --- /dev/null +++ b/pyutils/src/icon4py/f2ser/exceptions.py @@ -0,0 +1,20 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + + +class MissingDerivedTypeError(Exception): + ... + + +class ParsingError(Exception): + ... diff --git a/pyutils/src/icon4py/f2ser/parse.py b/pyutils/src/icon4py/f2ser/parse.py new file mode 100644 index 0000000000..5080ce9ba4 --- /dev/null +++ b/pyutils/src/icon4py/f2ser/parse.py @@ -0,0 +1,324 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later +import re +from copy import deepcopy +from dataclasses import dataclass +from enum import Enum +from pathlib import Path +from typing import Optional + +from numpy.f2py.crackfortran import crackfortran + +from icon4py.f2ser.exceptions import MissingDerivedTypeError, ParsingError + + +def crack(path: Path) -> dict: + return crackfortran(path)[0] + + +class SubroutineType(Enum): + RUN = "_run" + INIT = "_init" + + +@dataclass +class CodegenContext: + first_declaration_ln: int + last_declaration_ln: int + end_subroutine_ln: int + + +ParsedSubroutines = dict[str, dict[str, dict[str, any] | CodegenContext]] + + +@dataclass +class ParsedGranule: + subroutines: ParsedSubroutines + last_import_ln: int + + +class GranuleParser: + """Parses a Fortran source file and extracts information about its subroutines and variables. + + Attributes: + granule (Path): A path to the Fortran source file to be parsed. + dependencies (Optional[list[Path]]): A list of paths to any additional Fortran source files that the input file depends on. + + Example usage: + parser = GranuleParser(Path("my_file.f90"), dependencies=[Path("common.f90"), Path("constants.f90")]) + parsed_types = parser() + """ + + def __init__( + self, granule: Path, dependencies: Optional[list[Path]] = None + ) -> None: + self.granule_path = granule + self.dependencies = dependencies + + def __call__(self) -> ParsedGranule: + """Parse the granule and return the parsed data.""" + subroutines = self.parse_subroutines() + last_import_ln = self.find_last_fortran_use_statement() + return ParsedGranule(subroutines=subroutines, last_import_ln=last_import_ln) + + def parse_subroutines(self): + subroutines = self._extract_subroutines(crack(self.granule_path)) + variables_grouped_by_intent = { + name: self._extract_intent_vars(routine) + for name, routine in subroutines.items() + } + intrinsic_type_vars, derived_type_vars = self._parse_types( + variables_grouped_by_intent + ) + combined_type_vars = self._combine_types(derived_type_vars, intrinsic_type_vars) + with_lines = self._update_with_codegen_lines(combined_type_vars) + return with_lines + + def _extract_subroutines(self, parsed: dict[str, any]) -> dict[str, any]: + """Extract the _init and _run subroutines from the parsed granule. + + Args: + parsed: A dictionary representing the parsed granule. + + Returns: + A dictionary containing the extracted _init and _run subroutines. + """ + subroutines = {} + for elt in parsed["body"]: + name = elt["name"] + if SubroutineType.RUN.value in name or SubroutineType.INIT.value in name: + subroutines[name] = elt + + if len(subroutines) != 2: + raise ParsingError( + f"Did not find _init and _run subroutines in {self.granule_path}" + ) + + return subroutines + + @staticmethod + def _extract_intent_vars(subroutine: dict) -> dict: + """Extract variables grouped by their intent. + + Args: + subroutine (dict): A dictionary representing the subroutine. + + Returns: + A dictionary representing variables grouped by their intent. + """ + intents = ["in", "inout", "out"] + result: dict = {} + for var in subroutine["vars"]: + var_intent = subroutine["vars"][var]["intent"] + common_intents = list(set(intents).intersection(var_intent)) + for intent in common_intents: + if intent not in result: + result[intent] = {} + result[intent][var] = subroutine["vars"][var] + return result + + def _parse_types(self, parsed: dict) -> tuple[dict, dict]: + """Parse the intrinsic and derived type variables of each subroutine and intent from a parsed granule dictionary. + + Args: + parsed (dict): A dictionary containing the parsed information of a granule. + + Returns: + tuple[dict, dict]: A tuple containing two dictionaries. The first one maps each subroutine and intent to a dictionary of intrinsic type variables. The second one maps each subroutine and intent to a dictionary of derived type variables. + """ + intrinsic_types: dict = {} + derived_types: dict = {} + + for subroutine, subroutine_vars in parsed.items(): + intrinsic_types[subroutine] = {} + derived_types[subroutine] = {} + + for intent, intent_vars in subroutine_vars.items(): + intrinsic_vars = {} + derived_vars = {} + + for var_name, var_dict in intent_vars.items(): + if var_dict["typespec"] != "type": + intrinsic_vars[var_name] = var_dict + else: + derived_vars[var_name] = var_dict + + intrinsic_types[subroutine][intent] = intrinsic_vars + derived_types[subroutine][intent] = derived_vars + + return intrinsic_types, self._parse_derived_types(derived_types) + + def _parse_derived_types(self, derived_types: dict) -> dict: + """Parse the derived types defined in the input dictionary by adding their type definitions. + + Args: + derived_types (dict): A dictionary containing the derived types. + + Returns: + dict: A dictionary containing the parsed derived types with their type definitions. + + Raises: + MissingDerivedTypeError: If the type definition for a derived type could not be found in any of the dependency files. + """ + derived_type_defs = {} + if self.dependencies: + for dep in self.dependencies: + parsed = crack(dep) + for block in parsed["body"]: + if block["block"] == "type": + derived_type_defs[block["name"]] = block["vars"] + + for _, subroutine_vars in derived_types.items(): + for _, intent_vars in subroutine_vars.items(): + for _, var in intent_vars.items(): + if var["typespec"] == "type": + typename = var["typename"] + if typename in derived_type_defs: + var["typedef"] = derived_type_defs[typename] + var["decomposed"] = True + else: + raise MissingDerivedTypeError( + f"Could not find type definition for TYPE: {typename} in dependency files: {self.dependencies}" + ) + + return self._decompose_derived_types(derived_types) + + @staticmethod + def _decompose_derived_types(derived_types: dict) -> dict: + """Decompose derived types into individual subtypes. + + Args: + derived_types (dict): A dictionary containing the derived types to be decomposed. + + Returns: + dict: A dictionary containing the decomposed derived types, with each subtype represented by a separate entry. + """ + decomposed_vars: dict = {} + for subroutine, subroutine_vars in derived_types.items(): + decomposed_vars[subroutine] = {} + for intent, intent_vars in subroutine_vars.items(): + decomposed_vars[subroutine][intent] = {} + for var_name, var_dict in intent_vars.items(): + if "typedef" in var_dict: + typedef = var_dict["typedef"] + del var_dict["typedef"] + for subtype_name, subtype_spec in typedef.items(): + new_type_name = f"{var_name}_{subtype_name}" + new_var_dict = var_dict.copy() + new_var_dict.update(subtype_spec) + decomposed_vars[subroutine][intent][ + new_type_name + ] = new_var_dict + new_var_dict["ptr_var"] = subtype_name + else: + decomposed_vars[subroutine][intent][var_name] = var_dict + + return decomposed_vars + + @staticmethod + def _combine_types(derived_type_vars: dict, intrinsic_type_vars: dict) -> dict: + """Combine intrinsic and derived type variables and returns a dictionary with the combined result. + + Args: + derived_type_vars (dict): A dictionary with derived type variables. + intrinsic_type_vars (dict): A dictionary with intrinsic type variables. + + Returns: + dict: A dictionary with the combined intrinsic and derived type variables. + """ + combined = deepcopy(intrinsic_type_vars) + for subroutine_name in combined: + for intent in combined[subroutine_name]: + new_vars = derived_type_vars[subroutine_name][intent] + combined[subroutine_name][intent].update(new_vars) + return combined + + def _update_with_codegen_lines(self, parsed_types: dict) -> dict: + """Update the parsed_types dictionary with the line numbers for codegen. + + Args: + parsed_types (dict): A dictionary containing the parsed intrinsic and derived types. + + Returns: + dict: A dictionary containing the parsed intrinsic and derived types with line numbers for codegen. + """ + with_lines = deepcopy(parsed_types) + for subroutine in with_lines: + for intent in with_lines[subroutine]: + with_lines[subroutine][intent][ + "codegen_ctx" + ] = self.get_subroutine_lines(subroutine) + return with_lines + + def find_last_fortran_use_statement(self): + with open(self.granule_path) as f: + file_contents = f.readlines() + + # Reverse the order of the lines so we can search from the end + file_contents.reverse() + + # Look for the last USE statement + use_ln = None + for i, line in enumerate(file_contents): + if line.strip().lower().startswith("use"): + use_ln = len(file_contents) - i + if i > 0 and file_contents[i - 1].strip().lower() == "#endif": + # If the USE statement is preceded by an #endif statement, return the line number after the #endif statement + return use_ln + 1 + else: + return use_ln + return None + + def get_subroutine_lines(self, subroutine_name: str) -> CodegenContext: + """Return CodegenContext object containing line numbers of the last declaration statement and the code before the end of the given subroutine. + + Args: + subroutine_name (str): Name of the subroutine to look for in the code. + + Returns: + CodegenContext: Object containing the line number of the last declaration statement and the line number of the last line of the code before the end of the given subroutine. + """ + with open(self.granule_path) as f: + code = f.read() + + # Find the line number where the subroutine is defined + start_subroutine_pattern = r"SUBROUTINE\s+" + subroutine_name + r"\s*\(" + end_subroutine_pattern = r"END\s+SUBROUTINE\s+" + subroutine_name + r"\s*" + start_match = re.search(start_subroutine_pattern, code) + end_match = re.search(end_subroutine_pattern, code) + if start_match is None or end_match is None: + return None + start_subroutine_ln = code[: start_match.start()].count("\n") + 1 + end_subroutine_ln = code[: end_match.start()].count("\n") + 1 + + # Find the last intent statement line number in the subroutine + declaration_pattern = r".*::\s*(\w+\b)" + declaration_pattern_lines = [ + i + for i, line in enumerate( + code.splitlines()[start_subroutine_ln:end_subroutine_ln] + ) + if re.search(declaration_pattern, line) + ] + if not declaration_pattern_lines: + raise ParsingError(f"No declarations found in {self.granule_path}") + last_declaration_ln = declaration_pattern_lines[-1] + start_subroutine_ln + 1 + first_declaration_ln = declaration_pattern_lines[0] + start_subroutine_ln + + pre_end_subroutine_ln = ( + end_subroutine_ln - 1 + ) # we want to generate the code before the end of the subroutine + + return CodegenContext( + first_declaration_ln, last_declaration_ln, pre_end_subroutine_ln + ) diff --git a/pyutils/src/icon4py/icon4pygen/bindings/codegen/render/location.py b/pyutils/src/icon4py/icon4pygen/bindings/codegen/render/location.py index b90f32d846..56dd9cf70c 100644 --- a/pyutils/src/icon4py/icon4pygen/bindings/codegen/render/location.py +++ b/pyutils/src/icon4py/icon4pygen/bindings/codegen/render/location.py @@ -15,7 +15,6 @@ class LocationRenderer: - type_dispatcher = {"Cell": "Cells", "Edge": "Edges", "Vertex": "Vertices"} @classmethod diff --git a/pyutils/src/icon4py/icon4pygen/cli.py b/pyutils/src/icon4py/icon4pygen/cli.py index 382726abf5..b27c529d31 100644 --- a/pyutils/src/icon4py/icon4pygen/cli.py +++ b/pyutils/src/icon4py/icon4pygen/cli.py @@ -45,13 +45,20 @@ def shell_complete(self, ctx, param, incomplete): @click.argument("fencil", type=ModuleType()) @click.argument("block_size", type=int, default=128) @click.argument("levels_per_thread", type=int, default=4) -@click.option("--is_global", is_flag=True, type=bool) +@click.option( + "--is_global", is_flag=True, type=bool, help="Whether this is a global run." +) @click.argument( "outpath", type=click.Path(dir_okay=True, resolve_path=True, path_type=pathlib.Path), default=".", ) -@click.option("--imperative", is_flag=True, type=bool) +@click.option( + "--imperative", + is_flag=True, + type=bool, + help="Whether to use the imperative code generation backend.", +) def main( fencil: str, block_size: int, @@ -63,15 +70,11 @@ def main( """ Generate Gridtools C++ code for an icon4py fencil as well as all the associated C++ and Fortran bindings. - Args: - fencil: may be specified as :, where is the dotted name of the containing module - and is the name of the fencil. - - block_size: refers to the number of threads per block to use in a cuda kernel. - - levels_per_thread: how many k-levels to process per thread. - - outpath: represents a path to the folder in which to write all generated code. + Arguments: + FENCIL: may be specified as :, where is the dotted name of the containing module and is the name of the fencil. + BLOCK_SIZE: refers to the number of threads per block to use in a cuda kernel. + LEVELS_PER_THREAD: how many k-levels to process per thread. + OUTPATH: represents a path to the folder in which to write all generated code. """ from icon4py.icon4pygen.backend import GTHeader from icon4py.icon4pygen.bindings.workflow import PyBindGen diff --git a/pyutils/tests/f2ser/conftest.py b/pyutils/tests/f2ser/conftest.py new file mode 100644 index 0000000000..ff73a38d4f --- /dev/null +++ b/pyutils/tests/f2ser/conftest.py @@ -0,0 +1,43 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from pathlib import Path + +import pytest + +import icon4py.testutils as testutils + + +@pytest.fixture +def samples_path(): + return Path(testutils.__file__).parent / "fortran" + + +@pytest.fixture +def diffusion_granule(samples_path): + return samples_path / "diffusion_granule.f90" + + +@pytest.fixture +def diffusion_granule_deps(samples_path): + return [samples_path / "derived_types_example.f90"] + + +@pytest.fixture +def no_deps_source_file(samples_path): + return samples_path / "no_deps_subroutine_example.f90" + + +@pytest.fixture +def not_existing_diffusion_granule(samples_path): + return samples_path / "not_existing_file.f90" diff --git a/pyutils/tests/f2ser/test_f2ser_cli.py b/pyutils/tests/f2ser/test_f2ser_cli.py new file mode 100644 index 0000000000..e4f27dcdbd --- /dev/null +++ b/pyutils/tests/f2ser/test_f2ser_cli.py @@ -0,0 +1,77 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import pytest +from click.testing import CliRunner + +from icon4py.f2ser.cli import main +from icon4py.f2ser.exceptions import MissingDerivedTypeError + + +@pytest.fixture +def outfile(tmp_path): + return str(tmp_path / "gen.f90") + + +@pytest.fixture +def cli(): + return CliRunner() + + +@pytest.mark.parametrize("multinode", [False, True]) +def test_cli(diffusion_granule, diffusion_granule_deps, outfile, cli, multinode): + inp = str(diffusion_granule) + deps = [str(p) for p in diffusion_granule_deps] + args = [inp, outfile, "-d", ",".join(deps)] + if multinode: + args.append("--multinode") + result = cli.invoke(main, args) + assert result.exit_code == 0 + + +def test_cli_no_deps(no_deps_source_file, outfile, cli): + inp = str(no_deps_source_file) + args = [inp, outfile] + result = cli.invoke(main, args) + assert result.exit_code == 0 + + +def test_cli_wrong_deps(diffusion_granule, samples_path, outfile, cli): + inp = str(diffusion_granule) + deps = [str(samples_path / "wrong_derived_types_example.f90")] + args = [inp, outfile, "-d", *deps] + result = cli.invoke(main, args) + assert result.exit_code == 2 + assert "Invalid value for '--dependencies' / '-d'" in result.output + + +def test_cli_missing_deps(diffusion_granule, outfile, cli): + inp = str(diffusion_granule) + args = [inp, outfile] + result = cli.invoke(main, args) + assert isinstance(result.exception, MissingDerivedTypeError) + + +def test_cli_wrong_source(outfile, cli): + inp = str("foo.90") + args = [inp, outfile] + result = cli.invoke(main, args) + assert "Invalid value for 'GRANULE_PATH'" in result.output + + +def test_cli_missing_source(not_existing_diffusion_granule, outfile, cli): + inp = str(not_existing_diffusion_granule) + args = [inp, outfile] + result = cli.invoke(main, args) + assert isinstance(result.exception, SystemExit) + assert "Invalid value for 'GRANULE_PATH'" in result.output diff --git a/pyutils/tests/f2ser/test_f2ser_codegen.py b/pyutils/tests/f2ser/test_f2ser_codegen.py new file mode 100644 index 0000000000..23604e7d93 --- /dev/null +++ b/pyutils/tests/f2ser/test_f2ser_codegen.py @@ -0,0 +1,110 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import pytest + +from icon4py.f2ser.deserialise import ParsedGranuleDeserialiser +from icon4py.f2ser.parse import GranuleParser +from icon4py.liskov.codegen.serialisation.generate import ( + SerialisationCodeGenerator, +) +from icon4py.liskov.codegen.shared.types import GeneratedCode + + +def test_deserialiser_diffusion_codegen(diffusion_granule, diffusion_granule_deps): + parsed = GranuleParser(diffusion_granule, diffusion_granule_deps)() + interface = ParsedGranuleDeserialiser(parsed)() + generated = SerialisationCodeGenerator(interface)() + assert len(generated) == 3 + + +@pytest.fixture +def expected_no_deps_serialization_directives(): + serialization_directives = [ + GeneratedCode( + startln=12, + source="\n" + ' !$ser init directory="." prefix="f2ser"\n' + "\n" + " !$ser savepoint no_deps_init_in\n" + "\n" + " PRINT *, 'Serializing a=a'\n" + "\n" + " !$ser data a=a\n" + "\n" + " PRINT *, 'Serializing b=b'\n" + "\n" + " !$ser data b=b", + ), + GeneratedCode( + startln=14, + source="\n" + " !$ser savepoint no_deps_init_out\n" + "\n" + " PRINT *, 'Serializing c=c'\n" + "\n" + " !$ser data c=c\n" + "\n" + " PRINT *, 'Serializing b=b'\n" + "\n" + " !$ser data b=b", + ), + GeneratedCode( + startln=20, + source="\n" + " !$ser savepoint no_deps_run_in\n" + "\n" + " PRINT *, 'Serializing a=a'\n" + "\n" + " !$ser data a=a\n" + "\n" + " PRINT *, 'Serializing b=b'\n" + "\n" + " !$ser data b=b", + ), + GeneratedCode( + startln=22, + source="\n" + " !$ser savepoint no_deps_run_out\n" + "\n" + " PRINT *, 'Serializing c=c'\n" + "\n" + " !$ser data c=c\n" + "\n" + " PRINT *, 'Serializing b=b'\n" + "\n" + " !$ser data b=b", + ), + ] + return serialization_directives + + +def test_deserialiser_directives_no_deps_codegen( + no_deps_source_file, expected_no_deps_serialization_directives +): + parsed = GranuleParser(no_deps_source_file)() + interface = ParsedGranuleDeserialiser(parsed)() + generated = SerialisationCodeGenerator(interface)() + assert generated == expected_no_deps_serialization_directives + + +def test_deserialiser_directives_diffusion_codegen( + diffusion_granule, diffusion_granule_deps, samples_path +): + parsed = GranuleParser(diffusion_granule, diffusion_granule_deps)() + interface = ParsedGranuleDeserialiser(parsed)() + generated = SerialisationCodeGenerator(interface)() + reference_savepoint = ( + samples_path / "expected_diffusion_granule_savepoint.f90" + ).read_text() + assert generated[0].source == reference_savepoint.rstrip() diff --git a/pyutils/tests/f2ser/test_granule_deserialiser.py b/pyutils/tests/f2ser/test_granule_deserialiser.py new file mode 100644 index 0000000000..8c8a60b912 --- /dev/null +++ b/pyutils/tests/f2ser/test_granule_deserialiser.py @@ -0,0 +1,99 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later +import pytest + +from icon4py.f2ser.deserialise import ParsedGranuleDeserialiser +from icon4py.f2ser.parse import CodegenContext, GranuleParser, ParsedGranule +from icon4py.liskov.codegen.serialisation.interface import ( + FieldSerialisationData, + SavepointData, + SerialisationCodeInterface, +) + + +@pytest.fixture +def mock_parsed_granule(): + return ParsedGranule( + subroutines={ + "diffusion_init": { + "in": { + "jg": {"typespec": "integer", "attrspec": [], "intent": ["in"]}, + "vt": { + "typespec": "real", + "kindselector": {"kind": "vp"}, + "attrspec": [], + "intent": ["in"], + "dimension": [":", ":", ":"], + }, + "codegen_ctx": CodegenContext(432, 450, 600), + } + }, + "diffusion_run": { + "out": { + "vert_idx": { + "typespec": "logical", + "kindselector": {"kind": "vp"}, + "attrspec": [], + "intent": ["in"], + "dimension": [":", ":", ":"], + }, + "codegen_ctx": CodegenContext(800, 850, 1000), + }, + "in": { + "vn": {"typespec": "integer", "attrspec": [], "intent": ["out"]}, + "vert_idx": { + "typespec": "logical", + "kindselector": {"kind": "vp"}, + "attrspec": [], + "intent": ["in"], + "dimension": [":", ":", ":"], + }, + "codegen_ctx": CodegenContext(600, 690, 750), + }, + "inout": { + "vn": {"typespec": "integer", "attrspec": [], "intent": ["out"]}, + "vert_idx": { + "typespec": "logical", + "kindselector": {"kind": "vp"}, + "attrspec": [], + "intent": ["in"], + "dimension": [":", ":", ":"], + }, + }, + }, + }, + last_import_ln=59, + ) + + +def test_deserialiser_mock(mock_parsed_granule): + deserialiser = ParsedGranuleDeserialiser(mock_parsed_granule) + interface = deserialiser() + assert isinstance(interface, SerialisationCodeInterface) + assert len(interface.Savepoint) == 3 + assert all([isinstance(s, SavepointData) for s in interface.Savepoint]) + assert all( + [ + isinstance(f, FieldSerialisationData) + for s in interface.Savepoint + for f in s.fields + ] + ) + + +def test_deserialiser_diffusion_granule(diffusion_granule, diffusion_granule_deps): + parser = GranuleParser(diffusion_granule, diffusion_granule_deps) + parsed = parser() + deserialiser = ParsedGranuleDeserialiser(parsed) + interface = deserialiser() + assert len(interface.Savepoint) == 3 diff --git a/pyutils/tests/f2ser/test_parsing.py b/pyutils/tests/f2ser/test_parsing.py new file mode 100644 index 0000000000..2b606262b5 --- /dev/null +++ b/pyutils/tests/f2ser/test_parsing.py @@ -0,0 +1,63 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import pytest + +from icon4py.f2ser.exceptions import MissingDerivedTypeError, ParsingError +from icon4py.f2ser.parse import CodegenContext, GranuleParser + + +def test_granule_parsing(diffusion_granule, diffusion_granule_deps): + parser = GranuleParser(diffusion_granule, diffusion_granule_deps) + parsed_granule = parser() + + subroutines = parsed_granule.subroutines + + assert list(subroutines) == ["diffusion_init", "diffusion_run"] + + assert list(subroutines["diffusion_init"]) == ["in"] + assert len(subroutines["diffusion_init"]["in"]) == 107 + assert subroutines["diffusion_init"]["in"]["codegen_ctx"] == CodegenContext( + first_declaration_ln=190, last_declaration_ln=280, end_subroutine_ln=401 + ) + + assert list(subroutines["diffusion_run"]) == ["in", "inout", "out"] + assert len(subroutines["diffusion_run"]["in"]) == 5 + assert subroutines["diffusion_run"]["in"]["codegen_ctx"] == CodegenContext( + first_declaration_ln=417, last_declaration_ln=492, end_subroutine_ln=1965 + ) + + assert len(subroutines["diffusion_run"]["inout"]) == 8 + + assert len(subroutines["diffusion_run"]["out"]) == 5 + assert subroutines["diffusion_run"]["out"]["codegen_ctx"] == CodegenContext( + first_declaration_ln=417, last_declaration_ln=492, end_subroutine_ln=1965 + ) + + assert isinstance(subroutines, dict) + assert parsed_granule.last_import_ln == 60 + + +def test_granule_parsing_missing_derived_typedef(diffusion_granule, samples_path): + dependencies = [samples_path / "subroutine_example.f90"] + parser = GranuleParser(diffusion_granule, dependencies) + with pytest.raises( + MissingDerivedTypeError, match="Could not find type definition for TYPE" + ): + parser() + + +def test_granule_parsing_no_intent(samples_path): + parser = GranuleParser(samples_path / "subroutine_example.f90", []) + with pytest.raises(ParsingError): + parser() diff --git a/testutils/src/icon4py/testutils/fortran/derived_types_example.f90 b/testutils/src/icon4py/testutils/fortran/derived_types_example.f90 new file mode 100644 index 0000000000..afbf542e4f --- /dev/null +++ b/testutils/src/icon4py/testutils/fortran/derived_types_example.f90 @@ -0,0 +1,150 @@ +!> +!! Contains basic math types +!! +!! @par Revision History +!! +!! @par Copyright and License +!! +!! This code is subject to the DWD and MPI-M-Software-License-Agreement in +!! its most recent form. +!! Please see the file LICENSE in the root of the source tree for this code. +!! Where software is supplied by third parties, it is indicated in the +!! headers of the routines. + +MODULE mo_math_types_base + USE ISO_C_BINDING, ONLY: C_INT64_T, C_DOUBLE + USE mo_kind_base, ONLY: wp, sp, dp + IMPLICIT NONE + + PRIVATE + + PUBLIC :: t_cartesian_coordinates + PUBLIC :: t_geographical_coordinates + PUBLIC :: t_line + PUBLIC :: t_tangent_vectors + PUBLIC :: t_Statistics + + ! cartesian coordinate class + TYPE, BIND(C) :: t_cartesian_coordinates + REAL(C_DOUBLE) :: x(3) + END TYPE t_cartesian_coordinates + + ! geographical coordinate class + TYPE, BIND(C) :: t_geographical_coordinates + REAL(C_DOUBLE) :: lon + REAL(C_DOUBLE) :: lat + END TYPE t_geographical_coordinates + + ! the two coordinates on the tangent plane + TYPE, BIND(C) :: t_tangent_vectors + REAL(C_DOUBLE) :: v1 + REAL(C_DOUBLE) :: v2 + END TYPE t_tangent_vectors + + ! line class + TYPE t_line + TYPE(t_geographical_coordinates) :: p1 + TYPE(t_geographical_coordinates) :: p2 + END TYPE t_line + + TYPE :: t_Statistics + INTEGER(C_INT64_T) :: sampleCount + REAL(wp) :: MIN, mean, MAX + CONTAINS + PROCEDURE :: reset => statistics_reset + GENERIC :: add => addData_s1d, addData_d1d, addStatistics + ! scan a given array AND update the statistics accordingly + PROCEDURE :: addData_s1d => statistics_addData_s1d + PROCEDURE :: addData_d1d => statistics_addData_d1d + PROCEDURE :: addStatistics => statistics_addStatistics ! update the statistics with the contents of another t_Statistics object + END TYPE t_Statistics + +CONTAINS + + SUBROUTINE statistics_reset(me) + CLASS(t_Statistics), INTENT(INOUT) :: me + + me%sampleCount = 0_C_INT64_T + me%MIN = HUGE(me%MIN) + me%mean = 0.0_wp + me%MAX = -HUGE(me%MAX) + END SUBROUTINE statistics_reset + + SUBROUTINE statistics_addData_s1d(me, DATA) + CLASS(t_Statistics), INTENT(INOUT) :: me + REAL(sp), INTENT(IN) :: DATA(:) + INTEGER :: i, icount + REAL(wp) :: data_sum, data_max, data_min, data_wp + + TYPE(t_Statistics) :: newStatistics + + CALL newStatistics%reset() + + icount = 0 + data_sum = 0._wp + data_max = -HUGE(DATA) + data_min = HUGE(DATA) +!$OMP PARALLEL DO PRIVATE(data_wp), REDUCTION(+:data_sum,icount), REDUCTION(MAX:data_max), REDUCTION(MIN:data_min) + DO i = 1, SIZE(DATA) + icount = icount+1 + data_wp = REAL(DATA(i), wp) + data_sum = data_sum + data_wp + data_max = MAX(data_max, data_wp) + data_min = MIN(data_min, data_wp) + ENDDO +!$OMP END PARALLEL DO + newStatistics%sampleCount = icount + newStatistics%MIN = data_min + newStatistics%MAX = data_max + IF (icount > 0) THEN + newStatistics%mean = data_sum / REAL(icount,wp) + ENDIF + CALL me%add(newStatistics) + END SUBROUTINE statistics_addData_s1d + + SUBROUTINE statistics_addData_d1d(me, DATA) + CLASS(t_Statistics), INTENT(INOUT) :: me + REAL(dp), INTENT(IN) :: DATA(:) + INTEGER :: i, icount + REAL(wp) :: data_sum, data_max, data_min + + TYPE(t_Statistics) :: newStatistics + + CALL newStatistics%reset() + + icount = 0 + data_sum = 0._wp + data_max = -HUGE(DATA) + data_min = HUGE(DATA) +!$OMP PARALLEL DO REDUCTION(+:data_sum,icount), REDUCTION(MAX:data_max), REDUCTION(MIN:data_min) + DO i = 1, SIZE(DATA) + icount = icount+1 + data_sum = data_sum + DATA(i) + data_max = MAX(data_max, DATA(i)) + data_min = MIN(data_min, DATA(i)) + ENDDO +!$OMP END PARALLEL DO + newStatistics%sampleCount = icount + newStatistics%MIN = data_min + newStatistics%MAX = data_max + IF (icount > 0) THEN + newStatistics%mean = data_sum / REAL(icount,wp) + ENDIF + + CALL me%add(newStatistics) + END SUBROUTINE statistics_addData_d1d + + SUBROUTINE statistics_addStatistics(me, other) + CLASS(t_Statistics), INTENT(INOUT) :: me + CLASS(t_Statistics), INTENT(IN) :: other + + INTEGER(C_INT64_T) :: newSampleCount + + newSampleCount = me%sampleCount + other%sampleCount + me%MIN = MIN(me%MIN, other%MIN) + me%mean = (me%mean*REAL(me%sampleCount, wp) + other%mean*REAL(other%sampleCount, wp))/REAL(newSampleCount, wp) + me%MAX = MAX(me%MAX, other%MAX) + me%sampleCount = newSampleCount + END SUBROUTINE statistics_addStatistics + +END MODULE mo_math_types_base diff --git a/testutils/src/icon4py/testutils/fortran/diffusion_granule.f90 b/testutils/src/icon4py/testutils/fortran/diffusion_granule.f90 new file mode 100644 index 0000000000..84b5768b2c --- /dev/null +++ b/testutils/src/icon4py/testutils/fortran/diffusion_granule.f90 @@ -0,0 +1,1981 @@ +!> +!! mo_nh_diffusion_new +!! +!! Diffusion in the nonhydrostatic model +!! +!! @author Almut Gassmann, MPI-M +!! +!! +!! @par Revision History +!! Initial release by Almut Gassmann, MPI-M (2009-08.25) +!! Modification by William Sawyer, CSCS (2015-02-06) +!! - OpenACC implementation +!! Modification by William Sawyer, CSCS (2015-02-06) +!! - Turned into a granule +!! +!! @par Copyright and License +!! +!! This code is subject to the DWD and MPI-M-Software-License-Agreement in +!! its most recent form. +!! Please see the file LICENSE in the root of the source tree for this code. +!! Where software is supplied by third parties, it is indicated in the +!! headers of the routines. +!! + +!---------------------------- +#include "omp_definitions.inc" +!---------------------------- + +MODULE mo_nh_diffusion_new + +#ifdef __SX__ +! for strange reasons, this routine is faster without mixed precision on the NEC +#undef __MIXED_PRECISION + USE mo_kind_base, ONLY: wp, vp => wp +#else + USE mo_kind_base, ONLY: wp, vp +#endif + USE mo_math_types, ONLY: t_tangent_vectors ! to maintain compatibility w/ p_patch +#if 0 + USE mo_math_types_base, ONLY: t_tangent_vectors +#endif + USE mo_model_domain_advanced, ONLY: t_patch ! until GridManager available + USE mo_model_domain, ONLY: p_patch + USE mo_intp_rbf_math, ONLY: rbf_vec_interpol_vertex, rbf_vec_interpol_cell + USE mo_interpolation_scalar_math, ONLY: cells2verts_scalar + USE mo_interpolation_vector_math, ONLY: edges2cells_vector + USE mo_loopindices_advanced, ONLY: get_indices_e, get_indices_c + USE mo_impl_constants_base , ONLY: min_rledge, min_rlcell, min_rlvert, min_rledge_int, min_rlcell_int, min_rlvert_int + USE mo_impl_constants_grf_base, ONLY: grf_bdywidth_e, grf_bdywidth_c + USE mo_math_types_base, ONLY: t_geographical_coordinates + USE mo_math_laplace_math, ONLY: nabla4_vec + USE mo_math_constants_base, ONLY: dbl_eps + USE mo_sync, ONLY: SYNC_E, SYNC_C, SYNC_V, sync_patch_array, & + sync_patch_array_mult, sync_patch_array_mult_mp + USE mo_timer, ONLY: timer_nh_hdiffusion, timer_start, timer_stop + USE mo_exception_advanced, ONLY: finish, message, message_text + +#ifdef _OPENACC + USE mo_mpi_advanced, ONLY: i_am_accel_node +#endif + + IMPLICIT NONE + + PUBLIC :: t_diffusion, diffusion_alloc, diffusion_dealloc, diffusion_init, diffusion_run, diffusion_finalize + PRIVATE + + TYPE :: t_diffusion + LOGICAL :: lphys ! Is a run with physics? + LOGICAL :: ltimer ! Is the timer on? + LOGICAL :: l_limited_area ! Is a limited area run? + LOGICAL :: ltkeshs + LOGICAL :: lfeedback + LOGICAL :: l_zdiffu_t + LOGICAL :: lsmag_3d + LOGICAL :: lhdiff_rcf + LOGICAL :: lhdiff_w + LOGICAL :: lhdiff_temp + LOGICAL :: lvert_nest + LOGICAL :: p_test_run + LOGICAL :: ddt_vn_hdf_is_associated, ddt_vn_dyn_is_associated + INTEGER :: nproma, nlev, nlevp1, nblks_c, nblks_e, nblks_v, nrdmax, ndyn_substeps + INTEGER :: nshift, nshift_total + + INTEGER :: hdiff_order + INTEGER :: discr_vn, discr_t + INTEGER :: itype_sher + INTEGER :: itype_comm + + REAL(wp) :: grav + REAL(vp) :: cvd_o_rd + REAL(wp) :: nudge_max_coeff + REAL(wp) :: denom_diffu_v + REAL(wp) :: k4, k4w + REAL(wp) :: hdiff_smag_z, hdiff_smag_z2, hdiff_smag_z3, hdiff_smag_z4 + REAL(wp) :: hdiff_smag_fac, hdiff_smag_fac2, hdiff_smag_fac3, hdiff_smag_fac4 + REAL(wp) :: hdiff_efdt_ratio + + REAL(wp), POINTER :: vct_a(:) ! vertical coordinate part A + + REAL(wp), POINTER :: c_lin_e(:,:,:) ! p_int + REAL(wp), POINTER :: e_bln_c_s(:,:,:) ! : + REAL(wp), POINTER :: e_bln_c_u(:,:,:) ! : + REAL(wp), POINTER :: e_bln_c_v(:,:,:) ! : + REAL(wp), POINTER :: cells_aw_verts(:,:,:) + REAL(wp), POINTER :: geofac_div(:,:,:) + REAL(wp), POINTER :: geofac_rot(:,:,:) + REAL(wp), POINTER :: geofac_n2s(:,:,:) + REAL(wp), POINTER :: geofac_grg(:,:,:,:) + REAL(wp), POINTER :: nudgecoeff_e(:,:) + INTEGER, POINTER :: rbf_vec_idx_v(:,:,:) + INTEGER, POINTER :: rbf_vec_blk_v(:,:,:) + REAL(wp), POINTER :: rbf_vec_coeff_v(:,:,:,:) + + REAL(wp), POINTER :: enhfac_diffu(:) ! p_nh_metrics + REAL(wp), POINTER :: zd_intcoef(:,:) ! : + REAL(wp), POINTER :: zd_geofac(:,:) ! : + REAL(wp), POINTER :: zd_diffcoef(:) ! : + REAL(vp), POINTER :: ddqz_z_full_e(:,:,:) + REAL(vp), POINTER :: theta_ref_mc(:,:,:) + REAL(vp), POINTER :: wgtfac_c(:,:,:) + REAL(vp), POINTER :: wgtfac_e(:,:,:) + REAL(vp), POINTER :: wgtfacq_e(:,:,:) + REAL(vp), POINTER :: wgtfacq1_e(:,:,:) + INTEGER, POINTER :: zd_indlist(:,:) + INTEGER, POINTER :: zd_blklist(:,:) + INTEGER, POINTER :: zd_vertidx(:,:) + INTEGER :: zd_listdim + + END TYPE t_diffusion + + TYPE (t_diffusion), ALLOCATABLE :: diff_inst(:) + + TYPE (t_patch), ALLOCATABLE :: p_patch_diff(:) + +#ifndef __SX__ +#define __ENABLE_DDT_VN_XYZ__ +#endif + + CONTAINS + + SUBROUTINE diffusion_alloc( n_dom ) + INTEGER, INTENT(IN) :: n_dom + ALLOCATE( diff_inst( n_dom ) ) + ALLOCATE( p_patch_diff( n_dom ) ) + !$ACC ENTER DATA CREATE(diff_inst,p_patch_diff) + END SUBROUTINE diffusion_alloc + + SUBROUTINE diffusion_dealloc( ) + !$ACC EXIT DATA DELETE(diff_inst,p_patch_diff) + IF (ALLOCATED(diff_inst)) DEALLOCATE( diff_inst ) + IF (ALLOCATED(p_patch_diff)) DEALLOCATE( p_patch_diff ) + END SUBROUTINE diffusion_dealloc + + + !> + !! init_diffusion + !! + !! Prepares the horizontal diffusion of velocity and temperature + !! + !! @par Revision History + !! Initial release by William Sawyer, CSCS (2022-11-25) + !! + SUBROUTINE diffusion_init(cvd_o_rd, grav, & + jg, nproma, nlev, nblks_e, nblks_v, nblks_c, nshift, nshift_total, & + nrdmax, ndyn_substeps, nudge_max_coeff, denom_diffu_v, & + hdiff_smag_z, hdiff_smag_z2, hdiff_smag_z3, hdiff_smag_z4, & + hdiff_smag_fac, hdiff_smag_fac2, hdiff_smag_fac3, hdiff_smag_fac4, & + hdiff_order, hdiff_efdt_ratio, & + k4, k4w, itype_comm, itype_sher, itype_vn_diffu, itype_t_diffu, & + p_test_run, lphys, lhdiff_rcf, lhdiff_w, lhdiff_temp, l_limited_area,& + lfeedback, l_zdiffu_t, ltkeshs, lsmag_3d, lvert_nest, ltimer, & + ddt_vn_hdf_is_associated, ddt_vn_dyn_is_associated, & + vct_a, c_lin_e, e_bln_c_s, e_bln_c_u, e_bln_c_v, cells_aw_verts, & ! p_int + geofac_div, geofac_rot, geofac_n2s, geofac_grg, nudgecoeff_e, & ! p_int + rbf_vec_idx_v, rbf_vec_blk_v, rbf_vec_coeff_v, & ! p_int + enhfac_diffu, zd_intcoef, zd_geofac, zd_diffcoef, & ! p_nh_metrics + wgtfac_c, wgtfac_e, wgtfacq_e, wgtfacq1_e, & ! p_nh_metrics + ddqz_z_full_e, theta_ref_mc, & ! p_nh_metrics + zd_indlist, zd_blklist, zd_vertidx, zd_listdim, & ! p_nh_metrics + edges_start_block, edges_end_block, edges_start_index, edges_end_index,&! p_patch%edges + edges_vertex_idx, edges_vertex_blk, edges_cell_idx, edges_cell_blk, & ! p_patch%edges + edges_tangent_orientation, & ! p_patch%edges + edges_primal_normal_vert, edges_dual_normal_vert, & ! p_patch%edges + edges_primal_normal_cell, edges_dual_normal_cell, & ! p_patch%edges + edges_inv_vert_vert_length, edges_inv_primal_edge_length, & ! p_patch%edges + edges_inv_dual_edge_length, edges_area_edge, & ! p_patch%edges + cells_start_block, cells_end_block, cells_start_index, cells_end_index,&! p_patch%cells + cells_neighbor_idx, cells_neighbor_blk, & ! p_patch%cells + cells_edge_idx, cells_edge_blk, cells_area, & ! p_patch%cells + verts_start_block, verts_end_block, verts_start_index, verts_end_index )! p_patch%verts + REAL(wp), INTENT(IN) :: cvd_o_rd, grav ! Physical constants from central location + INTEGER, INTENT(IN) :: jg + INTEGER, INTENT(IN) :: nproma, nlev, nblks_e, nblks_v, nblks_c, nshift, nshift_total + INTEGER, INTENT(IN) :: nrdmax ! = nrdmax(jg) + INTEGER, INTENT(IN) :: ndyn_substeps + INTEGER, INTENT(IN) :: hdiff_order, itype_comm, itype_sher, itype_vn_diffu, itype_t_diffu + REAL(wp), INTENT(IN) :: hdiff_smag_z, hdiff_smag_z2, hdiff_smag_z3, hdiff_smag_z4 + REAL(wp), INTENT(IN) :: hdiff_smag_fac, hdiff_smag_fac2, hdiff_smag_fac3, hdiff_smag_fac4 + REAL(wp), INTENT(IN) :: hdiff_efdt_ratio + REAL(wp), INTENT(IN) :: k4, k4w + REAL(wp), INTENT(IN) :: nudge_max_coeff + REAL(wp), INTENT(IN) :: denom_diffu_v + LOGICAL, INTENT(IN) :: p_test_run + LOGICAL, INTENT(IN) :: lphys !< is a run with physics + LOGICAL, INTENT(IN) :: lhdiff_rcf + LOGICAL, INTENT(IN) :: lhdiff_w + LOGICAL, INTENT(IN) :: lhdiff_temp + LOGICAL, INTENT(IN) :: l_zdiffu_t + LOGICAL, INTENT(IN) :: l_limited_area + LOGICAL, INTENT(IN) :: lfeedback ! = lfeedback(jg) + LOGICAL, INTENT(IN) :: ltkeshs + LOGICAL, INTENT(IN) :: lsmag_3d + LOGICAL, INTENT(IN) :: lvert_nest + LOGICAL, INTENT(IN) :: ltimer + LOGICAL, INTENT(IN) :: ddt_vn_hdf_is_associated + LOGICAL, INTENT(IN) :: ddt_vn_dyn_is_associated + + REAL(wp), TARGET, INTENT(IN) :: vct_a(:) ! param. A of the vertical coordinte + + REAL(wp), TARGET, INTENT(IN) :: c_lin_e(:,:,:) ! p_int + REAL(wp), TARGET, INTENT(IN) :: e_bln_c_s(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: e_bln_c_u(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: e_bln_c_v(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: cells_aw_verts(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: geofac_div(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: geofac_rot(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: geofac_n2s(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: geofac_grg(:,:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: nudgecoeff_e(:,:) ! : + INTEGER, TARGET, INTENT(IN) :: rbf_vec_idx_v(:,:,:) + INTEGER, TARGET, INTENT(IN) :: rbf_vec_blk_v(:,:,:) + REAL(wp), TARGET, INTENT(IN) :: rbf_vec_coeff_v(:,:,:,:) + + REAL(wp), TARGET, INTENT(IN) :: enhfac_diffu(:) ! p_nh_metrics + REAL(wp), TARGET, INTENT(IN) :: zd_intcoef(:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: zd_geofac(:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: zd_diffcoef(:) ! : + REAL(vp), TARGET, INTENT(IN) :: wgtfac_c(:,:,:) ! : + REAL(vp), TARGET, INTENT(IN) :: wgtfac_e(:,:,:) ! : + REAL(vp), TARGET, INTENT(IN) :: wgtfacq_e(:,:,:) ! : + REAL(vp), TARGET, INTENT(IN) :: wgtfacq1_e(:,:,:) ! : + REAL(vp), TARGET, INTENT(IN) :: ddqz_z_full_e(:,:,:) ! : + REAL(vp), TARGET, INTENT(IN) :: theta_ref_mc(:,:,:) ! : + INTEGER, TARGET, INTENT(IN) :: zd_indlist(:,:) ! : + INTEGER, TARGET, INTENT(IN) :: zd_blklist(:,:) ! : + INTEGER, TARGET, INTENT(IN) :: zd_vertidx(:,:) ! : + INTEGER, INTENT(IN) :: zd_listdim ! : + + INTEGER, TARGET, INTENT(IN) :: edges_start_block(min_rledge:) ! p_patch%edges + INTEGER, TARGET, INTENT(IN) :: edges_end_block(min_rledge:) ! : + INTEGER, TARGET, INTENT(IN) :: edges_start_index(min_rledge:) ! p_patch%edges + INTEGER, TARGET, INTENT(IN) :: edges_end_index(min_rledge:) ! : + INTEGER, TARGET, INTENT(IN) :: edges_vertex_idx(:,:,:) ! : + INTEGER, TARGET, INTENT(IN) :: edges_vertex_blk(:,:,:) ! : + INTEGER, TARGET, INTENT(IN) :: edges_cell_idx(:,:,:) ! : + INTEGER, TARGET, INTENT(IN) :: edges_cell_blk(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: edges_tangent_orientation(:,:) + TYPE(t_tangent_vectors), TARGET, INTENT(IN) :: edges_primal_normal_vert(:,:,:) + TYPE(t_tangent_vectors), TARGET, INTENT(IN) :: edges_dual_normal_vert(:,:,:) + TYPE(t_tangent_vectors), TARGET, INTENT(IN) :: edges_primal_normal_cell(:,:,:) + TYPE(t_tangent_vectors), TARGET, INTENT(IN) :: edges_dual_normal_cell(:,:,:) + REAL(wp), TARGET, INTENT(IN) :: edges_inv_vert_vert_length(:,:) + REAL(wp), TARGET, INTENT(IN) :: edges_inv_primal_edge_length(:,:) + REAL(wp), TARGET, INTENT(IN) :: edges_inv_dual_edge_length(:,:) + REAL(wp), TARGET, INTENT(IN) :: edges_area_edge(:,:) + + INTEGER, TARGET, INTENT(IN) :: cells_start_block(min_rlcell:) ! p_patch%cells + INTEGER, TARGET, INTENT(IN) :: cells_end_block(min_rlcell:) ! : + INTEGER, TARGET, INTENT(IN) :: cells_start_index(min_rlcell:) ! p_patch%cells + INTEGER, TARGET, INTENT(IN) :: cells_end_index(min_rlcell:) ! : + INTEGER, TARGET, INTENT(IN) :: cells_neighbor_idx(:,:,:) ! : + INTEGER, TARGET, INTENT(IN) :: cells_neighbor_blk(:,:,:) ! : + INTEGER, TARGET, INTENT(IN) :: cells_edge_idx(:,:,:) ! : + INTEGER, TARGET, INTENT(IN) :: cells_edge_blk(:,:,:) ! : + REAL(wp), TARGET, INTENT(IN) :: cells_area(:,:) + + INTEGER, TARGET, INTENT(IN) :: verts_start_block(min_rlvert:) ! p_patch%verts + INTEGER, TARGET, INTENT(IN) :: verts_end_block(min_rlvert:) ! : + INTEGER, TARGET, INTENT(IN) :: verts_start_index(min_rlvert:) ! p_patch%verts + INTEGER, TARGET, INTENT(IN) :: verts_end_index(min_rlvert:) ! : + !-------------------------------------------------------------------------- + + diff_inst(jg)%vct_a => vct_a + + diff_inst(jg)%c_lin_e => c_lin_e ! p_int + diff_inst(jg)%e_bln_c_s => e_bln_c_s ! : + diff_inst(jg)%e_bln_c_u => e_bln_c_u ! : + diff_inst(jg)%e_bln_c_v => e_bln_c_v ! : + diff_inst(jg)%cells_aw_verts => cells_aw_verts + diff_inst(jg)%geofac_div => geofac_div + diff_inst(jg)%geofac_rot => geofac_rot + diff_inst(jg)%geofac_n2s => geofac_n2s + diff_inst(jg)%geofac_grg => geofac_grg + diff_inst(jg)%nudgecoeff_e => nudgecoeff_e + diff_inst(jg)%rbf_vec_idx_v => rbf_vec_idx_v + diff_inst(jg)%rbf_vec_blk_v => rbf_vec_blk_v + diff_inst(jg)%rbf_vec_coeff_v => rbf_vec_coeff_v + + diff_inst(jg)%enhfac_diffu => enhfac_diffu ! p_nh_metrics + diff_inst(jg)%zd_intcoef => zd_intcoef ! : + diff_inst(jg)%zd_geofac => zd_geofac + diff_inst(jg)%zd_diffcoef => zd_diffcoef + diff_inst(jg)%wgtfac_c => wgtfac_c + diff_inst(jg)%wgtfac_e => wgtfac_e + diff_inst(jg)%wgtfacq_e => wgtfacq_e + diff_inst(jg)%wgtfacq1_e => wgtfacq1_e + diff_inst(jg)%ddqz_z_full_e => ddqz_z_full_e + diff_inst(jg)%theta_ref_mc => theta_ref_mc + diff_inst(jg)%zd_indlist => zd_indlist + diff_inst(jg)%zd_blklist => zd_blklist + diff_inst(jg)%zd_vertidx => zd_vertidx + diff_inst(jg)%zd_listdim = zd_listdim + + p_patch_diff(jg)%edges%start_block => edges_start_block ! p_patch%edges + p_patch_diff(jg)%edges%end_block => edges_end_block ! : + p_patch_diff(jg)%edges%start_index => edges_start_index ! p_patch%edges + p_patch_diff(jg)%edges%end_index => edges_end_index ! : + p_patch_diff(jg)%edges%vertex_idx => edges_vertex_idx + p_patch_diff(jg)%edges%vertex_blk => edges_vertex_blk + p_patch_diff(jg)%edges%cell_idx => edges_cell_idx + p_patch_diff(jg)%edges%cell_blk => edges_cell_blk + p_patch_diff(jg)%edges%tangent_orientation => edges_tangent_orientation + p_patch_diff(jg)%edges%primal_normal_vert => edges_primal_normal_vert + p_patch_diff(jg)%edges%dual_normal_vert => edges_dual_normal_vert + p_patch_diff(jg)%edges%primal_normal_cell => edges_primal_normal_cell + p_patch_diff(jg)%edges%dual_normal_cell => edges_dual_normal_cell + p_patch_diff(jg)%edges%inv_vert_vert_length => edges_inv_vert_vert_length + p_patch_diff(jg)%edges%inv_primal_edge_length => edges_inv_primal_edge_length + p_patch_diff(jg)%edges%inv_dual_edge_length => edges_inv_dual_edge_length + p_patch_diff(jg)%edges%area_edge => edges_area_edge + + p_patch_diff(jg)%cells%start_block => cells_start_block ! p_patch%cells + p_patch_diff(jg)%cells%end_block => cells_end_block ! : + p_patch_diff(jg)%cells%start_index => cells_start_index ! p_patch%cells + p_patch_diff(jg)%cells%end_index => cells_end_index ! : + p_patch_diff(jg)%cells%neighbor_idx => cells_neighbor_idx + p_patch_diff(jg)%cells%neighbor_blk => cells_neighbor_blk + p_patch_diff(jg)%cells%edge_idx => cells_edge_idx + p_patch_diff(jg)%cells%edge_blk => cells_edge_blk + p_patch_diff(jg)%cells%area => cells_area + + p_patch_diff(jg)%verts%start_block => verts_start_block ! p_patch%cells + p_patch_diff(jg)%verts%end_block => verts_end_block ! : + p_patch_diff(jg)%verts%start_index => verts_start_index ! p_patch%cells + p_patch_diff(jg)%verts%end_index => verts_end_index ! : + + diff_inst(jg)%nrdmax = nrdmax + diff_inst(jg)%grav = grav ! from central location + diff_inst(jg)%cvd_o_rd = cvd_o_rd ! " + diff_inst(jg)%nudge_max_coeff= nudge_max_coeff + diff_inst(jg)%denom_diffu_v = denom_diffu_v + diff_inst(jg)%k4 = k4 + diff_inst(jg)%k4w = k4w + + ! number of vertical levels, blocks for edges, vertices and cells + diff_inst(jg)%nproma = nproma + diff_inst(jg)%nlev = nlev + diff_inst(jg)%nlevp1 = nlev+1 + diff_inst(jg)%nblks_e = nblks_e + diff_inst(jg)%nblks_v = nblks_v + diff_inst(jg)%nblks_c = nblks_c + + diff_inst(jg)%ndyn_substeps = ndyn_substeps + diff_inst(jg)%nshift = nshift ! p_patch%nshift + diff_inst(jg)%nshift_total = nshift_total ! p_patch%nshift_total + + diff_inst(jg)%itype_sher = itype_sher + diff_inst(jg)%itype_comm = itype_comm + + diff_inst(jg)%p_test_run = p_test_run + diff_inst(jg)%lphys = lphys + diff_inst(jg)%l_zdiffu_t = l_zdiffu_t + diff_inst(jg)%l_limited_area = l_limited_area + diff_inst(jg)%lfeedback = lfeedback + diff_inst(jg)%ltkeshs = ltkeshs + diff_inst(jg)%lsmag_3d = lsmag_3d + diff_inst(jg)%lhdiff_w = lhdiff_w + diff_inst(jg)%lhdiff_rcf = lhdiff_rcf + diff_inst(jg)%lhdiff_temp = lhdiff_temp + diff_inst(jg)%ltimer = ltimer + diff_inst(jg)%lvert_nest = lvert_nest + diff_inst(jg)%ddt_vn_hdf_is_associated = ddt_vn_hdf_is_associated + diff_inst(jg)%ddt_vn_dyn_is_associated = ddt_vn_dyn_is_associated + + diff_inst(jg)%hdiff_order = hdiff_order + + diff_inst(jg)%discr_vn = itype_vn_diffu + diff_inst(jg)%discr_t = itype_t_diffu + + diff_inst(jg)%hdiff_smag_z = hdiff_smag_z + diff_inst(jg)%hdiff_smag_z2 = hdiff_smag_z2 + diff_inst(jg)%hdiff_smag_z3 = hdiff_smag_z3 + diff_inst(jg)%hdiff_smag_z4 = hdiff_smag_z4 + diff_inst(jg)%hdiff_smag_fac = hdiff_smag_fac + diff_inst(jg)%hdiff_smag_fac2= hdiff_smag_fac2 + diff_inst(jg)%hdiff_smag_fac3= hdiff_smag_fac3 + diff_inst(jg)%hdiff_smag_fac4= hdiff_smag_fac4 + diff_inst(jg)%hdiff_efdt_ratio = hdiff_efdt_ratio + + !$ACC ENTER DATA COPYIN(diff_inst(jg)) + + END SUBROUTINE diffusion_init + + !> + !! diffusion + !! + !! Computes the horizontal diffusion of velocity and temperature + !! + !! @par Revision History + !! Initial release by Guenther Zaengl, DWD (2010-10-13), based on an earlier + !! version initially developed by Almut Gassmann, MPI-M + !! + SUBROUTINE diffusion_run(jg, dtime, linit, & + vn, w, theta_v, exner, & ! p_nh_prog + vt, theta_v_ic, div_ic, hdef_ic, dwdx, dwdy, & ! p_nh_diag + ddt_vn_dyn, ddt_vn_hdf ) ! p_nh_diag optional + + INTEGER, INTENT(IN) :: jg ! patch ID + REAL(wp), INTENT(IN) :: dtime !< time step + LOGICAL, INTENT(IN) :: linit !< initial call or runtime call + REAL(wp), INTENT(INOUT) :: vn(:,:,:) ! orthogonal normal wind (nproma,nlev,nblks_e) [m/s] + REAL(wp), INTENT(INOUT) :: w(:,:,:) ! orthogonal vertical wind (nproma,nlevp1,nblks_c) [m/s] + REAL(wp), INTENT(INOUT) :: theta_v(:,:,:) ! virtual potential temperature (nproma,nlev,nblks_c) [K] + REAL(wp), INTENT(INOUT) :: exner(:,:,:) ! Exner pressure (nproma,nlev,nblks_c) [-] + REAL(wp), INTENT(INOUT) :: theta_v_ic(:,:,:) ! theta_v at half levels (nproma,nlevp1,nblks_c) [K] + REAL(vp), INTENT(IN) :: vt(:,:,:) ! tangential wind (nproma,nlev,nblks_e) [m/s] + REAL(vp), INTENT(OUT) :: div_ic(:,:,:) ! divergence at half levels(nproma,nlevp1,nblks_c) [1/s] + REAL(vp), INTENT(OUT) :: hdef_ic(:,:,:) ! horizontal wind field deformation (nproma,nlevp1,nblks_c) [1/s^2] + REAL(vp), INTENT(OUT) :: dwdx(:,:,:) ! divergence at half levels(nproma,nlevp1,nblks_c) [1/s] + REAL(vp), INTENT(OUT) :: dwdy(:,:,:) ! horizontal wind field deformation (nproma,nlevp1,nblks_c) [1/s^2] + REAL(wp), INTENT(INOUT), OPTIONAL :: ddt_vn_dyn(:,:,:) ! d vn / dt (sum of all contributions) + REAL(wp), INTENT(INOUT), OPTIONAL :: ddt_vn_hdf(:,:,:) ! d vn / dt (horizontal diffusion only) + + ! local variables - vp means variable precision depending on the __MIXED_PRECISION cpp flag + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_c) :: z_temp + REAL(wp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_e) :: z_nabla2_e + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_c) :: z_nabla2_c + REAL(wp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_e) :: z_nabla4_e + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev) :: z_nabla4_e2 + + REAL(wp):: diff_multfac_vn(diff_inst(jg)%nlev), diff_multfac_w, diff_multfac_n2w(diff_inst(jg)%nlev) + INTEGER :: i_startblk, i_endblk, i_startidx, i_endidx + INTEGER :: rl_start, rl_end + INTEGER :: jk, jb, jc, je, ic, ishift, nshift, jk1 + INTEGER :: nlev, nlevp1 !< number of full and half levels + + ! start index levels and diffusion coefficient for boundary diffusion + INTEGER, PARAMETER :: start_bdydiff_e = 5 ! refin_ctrl level at which boundary diffusion starts + REAL(wp):: fac_bdydiff_v + + ! For Smagorinsky diffusion - vp means variable precision depending on the __MIXED_PRECISION cpp flag + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_e) :: kh_smag_e + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_e) :: kh_smag_ec + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_v) :: u_vert + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_v) :: v_vert + REAL(wp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_c) :: u_cell + REAL(wp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev,diff_inst(jg)%nblks_c) :: v_cell + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev) :: kh_c, div + + REAL(vp) :: dvt_norm, dvt_tang, vn_vert1, vn_vert2, vn_vert3, vn_vert4, vn_cell1, vn_cell2 + + REAL(vp) :: smag_offset, nabv_tang, nabv_norm, rd_o_cvd, nudgezone_diff, bdy_diff, enh_diffu + REAL(vp), DIMENSION(diff_inst(jg)%nlev) :: smag_limit, diff_multfac_smag, enh_smag_fac + INTEGER :: nblks_zdiffu, nproma_zdiffu, npromz_zdiffu, nlen_zdiffu + + REAL(wp) :: alin, dz32, df32, dz42, df42, bqdr, aqdr, zf, dzlin, dzqdr + + ! Additional variables for 3D Smagorinsky coefficient + REAL(wp):: z_w_v(diff_inst(jg)%nproma,diff_inst(jg)%nlevp1,diff_inst(jg)%nblks_v) + REAL(wp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlevp1) :: z_vn_ie, z_vt_ie + REAL(wp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev) :: dvndz, dvtdz, dwdz, dthvdz, dwdn, dwdt, kh_smag3d_e + + ! Variables for provisional fix against runaway cooling in local topography depressions + INTEGER :: icount(diff_inst(jg)%nblks_c), iclist(2*diff_inst(jg)%nproma,diff_inst(jg)%nblks_c), iklist(2*diff_inst(jg)%nproma,diff_inst(jg)%nblks_c) + REAL(wp) :: tdlist(2*diff_inst(jg)%nproma,diff_inst(jg)%nblks_c), tdiff, trefdiff, thresh_tdiff, z_theta, fac2d + + + INTEGER, DIMENSION(:,:,:), POINTER :: icidx, icblk, ieidx, ieblk, ividx, ivblk, iecidx, iecblk + INTEGER, DIMENSION(:,:), POINTER :: icell, ilev, iblk !, iedge, iedblk + REAL(wp), DIMENSION(:,:), POINTER :: vcoef, zd_geofac !, blcoef + LOGICAL :: ltemp_diffu + INTEGER :: diffu_type + +#ifdef _OPENACC + REAL(vp), DIMENSION(diff_inst(jg)%nproma,diff_inst(jg)%nlev-1:diff_inst(jg)%nlev,diff_inst(jg)%nblks_c) :: enh_diffu_3d +#endif + + ! Variables for tendency diagnostics + REAL(wp) :: z_d_vn_hdf + REAL(wp) :: r_dtimensubsteps + + CHARACTER(*), PARAMETER :: routine = "diffusion_run" + + !-------------------------------------------------------------------------- + + ividx => p_patch_diff(jg)%edges%vertex_idx + ivblk => p_patch_diff(jg)%edges%vertex_blk + + iecidx => p_patch_diff(jg)%edges%cell_idx + iecblk => p_patch_diff(jg)%edges%cell_blk + + icidx => p_patch_diff(jg)%cells%neighbor_idx + icblk => p_patch_diff(jg)%cells%neighbor_blk + + ieidx => p_patch_diff(jg)%cells%edge_idx + ieblk => p_patch_diff(jg)%cells%edge_blk + + ! prepare for tendency diagnostics + IF (diff_inst(jg)%lhdiff_rcf) THEN + r_dtimensubsteps = 1._wp/dtime ! without substepping, no averaging is necessary + ELSE + r_dtimensubsteps = 1._wp/(dtime*REAL(diff_inst(jg)%ndyn_substeps,wp)) ! with substepping the tendency is averaged over the substeps + END IF + + ! number of vertical levels + nlev = diff_inst(jg)%nlev + nlevp1 = nlev+1 + + ! Normalized diffusion coefficient for boundary diffusion + IF (diff_inst(jg)%lhdiff_rcf) THEN + fac_bdydiff_v = SQRT(REAL(diff_inst(jg)%ndyn_substeps,wp))/diff_inst(jg)%denom_diffu_v + ELSE + fac_bdydiff_v = 1._wp/diff_inst(jg)%denom_diffu_v + ENDIF + + ! scaling factor for enhanced diffusion in nudging zone (if present, i.e. for + ! limited-area runs and one-way nesting) + nudgezone_diff = 0.04_wp/(diff_inst(jg)%nudge_max_coeff + dbl_eps) + + ! scaling factor for enhanced near-boundary diffusion for + ! two-way nesting (used with Smagorinsky diffusion only; not needed otherwise) + bdy_diff = 0.015_wp/(diff_inst(jg)%nudge_max_coeff + dbl_eps) + + ! threshold temperature deviation from neighboring grid points + ! that activates extra diffusion against runaway cooling + thresh_tdiff = - 5._wp + + rd_o_cvd = 1._wp/diff_inst(jg)%cvd_o_rd + diffu_type = diff_inst(jg)%hdiff_order + + + IF (linit) THEN ! enhanced diffusion at all levels for initial velocity filtering call + diff_multfac_vn(:) = diff_inst(jg)%k4/3._wp*diff_inst(jg)%hdiff_efdt_ratio + smag_offset = 0.0_vp + diffu_type = 5 ! always combine nabla4 background diffusion with Smagorinsky diffusion for initial filtering call + smag_limit(:) = 0.125_wp-4._wp*diff_multfac_vn(:) + ELSE IF (diff_inst(jg)%lhdiff_rcf) THEN ! combination with divergence damping inside the dynamical core + IF (diffu_type == 4) THEN + diff_multfac_vn(:) = MIN(1._wp/128._wp,diff_inst(jg)%k4*REAL(diff_inst(jg)%ndyn_substeps,wp)/ & + 3._wp*diff_inst(jg)%enhfac_diffu(:)) + ELSE ! For Smagorinsky diffusion, the Smagorinsky coefficient rather than the background + ! diffusion coefficient is enhanced near the model top (see below) + diff_multfac_vn(:) = MIN(1._wp/128._wp,diff_inst(jg)%k4*REAL(diff_inst(jg)%ndyn_substeps,wp)/3._wp) + ENDIF + IF (diffu_type == 3) THEN + smag_offset = 0._vp + smag_limit(:) = 0.125_vp + ELSE + smag_offset = 0.25_wp*diff_inst(jg)%k4*REAL(diff_inst(jg)%ndyn_substeps,wp) + smag_limit(:) = 0.125_wp-4._wp*diff_multfac_vn(:) + ENDIF + ELSE ! enhanced diffusion near model top only + IF (diffu_type == 4) THEN + diff_multfac_vn(:) = diff_inst(jg)%k4/3._wp*diff_inst(jg)%enhfac_diffu(:) + ELSE ! For Smagorinsky diffusion, the Smagorinsky coefficient rather than the background + ! diffusion coefficient is enhanced near the model top (see below) + diff_multfac_vn(:) = diff_inst(jg)%k4/3._wp + ENDIF + smag_offset = 0.25_wp*diff_inst(jg)%k4 + smag_limit(:) = 0.125_wp-4._wp*diff_multfac_vn(:) + ! pure Smagorinsky diffusion does not work without divergence damping + IF (diff_inst(jg)%hdiff_order == 3) diffu_type = 5 + ENDIF + + ! Multiplication factor for nabla4 diffusion on vertical wind speed + diff_multfac_w = MIN(1._wp/48._wp,diff_inst(jg)%k4w*REAL(diff_inst(jg)%ndyn_substeps,wp)) + + ! Factor for additional nabla2 diffusion in upper damping zone + diff_multfac_n2w(:) = 0._wp + IF (diff_inst(jg)%nrdmax > 1) THEN ! seems to be redundant, but the NEC issues invalid operations otherwise + DO jk = 2, diff_inst(jg)%nrdmax + jk1 = jk + diff_inst(jg)%nshift_total + diff_multfac_n2w(jk) = 1._wp/12._wp*((diff_inst(jg)%vct_a(jk1)-diff_inst(jg)%vct_a(diff_inst(jg)%nshift_total+diff_inst(jg)%nrdmax+1))/ & + (diff_inst(jg)%vct_a(2)-diff_inst(jg)%vct_a(diff_inst(jg)%nshift_total+diff_inst(jg)%nrdmax+1)))**4 + ENDDO + ENDIF + + IF (diffu_type == 3 .OR. diffu_type == 5) THEN + + ! temperature diffusion is used only in combination with Smagorinsky diffusion + ltemp_diffu = diff_inst(jg)%lhdiff_temp + + ! The Smagorinsky diffusion factor enh_divdamp_fac is defined as a profile in height z + ! above sea level with 4 height sections: + ! + ! enh_smag_fac(z) = hdiff_smag_fac ! z <= hdiff_smag_z + ! enh_smag_fac(z) = hdiff_smag_fac + (z-hdiff_smag_z )* alin ! hdiff_smag_z <= z <= hdiff_smag_z2 + ! enh_smag_fac(z) = hdiff_smag_fac2 + (z-hdiff_smag_z2)*(aqdr+(z-hdiff_smag_z2)*bqdr) ! hdiff_smag_z2 <= z <= hdiff_smag_z4 + ! enh_smag_fac(z) = hdiff_smag_fac4 ! hdiff_smag_z4 <= z + ! + alin = (diff_inst(jg)%hdiff_smag_fac2-diff_inst(jg)%hdiff_smag_fac)/ & + & (diff_inst(jg)%hdiff_smag_z2 -diff_inst(jg)%hdiff_smag_z) + ! + df32 = diff_inst(jg)%hdiff_smag_fac3-diff_inst(jg)%hdiff_smag_fac2 + df42 = diff_inst(jg)%hdiff_smag_fac4-diff_inst(jg)%hdiff_smag_fac2 + ! + dz32 = diff_inst(jg)%hdiff_smag_z3-diff_inst(jg)%hdiff_smag_z2 + dz42 = diff_inst(jg)%hdiff_smag_z4-diff_inst(jg)%hdiff_smag_z2 + ! + bqdr = (df42*dz32-df32*dz42)/(dz32*dz42*(dz42-dz32)) + aqdr = df32/dz32-bqdr*dz32 + ! + DO jk = 1, nlev + jk1 = jk + diff_inst(jg)%nshift_total + ! + zf = 0.5_wp*(diff_inst(jg)%vct_a(jk1)+diff_inst(jg)%vct_a(jk1+1)) + dzlin = MIN( diff_inst(jg)%hdiff_smag_z2-diff_inst(jg)%hdiff_smag_z , & + & MAX( 0._wp, zf-diff_inst(jg)%hdiff_smag_z ) ) + dzqdr = MIN( diff_inst(jg)%hdiff_smag_z4-diff_inst(jg)%hdiff_smag_z2, & + & MAX( 0._wp, zf-diff_inst(jg)%hdiff_smag_z2) ) + ! + enh_smag_fac(jk) = REAL(diff_inst(jg)%hdiff_smag_fac + dzlin*alin + dzqdr*(aqdr+dzqdr*bqdr),vp) + ! + ENDDO + + ! Smagorinsky coefficient is also enhanced in the six model levels beneath a vertical nest interface + IF (diff_inst(jg)%lvert_nest .AND. (diff_inst(jg)%nshift > 0)) THEN + enh_smag_fac(1) = MAX(0.333_vp, enh_smag_fac(1)) + enh_smag_fac(2) = MAX(0.25_vp, enh_smag_fac(2)) + enh_smag_fac(3) = MAX(0.20_vp, enh_smag_fac(3)) + enh_smag_fac(4) = MAX(0.16_vp, enh_smag_fac(4)) + enh_smag_fac(5) = MAX(0.12_vp, enh_smag_fac(5)) + enh_smag_fac(6) = MAX(0.08_vp, enh_smag_fac(6)) + ENDIF + + ! empirically determined scaling factor + diff_multfac_smag(:) = enh_smag_fac(:)*REAL(dtime,vp) + + ELSE + ltemp_diffu = .FALSE. + ENDIF + + !$ACC DATA CREATE(div, kh_c, kh_smag_e, kh_smag_ec, u_vert, v_vert, u_cell, v_cell, z_w_v, z_temp) & + !$ACC CREATE(z_nabla4_e, z_nabla4_e2, z_nabla2_e, z_nabla2_c, enh_diffu_3d, icount) & + !$ACC CREATE(z_vn_ie, z_vt_ie, dvndz, dvtdz, dwdz, dthvdz, dwdn, dwdt, kh_smag3d_e) & + !$ACC COPYIN(diff_multfac_vn, diff_multfac_n2w, diff_multfac_smag, smag_limit) & + !$ACC PRESENT(diff_inst, p_patch_diff) & + !$ACC PRESENT(ividx, ivblk, iecidx, iecblk, icidx, icblk, ieidx, ieblk) & + !$ACC IF(i_am_accel_node) + + !!! Following variables may be present in certain situations, but we don't want it to fail in the general case. + !!! Should actually be in a separate data region with correct IF condition. + !!! !$ACC div_ic, dwdx, dwdy, hdef_ic, & + + ! The diffusion is an intrinsic part of the NH solver, thus it is added to the timer + IF (diff_inst(jg)%ltimer) CALL timer_start(timer_nh_hdiffusion) + + IF (diffu_type == 4) THEN + + CALL nabla4_vec( vn, jg, p_patch_diff(jg), diff_inst(jg)%geofac_div, diff_inst(jg)%geofac_rot, & + z_nabla4_e, opt_rlstart=7,opt_nabla2=z_nabla2_e ) + ELSE IF ((diffu_type == 3 .OR. diffu_type == 5) & + .AND. diff_inst(jg)%discr_vn == 1 .AND. .NOT. diff_inst(jg)%lsmag_3d) THEN + + IF (diff_inst(jg)%p_test_run) THEN + !$ACC KERNELS PRESENT(u_vert, v_vert) ASYNC(1) IF(i_am_accel_node) + u_vert = 0._vp + v_vert = 0._vp + !$ACC END KERNELS + ENDIF + + ! RBF reconstruction of velocity at vertices + CALL rbf_vec_interpol_vertex( vn, p_patch_diff(jg), & + diff_inst(jg)%rbf_vec_idx_v, diff_inst(jg)%rbf_vec_blk_v, & + diff_inst(jg)%rbf_vec_coeff_v, u_vert, v_vert, & + opt_rlend=min_rlvert_int, opt_acc_async=.TRUE. ) + rl_start = start_bdydiff_e + rl_end = min_rledge_int - 2 + + IF (diff_inst(jg)%itype_comm == 1 .OR. diff_inst(jg)%itype_comm == 3) THEN +#ifdef __MIXED_PRECISION + CALL sync_patch_array_mult_mp(SYNC_V,p_patch(jg),0,2,f3din1_sp=u_vert,f3din2_sp=v_vert, & + opt_varname="diffusion: u_vert and v_vert") +#else + CALL sync_patch_array_mult(SYNC_V,p_patch(jg),2,u_vert,v_vert, & + opt_varname="diffusion: u_vert and v_vert") +#endif + ENDIF + +!$OMP PARALLEL PRIVATE(i_startblk,i_endblk) + + i_startblk = p_patch_diff(jg)%edges%start_block(rl_start) + i_endblk = p_patch_diff(jg)%edges%end_block(rl_end) + +!$OMP DO PRIVATE(jb,i_startidx,i_endidx,jk,je,vn_vert1,vn_vert2,vn_vert3,vn_vert4, & +!$OMP dvt_norm,dvt_tang), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + ! Computation of wind field deformation + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO je = i_startidx, i_endidx +!DIR$ IVDEP + DO jk = 1, nlev +#else +!$NEC outerloop_unroll(4) + DO jk = 1, nlev + DO je = i_startidx, i_endidx +#endif + + vn_vert1 = u_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v1 + & + v_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v2 + + vn_vert2 = u_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v1 + & + v_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v2 + + dvt_tang = p_patch_diff(jg)%edges%tangent_orientation(je,jb)* ( & + u_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,2)%v1 + & + v_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,2)%v2 - & + (u_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,1)%v1 + & + v_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,1)%v2) ) + + vn_vert3 = u_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,3)%v1 + & + v_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,3)%v2 + + vn_vert4 = u_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,4)%v1 + & + v_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,4)%v2 + + dvt_norm = u_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,4)%v1 + & + v_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,4)%v2 - & + (u_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,3)%v1 + & + v_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,3)%v2) + ! Smagorinsky diffusion coefficient + kh_smag_e(je,jk,jb) = diff_multfac_smag(jk)*SQRT( ( & + (vn_vert4-vn_vert3)*p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb)- & + dvt_tang*p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) )**2 + ( & + (vn_vert2-vn_vert1)*p_patch_diff(jg)%edges%tangent_orientation(je,jb)* & + p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) + & + dvt_norm*p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb))**2 ) + + ! The factor of 4 comes from dividing by twice the "correct" length + z_nabla2_e(je,jk,jb) = 4._wp * ( & + (vn_vert4 + vn_vert3 - 2._wp*vn(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb)**2 + & + (vn_vert2 + vn_vert1 - 2._wp*vn(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb)**2 ) + +#if defined (__LOOP_EXCHANGE) && !defined (_OPENACC) + ENDDO + ENDDO + + DO jk = 1, nlev + DO je = i_startidx, i_endidx +#endif + kh_smag_ec(je,jk,jb) = kh_smag_e(je,jk,jb) + ! Subtract part of the fourth-order background diffusion coefficient + kh_smag_e(je,jk,jb) = MAX(0._vp,kh_smag_e(je,jk,jb) - smag_offset) + ! Limit diffusion coefficient to the theoretical CFL stability threshold + kh_smag_e(je,jk,jb) = MIN(kh_smag_e(je,jk,jb),smag_limit(jk)) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO ! block jb +!$OMP END DO NOWAIT +!$OMP END PARALLEL + + ELSE IF ((diffu_type == 3 .OR. diffu_type == 5) .AND. diff_inst(jg)%discr_vn == 1) THEN + ! 3D Smagorinsky diffusion + IF (diff_inst(jg)%p_test_run) THEN + !$ACC KERNELS PRESENT(u_vert, v_vert, z_w_v) ASYNC(1) IF(i_am_accel_node) + u_vert = 0._vp + v_vert = 0._vp + z_w_v = 0._wp + !$ACC END KERNELS + ENDIF + + ! RBF reconstruction of velocity at vertices + CALL rbf_vec_interpol_vertex( vn, p_patch_diff(jg), & + diff_inst(jg)%rbf_vec_idx_v, diff_inst(jg)%rbf_vec_blk_v, & + diff_inst(jg)%rbf_vec_coeff_v, u_vert, v_vert, & + opt_rlend=min_rlvert_int, opt_acc_async=.TRUE. ) + + rl_start = start_bdydiff_e + rl_end = min_rledge_int - 2 + + IF (diff_inst(jg)%itype_comm == 1 .OR. diff_inst(jg)%itype_comm == 3) THEN +#ifdef __MIXED_PRECISION + CALL sync_patch_array_mult_mp(SYNC_V,p_patch(jg),0,2,f3din1_sp=u_vert,f3din2_sp=v_vert, & + opt_varname="diffusion: u_vert and v_vert 2") +#else + CALL sync_patch_array_mult(SYNC_V,p_patch(jg),2,u_vert,v_vert, & + opt_varname="diffusion: u_vert and v_vert 2") +#endif + ENDIF + CALL cells2verts_scalar(w, p_patch_diff(jg), diff_inst(jg)%cells_aw_verts, z_w_v, opt_rlend=min_rlvert_int) + CALL sync_patch_array(SYNC_V,p_patch(jg),z_w_v,opt_varname="diffusion: z_w_v") + CALL sync_patch_array(SYNC_C,p_patch(jg),theta_v_ic,opt_varname="diffusion: theta_v_ic") + + fac2d = 0.0625_wp ! Factor of the 2D deformation field which is used as minimum of the 3D def field + +!$OMP PARALLEL PRIVATE(i_startblk,i_endblk) + + i_startblk = p_patch_diff(jg)%edges%start_block(rl_start) + i_endblk = p_patch_diff(jg)%edges%end_block(rl_end) + +!$OMP DO PRIVATE(jb,i_startidx,i_endidx,jk,je,vn_vert1,vn_vert2,vn_vert3,vn_vert4,dvt_norm,dvt_tang, & +!$OMP z_vn_ie,z_vt_ie,dvndz,dvtdz,dwdz,dthvdz,dwdn,dwdt,kh_smag3d_e), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 2, nlev + DO je = i_startidx, i_endidx + z_vn_ie(je,jk) = diff_inst(jg)%wgtfac_e(je,jk,jb)*vn(je,jk,jb) + & + (1._wp - diff_inst(jg)%wgtfac_e(je,jk,jb))*vn(je,jk-1,jb) + z_vt_ie(je,jk) = diff_inst(jg)%wgtfac_e(je,jk,jb)*vt(je,jk,jb) + & + (1._wp - diff_inst(jg)%wgtfac_e(je,jk,jb))*vt(je,jk-1,jb) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR ASYNC(1) IF(i_am_accel_node) + DO je = i_startidx, i_endidx + z_vn_ie(je,1) = & + diff_inst(jg)%wgtfacq1_e(je,1,jb)*vn(je,1,jb) + & + diff_inst(jg)%wgtfacq1_e(je,2,jb)*vn(je,2,jb) + & + diff_inst(jg)%wgtfacq1_e(je,3,jb)*vn(je,3,jb) + z_vn_ie(je,nlevp1) = & + diff_inst(jg)%wgtfacq_e(je,1,jb)*vn(je,nlev,jb) + & + diff_inst(jg)%wgtfacq_e(je,2,jb)*vn(je,nlev-1,jb) + & + diff_inst(jg)%wgtfacq_e(je,3,jb)*vn(je,nlev-2,jb) + z_vt_ie(je,1) = & + diff_inst(jg)%wgtfacq1_e(je,1,jb)*vt(je,1,jb) + & + diff_inst(jg)%wgtfacq1_e(je,2,jb)*vt(je,2,jb) + & + diff_inst(jg)%wgtfacq1_e(je,3,jb)*vt(je,3,jb) + z_vt_ie(je,nlevp1) = & + diff_inst(jg)%wgtfacq_e(je,1,jb)*vt(je,nlev,jb) + & + diff_inst(jg)%wgtfacq_e(je,2,jb)*vt(je,nlev-1,jb) + & + diff_inst(jg)%wgtfacq_e(je,3,jb)*vt(je,nlev-2,jb) + ENDDO + !$ACC END PARALLEL LOOP + + ! Computation of wind field deformation + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO je = i_startidx, i_endidx +!DIR$ IVDEP + DO jk = 1, nlev +#else +!$NEC outerloop_unroll(4) + DO jk = 1, nlev + DO je = i_startidx, i_endidx +#endif + + vn_vert1 = u_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v1 + & + v_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v2 + + vn_vert2 = u_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v1 + & + v_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v2 + + dvt_tang = p_patch_diff(jg)%edges%tangent_orientation(je,jb)* ( & + u_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,2)%v1 + & + v_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,2)%v2 - & + (u_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,1)%v1 + & + v_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,1)%v2) ) + + vn_vert3 = u_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,3)%v1 + & + v_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,3)%v2 + + vn_vert4 = u_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,4)%v1 + & + v_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,4)%v2 + + dvt_norm = u_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,4)%v1 + & + v_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,4)%v2 - & + (u_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,3)%v1 + & + v_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,3)%v2) + + dvndz(je,jk) = (z_vn_ie(je,jk) - z_vn_ie(je,jk+1)) / diff_inst(jg)%ddqz_z_full_e(je,jk,jb) + dvtdz(je,jk) = (z_vt_ie(je,jk) - z_vt_ie(je,jk+1)) / diff_inst(jg)%ddqz_z_full_e(je,jk,jb) + + dwdz (je,jk) = & + (diff_inst(jg)%c_lin_e(je,1,jb)*(w(iecidx(je,jb,1),jk, iecblk(je,jb,1)) - & + w(iecidx(je,jb,1),jk+1,iecblk(je,jb,1)) ) + & + diff_inst(jg)%c_lin_e(je,2,jb)*(w(iecidx(je,jb,2),jk, iecblk(je,jb,2)) - & + w(iecidx(je,jb,2),jk+1,iecblk(je,jb,2)) ) ) / & + diff_inst(jg)%ddqz_z_full_e(je,jk,jb) + + dthvdz(je,jk) = & + (diff_inst(jg)%c_lin_e(je,1,jb)*(theta_v_ic(iecidx(je,jb,1),jk, iecblk(je,jb,1)) - & + theta_v_ic(iecidx(je,jb,1),jk+1,iecblk(je,jb,1)) ) + & + diff_inst(jg)%c_lin_e(je,2,jb)*(theta_v_ic(iecidx(je,jb,2),jk, iecblk(je,jb,2)) - & + theta_v_ic(iecidx(je,jb,2),jk+1,iecblk(je,jb,2)) ) ) / & + diff_inst(jg)%ddqz_z_full_e(je,jk,jb) + + dwdn (je,jk) = p_patch_diff(jg)%edges%inv_dual_edge_length(je,jb)* ( & + 0.5_wp*(w(iecidx(je,jb,1),jk, iecblk(je,jb,1)) + & + w(iecidx(je,jb,1),jk+1,iecblk(je,jb,1))) - & + 0.5_wp*(w(iecidx(je,jb,2),jk, iecblk(je,jb,2)) + & + w(iecidx(je,jb,2),jk+1,iecblk(je,jb,2))) ) + + dwdt (je,jk) = p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) * & + p_patch_diff(jg)%edges%tangent_orientation(je,jb) * ( & + 0.5_wp*(z_w_v(ividx(je,jb,1),jk,ivblk(je,jb,1))+z_w_v(ividx(je,jb,1),jk+1,ivblk(je,jb,1))) - & + 0.5_wp*(z_w_v(ividx(je,jb,2),jk,ivblk(je,jb,2))+z_w_v(ividx(je,jb,2),jk+1,ivblk(je,jb,2))) ) + + kh_smag3d_e(je,jk) = 2._wp*( & + ( (vn_vert4-vn_vert3)*p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb) )**2 + & + (dvt_tang*p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb))**2 + dwdz(je,jk)**2) + & + 0.5_wp *( (p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) * & + p_patch_diff(jg)%edges%tangent_orientation(je,jb)*(vn_vert2-vn_vert1) + & + p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb)*dvt_norm )**2 + & + (dvndz(je,jk) + dwdn(je,jk))**2 + (dvtdz(je,jk) + dwdt(je,jk))**2 ) - & + 3._wp*diff_inst(jg)%grav * dthvdz(je,jk) / ( & + diff_inst(jg)%c_lin_e(je,1,jb)*theta_v(iecidx(je,jb,1),jk,iecblk(je,jb,1)) + & + diff_inst(jg)%c_lin_e(je,2,jb)*theta_v(iecidx(je,jb,2),jk,iecblk(je,jb,2)) ) + + ! 2D Smagorinsky diffusion coefficient + kh_smag_e(je,jk,jb) = diff_multfac_smag(jk)*SQRT( MAX(kh_smag3d_e(je,jk), fac2d*( & + ((vn_vert4-vn_vert3)*p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb)- & + dvt_tang*p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) )**2 + ( & + (vn_vert2-vn_vert1)*p_patch_diff(jg)%edges%tangent_orientation(je,jb)* & + p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) + & + dvt_norm*p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb) )**2 ) ) ) + + ! The factor of 4 comes from dividing by twice the "correct" length + z_nabla2_e(je,jk,jb) = 4._wp * ( & + (vn_vert4 + vn_vert3 - 2._wp*vn(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb)**2 + & + (vn_vert2 + vn_vert1 - 2._wp*vn(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb)**2 ) + + kh_smag_ec(je,jk,jb) = kh_smag_e(je,jk,jb) + ! Subtract part of the fourth-order background diffusion coefficient + kh_smag_e(je,jk,jb) = MAX(0._vp,kh_smag_e(je,jk,jb) - smag_offset) + ! Limit diffusion coefficient to the theoretical CFL stability threshold + kh_smag_e(je,jk,jb) = MIN(kh_smag_e(je,jk,jb),smag_limit(jk)) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ENDDO +!$OMP END DO NOWAIT +!$OMP END PARALLEL + + ELSE IF ((diffu_type == 3 .OR. diffu_type == 5) .AND. diff_inst(jg)%discr_vn >= 2) THEN + + ! RBF reconstruction of velocity at vertices and cells + CALL rbf_vec_interpol_vertex( vn, p_patch_diff(jg), & + diff_inst(jg)%rbf_vec_idx_v, diff_inst(jg)%rbf_vec_blk_v, & + diff_inst(jg)%rbf_vec_coeff_v, u_vert, v_vert, & + opt_rlend=min_rlvert_int-1, opt_acc_async=.TRUE. ) + + ! DA: This wait ideally should be removed + !$ACC WAIT + + IF (diff_inst(jg)%discr_vn == 2) THEN + + CALL rbf_vec_interpol_cell( vn, p_patch_diff(jg), & + diff_inst(jg)%rbf_vec_idx_v, diff_inst(jg)%rbf_vec_blk_v, & + diff_inst(jg)%rbf_vec_coeff_v, u_cell, v_cell, & + opt_rlend=min_rlcell_int-1 ) + ELSE + + CALL edges2cells_vector( vn, vt, p_patch_diff(jg), diff_inst(jg)%e_bln_c_u, diff_inst(jg)%e_bln_c_v, & + u_cell, v_cell, opt_rlend=min_rlcell_int-1 ) + ENDIF + + IF (diff_inst(jg)%p_test_run) THEN + !$ACC KERNELS IF(i_am_accel_node) ASYNC(1) + z_nabla2_e = 0._wp + !$ACC END KERNELS + ENDIF + +!$OMP PARALLEL PRIVATE(i_startblk,i_endblk,rl_start,rl_end) + + rl_start = start_bdydiff_e + rl_end = min_rledge_int - 1 + + i_startblk = p_patch_diff(jg)%edges%start_block(rl_start) + i_endblk = p_patch_diff(jg)%edges%end_block(rl_end) + +!$OMP DO PRIVATE(jb,i_startidx,i_endidx,jk,je,vn_vert1,vn_vert2,vn_cell1,vn_cell2,& +!$OMP dvt_norm,dvt_tang), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + ! Computation of wind field deformation + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO je = i_startidx, i_endidx +!DIR$ IVDEP + DO jk = 1, nlev +#else +!$NEC outerloop_unroll(4) + DO jk = 1, nlev + DO je = i_startidx, i_endidx +#endif + + vn_vert1 = u_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v1 + & + v_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v2 + + vn_vert2 = u_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v1 + & + v_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v2 + + dvt_tang = p_patch_diff(jg)%edges%tangent_orientation(je,jb)* ( & + u_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,2)%v1 + & + v_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,2)%v2 - & + (u_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,1)%v1 + & + v_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_vert(je,jb,1)%v2) ) + + vn_cell1 = u_cell(iecidx(je,jb,1),jk,iecblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_cell(je,jb,1)%v1 + & + v_cell(iecidx(je,jb,1),jk,iecblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_cell(je,jb,1)%v2 + + vn_cell2 = u_cell(iecidx(je,jb,2),jk,iecblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_cell(je,jb,2)%v1 + & + v_cell(iecidx(je,jb,2),jk,iecblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_cell(je,jb,2)%v2 + + dvt_norm = u_cell(iecidx(je,jb,2),jk,iecblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_cell(je,jb,2)%v1 + & + v_cell(iecidx(je,jb,2),jk,iecblk(je,jb,2)) * & + p_patch_diff(jg)%edges%dual_normal_cell(je,jb,2)%v2 - & + (u_cell(iecidx(je,jb,1),jk,iecblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_cell(je,jb,1)%v1 + & + v_cell(iecidx(je,jb,1),jk,iecblk(je,jb,1)) * & + p_patch_diff(jg)%edges%dual_normal_cell(je,jb,1)%v2) + + + ! Smagorinsky diffusion coefficient + kh_smag_e(je,jk,jb) = diff_multfac_smag(jk)*SQRT( ( & + (vn_cell2-vn_cell1)*p_patch_diff(jg)%edges%inv_dual_edge_length(je,jb)- & + dvt_tang*p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) )**2 + ( & + (vn_vert2-vn_vert1)*p_patch_diff(jg)%edges%tangent_orientation(je,jb)* & + p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb) + & + dvt_norm*p_patch_diff(jg)%edges%inv_dual_edge_length(je,jb))**2 ) + + ! The factor of 4 comes from dividing by twice the "correct" length + z_nabla2_e(je,jk,jb) = 4._wp * ( & + (vn_cell2 + vn_cell1 - 2._wp*vn(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_dual_edge_length(je,jb)**2 + & + (vn_vert2 + vn_vert1 - 2._wp*vn(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb)**2 ) + +#ifndef _OPENACC + ENDDO + ENDDO + + DO jk = 1, nlev + DO je = i_startidx, i_endidx +#endif + + kh_smag_ec(je,jk,jb) = kh_smag_e(je,jk,jb) + ! Subtract part of the fourth-order background diffusion coefficient + kh_smag_e(je,jk,jb) = MAX(0._vp,kh_smag_e(je,jk,jb) - smag_offset) + ! Limit diffusion coefficient to the theoretical CFL stability threshold + kh_smag_e(je,jk,jb) = MIN(kh_smag_e(je,jk,jb),smag_limit(jk)) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ENDDO +!$OMP END DO NOWAIT +!$OMP END PARALLEL + + ENDIF + + ! Compute input quantities for turbulence scheme + IF ((diffu_type == 3 .OR. diffu_type == 5) .AND. & + (diff_inst(jg)%itype_sher >= 1 .OR. diff_inst(jg)%ltkeshs)) THEN + +!$OMP PARALLEL PRIVATE(i_startblk,i_endblk) + rl_start = grf_bdywidth_c+1 + rl_end = min_rlcell_int + + i_startblk = p_patch_diff(jg)%cells%start_block(rl_start) + i_endblk = p_patch_diff(jg)%cells%end_block(rl_end) + +!$OMP DO PRIVATE(jk,jc,jb,i_startidx,i_endidx,kh_c,div), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_c(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO jc = i_startidx, i_endidx + DO jk = 1, nlev +#else + DO jk = 1, nlev + DO jc = i_startidx, i_endidx +#endif + + kh_c(jc,jk) = (kh_smag_ec(ieidx(jc,jb,1),jk,ieblk(jc,jb,1))*diff_inst(jg)%e_bln_c_s(jc,1,jb) + & + kh_smag_ec(ieidx(jc,jb,2),jk,ieblk(jc,jb,2))*diff_inst(jg)%e_bln_c_s(jc,2,jb) + & + kh_smag_ec(ieidx(jc,jb,3),jk,ieblk(jc,jb,3))*diff_inst(jg)%e_bln_c_s(jc,3,jb))/ & + diff_multfac_smag(jk) + + div(jc,jk) = vn(ieidx(jc,jb,1),jk,ieblk(jc,jb,1))*diff_inst(jg)%geofac_div(jc,1,jb) + & + vn(ieidx(jc,jb,2),jk,ieblk(jc,jb,2))*diff_inst(jg)%geofac_div(jc,2,jb) + & + vn(ieidx(jc,jb,3),jk,ieblk(jc,jb,3))*diff_inst(jg)%geofac_div(jc,3,jb) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 2, nlev ! levels 1 and nlevp1 are unused +!DIR$ IVDEP + DO jc = i_startidx, i_endidx + + div_ic(jc,jk,jb) = diff_inst(jg)%wgtfac_c(jc,jk,jb)*div(jc,jk) + & + (1._wp-diff_inst(jg)%wgtfac_c(jc,jk,jb))*div(jc,jk-1) + + hdef_ic(jc,jk,jb) = (diff_inst(jg)%wgtfac_c(jc,jk,jb)*kh_c(jc,jk) + & + (1._wp-diff_inst(jg)%wgtfac_c(jc,jk,jb))*kh_c(jc,jk-1))**2 + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO +!$OMP END PARALLEL + + ENDIF + + IF (diffu_type == 5) THEN ! Add fourth-order background diffusion + + IF (diff_inst(jg)%discr_vn > 1) THEN + CALL sync_patch_array(SYNC_E,p_patch(jg),z_nabla2_e, & + opt_varname="diffusion: nabla2_e") + END IF + + ! Interpolate nabla2(v) to vertices in order to compute nabla2(nabla2(v)) + + IF (diff_inst(jg)%p_test_run) THEN + !$ACC KERNELS IF(i_am_accel_node) + u_vert = 0._wp + v_vert = 0._wp + !$ACC END KERNELS + ENDIF + + CALL rbf_vec_interpol_vertex( z_nabla2_e, p_patch_diff(jg), & + diff_inst(jg)%rbf_vec_idx_v, diff_inst(jg)%rbf_vec_blk_v, & + diff_inst(jg)%rbf_vec_coeff_v, u_vert, v_vert, & + opt_rlstart=4, opt_rlend=min_rlvert_int, opt_acc_async=.TRUE. ) + rl_start = grf_bdywidth_e+1 + rl_end = min_rledge_int + + IF (diff_inst(jg)%itype_comm == 1 .OR. diff_inst(jg)%itype_comm == 3) THEN +#ifdef __MIXED_PRECISION + CALL sync_patch_array_mult_mp(SYNC_V,p_patch(jg),0,2,f3din1_sp=u_vert,f3din2_sp=v_vert, & + opt_varname="diffusion: u_vert and v_vert 3") +#else + CALL sync_patch_array_mult(SYNC_V,p_patch(jg),2,u_vert,v_vert, & + opt_varname="diffusion: u_vert and v_vert 3") +#endif + ENDIF + +!$OMP PARALLEL PRIVATE(i_startblk,i_endblk) + + i_startblk = p_patch_diff(jg)%edges%start_block(rl_start) + i_endblk = p_patch_diff(jg)%edges%end_block(rl_end) + +!$OMP DO PRIVATE(jb,i_startidx,i_endidx,jk,je,nabv_tang,nabv_norm,z_nabla4_e2,z_d_vn_hdf), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + ! Compute nabla4(v) + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO je = i_startidx, i_endidx + DO jk = 1, nlev +#else +!$NEC outerloop_unroll(4) + DO jk = 1, nlev + DO je = i_startidx, i_endidx +#endif + + nabv_tang = u_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v1 + & + v_vert(ividx(je,jb,1),jk,ivblk(je,jb,1)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,1)%v2 + & + u_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v1 + & + v_vert(ividx(je,jb,2),jk,ivblk(je,jb,2)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,2)%v2 + + nabv_norm = u_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,3)%v1 + & + v_vert(ividx(je,jb,3),jk,ivblk(je,jb,3)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,3)%v2 + & + u_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,4)%v1 + & + v_vert(ividx(je,jb,4),jk,ivblk(je,jb,4)) * & + p_patch_diff(jg)%edges%primal_normal_vert(je,jb,4)%v2 + + ! The factor of 4 comes from dividing by twice the "correct" length + z_nabla4_e2(je,jk) = 4._wp * ( & + (nabv_norm - 2._wp*z_nabla2_e(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_vert_vert_length(je,jb)**2 + & + (nabv_tang - 2._wp*z_nabla2_e(je,jk,jb)) & + *p_patch_diff(jg)%edges%inv_primal_edge_length(je,jb)**2 ) + + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ! Apply diffusion for the case of diffu_type = 5 + IF ( jg == 1 .AND. diff_inst(jg)%l_limited_area .OR. jg > 1 .AND. .NOT. diff_inst(jg)%lfeedback ) THEN + ! + ! Domains with lateral boundary and nests without feedback + ! + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(z_d_vn_hdf) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO je = i_startidx, i_endidx + ! + z_d_vn_hdf = p_patch_diff(jg)%edges%area_edge(je,jb) & + & * ( MAX(nudgezone_diff*diff_inst(jg)%nudgecoeff_e(je,jb), & + & REAL(kh_smag_e(je,jk,jb),wp)) * z_nabla2_e(je,jk,jb) & + & - p_patch_diff(jg)%edges%area_edge(je,jb) & + & * diff_multfac_vn(jk) * z_nabla4_e2(je,jk) ) + ! + vn(je,jk,jb) = vn(je,jk,jb) + z_d_vn_hdf + ! +#ifdef __ENABLE_DDT_VN_XYZ__ + IF ( diff_inst(jg)%ddt_vn_hdf_is_associated) THEN + ddt_vn_hdf(je,jk,jb) = ddt_vn_hdf(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF + ! + IF ( diff_inst(jg)%ddt_vn_dyn_is_associated) THEN + ddt_vn_dyn(je,jk,jb) = ddt_vn_dyn(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF +#endif + ! + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ELSE IF (jg > 1) THEN + ! + ! Nests with feedback + ! + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(z_d_vn_hdf) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO je = i_startidx, i_endidx + ! + z_d_vn_hdf = p_patch_diff(jg)%edges%area_edge(je,jb) & + & * ( kh_smag_e(je,jk,jb) * z_nabla2_e(je,jk,jb) & + & - p_patch_diff(jg)%edges%area_edge(je,jb) & + & * MAX(diff_multfac_vn(jk),bdy_diff*diff_inst(jg)%nudgecoeff_e(je,jb)) * z_nabla4_e2(je,jk) ) + ! + vn(je,jk,jb) = vn(je,jk,jb) + z_d_vn_hdf + ! +#ifdef __ENABLE_DDT_VN_XYZ__ + IF ( diff_inst(jg)%ddt_vn_hdf_is_associated) THEN + ddt_vn_hdf(je,jk,jb) = ddt_vn_hdf(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF + ! + IF ( diff_inst(jg)%ddt_vn_dyn_is_associated) THEN + ddt_vn_dyn(je,jk,jb) = ddt_vn_dyn(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF +#endif + ! + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ELSE + + ! + ! Global domains + ! + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(z_d_vn_hdf) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO je = i_startidx, i_endidx + ! + z_d_vn_hdf = p_patch_diff(jg)%edges%area_edge(je,jb) & + & * ( kh_smag_e(je,jk,jb) * z_nabla2_e(je,jk,jb) & + & - p_patch_diff(jg)%edges%area_edge(je,jb) & + & * diff_multfac_vn(jk) * z_nabla4_e2(je,jk) ) + + ! + vn(je,jk,jb) = vn(je,jk,jb) + z_d_vn_hdf + ! +#ifdef __ENABLE_DDT_VN_XYZ__ + IF ( diff_inst(jg)%ddt_vn_hdf_is_associated) THEN + ddt_vn_hdf(je,jk,jb) = ddt_vn_hdf(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF + ! + IF ( diff_inst(jg)%ddt_vn_dyn_is_associated) THEN + ddt_vn_dyn(je,jk,jb) = ddt_vn_dyn(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF +#endif + ! + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDIF + + ENDDO +!$OMP END DO NOWAIT +!$OMP END PARALLEL + + ENDIF + + ! Apply diffusion for the cases of diffu_type = 3 or 4 + +!$OMP PARALLEL PRIVATE(i_startblk,i_endblk,rl_start,rl_end) + + rl_start = grf_bdywidth_e+1 + rl_end = min_rledge_int + + i_startblk = p_patch_diff(jg)%edges%start_block(rl_start) + i_endblk = p_patch_diff(jg)%edges%end_block(rl_end) + + IF (diffu_type == 3) THEN ! Only Smagorinsky diffusion + IF ( jg == 1 .AND. diff_inst(jg)%l_limited_area .OR. jg > 1 .AND. .NOT. diff_inst(jg)%lfeedback ) THEN + +!$OMP DO PRIVATE(jb,i_startidx,i_endidx,jk,je,z_d_vn_hdf) ICON_OMP_DEFAULT_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(z_d_vn_hdf) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO je = i_startidx, i_endidx + ! + z_d_vn_hdf = p_patch_diff(jg)%edges%area_edge(je,jb) & + & * MAX(nudgezone_diff*diff_inst(jg)%nudgecoeff_e(je,jb),REAL(kh_smag_e(je,jk,jb),wp)) & + & * z_nabla2_e(je,jk,jb) + ! + vn(je,jk,jb) = vn(je,jk,jb) + z_d_vn_hdf + ! +#ifdef __ENABLE_DDT_VN_XYZ__ + IF ( diff_inst(jg)%ddt_vn_hdf_is_associated) THEN + ddt_vn_hdf(je,jk,jb) = ddt_vn_hdf(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF + ! + IF ( diff_inst(jg)%ddt_vn_dyn_is_associated) THEN + ddt_vn_dyn(je,jk,jb) = ddt_vn_dyn(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF +#endif + ! + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + ELSE + +!$OMP DO PRIVATE(jb,i_startidx,i_endidx,jk,je,z_d_vn_hdf) ICON_OMP_DEFAULT_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(z_d_vn_hdf) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO je = i_startidx, i_endidx + ! + z_d_vn_hdf = p_patch_diff(jg)%edges%area_edge(je,jb) * kh_smag_e(je,jk,jb) * z_nabla2_e(je,jk,jb) + ! + vn(je,jk,jb) = vn(je,jk,jb) + z_d_vn_hdf + ! +#ifdef __ENABLE_DDT_VN_XYZ__ + IF ( diff_inst(jg)%ddt_vn_hdf_is_associated) THEN + ddt_vn_hdf(je,jk,jb) = ddt_vn_hdf(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF + ! + IF ( diff_inst(jg)%ddt_vn_dyn_is_associated) THEN + ddt_vn_dyn(je,jk,jb) = ddt_vn_dyn(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF +#endif + ! + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + ENDIF + + ELSE IF (diffu_type == 4) THEN + +!$OMP DO PRIVATE(jb,i_startidx,i_endidx,jk,je,z_d_vn_hdf) ICON_OMP_DEFAULT_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(z_d_vn_hdf) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO je = i_startidx, i_endidx + ! + z_d_vn_hdf = - p_patch_diff(jg)%edges%area_edge(je,jb)*p_patch_diff(jg)%edges%area_edge(je,jb) & + & * diff_multfac_vn(jk) * z_nabla4_e(je,jk,jb) + ! + vn(je,jk,jb) = vn(je,jk,jb) + z_d_vn_hdf + ! +#ifdef __ENABLE_DDT_VN_XYZ__ + IF ( diff_inst(jg)%ddt_vn_hdf_is_associated) THEN + ddt_vn_hdf(je,jk,jb) = ddt_vn_hdf(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF + ! + IF ( diff_inst(jg)%ddt_vn_dyn_is_associated) THEN + ddt_vn_dyn(je,jk,jb) = ddt_vn_dyn(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF +#endif + ! + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + ENDIF + + IF (diff_inst(jg)%l_limited_area .OR. jg > 1) THEN + + ! Lateral boundary diffusion for vn + i_startblk = p_patch_diff(jg)%edges%start_block(start_bdydiff_e) + i_endblk = p_patch_diff(jg)%edges%end_block(grf_bdywidth_e) + +!$OMP DO PRIVATE(je,jk,jb,i_startidx,i_endidx,z_d_vn_hdf) ICON_OMP_DEFAULT_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, start_bdydiff_e, grf_bdywidth_e) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(z_d_vn_hdf) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO je = i_startidx, i_endidx + ! + z_d_vn_hdf = p_patch_diff(jg)%edges%area_edge(je,jb) * fac_bdydiff_v * z_nabla2_e(je,jk,jb) + ! + vn(je,jk,jb) = vn(je,jk,jb) + z_d_vn_hdf + ! +#ifdef __ENABLE_DDT_VN_XYZ__ + IF ( diff_inst(jg)%ddt_vn_hdf_is_associated) THEN + ddt_vn_hdf(je,jk,jb) = ddt_vn_hdf(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF + ! + IF ( diff_inst(jg)%ddt_vn_dyn_is_associated) THEN + ddt_vn_dyn(je,jk,jb) = ddt_vn_dyn(je,jk,jb) + z_d_vn_hdf * r_dtimensubsteps + END IF +#endif + ! + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + ENDIF ! vn boundary diffusion + + IF (diff_inst(jg)%lhdiff_rcf .AND. diff_inst(jg)%lhdiff_w) THEN ! add diffusion on vertical wind speed + ! remark: the surface level (nlevp1) is excluded because w is diagnostic there + + IF (diff_inst(jg)%l_limited_area .AND. jg == 1) THEN + rl_start = grf_bdywidth_c+1 + ELSE + rl_start = grf_bdywidth_c + ENDIF + rl_end = min_rlcell_int-1 + + i_startblk = p_patch_diff(jg)%cells%start_block(rl_start) + i_endblk = p_patch_diff(jg)%cells%end_block(rl_end) + +!$OMP DO PRIVATE(jk,jc,jb,i_startidx,i_endidx), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_c(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO jc = i_startidx, i_endidx +!DIR$ IVDEP +#ifdef _CRAYFTN +!DIR$ PREFERVECTOR +#endif + DO jk = 1, nlev +#else + DO jk = 1, nlev + DO jc = i_startidx, i_endidx +#endif + z_nabla2_c(jc,jk,jb) = & + w(jc,jk,jb) *diff_inst(jg)%geofac_n2s(jc,1,jb) + & + w(icidx(jc,jb,1),jk,icblk(jc,jb,1))*diff_inst(jg)%geofac_n2s(jc,2,jb) + & + w(icidx(jc,jb,2),jk,icblk(jc,jb,2))*diff_inst(jg)%geofac_n2s(jc,3,jb) + & + w(icidx(jc,jb,3),jk,icblk(jc,jb,3))*diff_inst(jg)%geofac_n2s(jc,4,jb) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + IF (diff_inst(jg)%itype_sher == 2) THEN ! compute horizontal gradients of w + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO jc = i_startidx, i_endidx +!DIR$ IVDEP + DO jk = 2, nlev +#else + DO jk = 2, nlev + DO jc = i_startidx, i_endidx +#endif + dwdx(jc,jk,jb) = diff_inst(jg)%geofac_grg(jc,1,jb,1)*w(jc,jk,jb) + & + diff_inst(jg)%geofac_grg(jc,2,jb,1)*w(icidx(jc,jb,1),jk,icblk(jc,jb,1)) + & + diff_inst(jg)%geofac_grg(jc,3,jb,1)*w(icidx(jc,jb,2),jk,icblk(jc,jb,2)) + & + diff_inst(jg)%geofac_grg(jc,4,jb,1)*w(icidx(jc,jb,3),jk,icblk(jc,jb,3)) + + dwdy(jc,jk,jb) = diff_inst(jg)%geofac_grg(jc,1,jb,2)*w(jc,jk,jb) + & + diff_inst(jg)%geofac_grg(jc,2,jb,2)*w(icidx(jc,jb,1),jk,icblk(jc,jb,1)) + & + diff_inst(jg)%geofac_grg(jc,3,jb,2)*w(icidx(jc,jb,2),jk,icblk(jc,jb,2)) + & + diff_inst(jg)%geofac_grg(jc,4,jb,2)*w(icidx(jc,jb,3),jk,icblk(jc,jb,3)) + + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDIF + + ENDDO +!$OMP END DO + + IF (diff_inst(jg)%l_limited_area .AND. jg == 1) THEN + rl_start = 0 + ELSE + rl_start = grf_bdywidth_c+1 + ENDIF + rl_end = min_rlcell_int + + i_startblk = p_patch_diff(jg)%cells%start_block(rl_start) + i_endblk = p_patch_diff(jg)%cells%end_block(rl_end) + +!$OMP DO PRIVATE(jk,jc,jb,i_startidx,i_endidx), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_c(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO jc = i_startidx, i_endidx +!DIR$ IVDEP + DO jk = 1, nlev +#else + DO jk = 1, nlev + DO jc = i_startidx, i_endidx +#endif + w(jc,jk,jb) = w(jc,jk,jb) - diff_multfac_w * p_patch_diff(jg)%cells%area(jc,jb)**2 * & + (z_nabla2_c(jc,jk,jb) *diff_inst(jg)%geofac_n2s(jc,1,jb) + & + z_nabla2_c(icidx(jc,jb,1),jk,icblk(jc,jb,1))*diff_inst(jg)%geofac_n2s(jc,2,jb) + & + z_nabla2_c(icidx(jc,jb,2),jk,icblk(jc,jb,2))*diff_inst(jg)%geofac_n2s(jc,3,jb) + & + z_nabla2_c(icidx(jc,jb,3),jk,icblk(jc,jb,3))*diff_inst(jg)%geofac_n2s(jc,4,jb)) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ! Add nabla2 diffusion in upper damping layer (if present) + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 2, diff_inst(jg)%nrdmax +!DIR$ IVDEP + DO jc = i_startidx, i_endidx + w(jc,jk,jb) = w(jc,jk,jb) + & + diff_multfac_n2w(jk) * p_patch_diff(jg)%cells%area(jc,jb) * z_nabla2_c(jc,jk,jb) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ENDDO +!$OMP END DO + + ENDIF ! w diffusion + +!$OMP END PARALLEL + + IF (diff_inst(jg)%itype_comm == 1 .OR. diff_inst(jg)%itype_comm == 3) THEN + CALL sync_patch_array(SYNC_E, p_patch(jg), vn,opt_varname="diffusion: vn sync") + ENDIF + + IF (ltemp_diffu) THEN ! Smagorinsky temperature diffusion + + IF (diff_inst(jg)%l_zdiffu_t) THEN + icell => diff_inst(jg)%zd_indlist + iblk => diff_inst(jg)%zd_blklist + ilev => diff_inst(jg)%zd_vertidx + ! iedge => diff_inst(jg)%zd_edgeidx + ! iedblk => diff_inst(jg)%zd_edgeblk + vcoef => diff_inst(jg)%zd_intcoef + ! blcoef => diff_inst(jg)%zd_e2cell + zd_geofac => diff_inst(jg)%zd_geofac + +!!! nproma_zdiffu = cpu_min_nproma(diff_inst(jg)%nproma,256) +#ifdef _OPENACC + nproma_zdiffu = diff_inst(jg)%nproma +#else + nproma_zdiffu = MIN(diff_inst(jg)%nproma,256) +#endif + nblks_zdiffu = INT(diff_inst(jg)%zd_listdim/nproma_zdiffu) + npromz_zdiffu = MOD(diff_inst(jg)%zd_listdim,nproma_zdiffu) + IF (npromz_zdiffu > 0) THEN + nblks_zdiffu = nblks_zdiffu + 1 + ELSE + npromz_zdiffu = nproma_zdiffu + ENDIF + ENDIF + +!$OMP PARALLEL PRIVATE(rl_start,rl_end,i_startblk,i_endblk) + + ! Enhance Smagorinsky diffusion coefficient in the presence of excessive grid-point cold pools + ! This is restricted to the two lowest model levels + ! + rl_start = grf_bdywidth_c + rl_end = min_rlcell_int-1 + + i_startblk = p_patch_diff(jg)%cells%start_block(rl_start) + i_endblk = p_patch_diff(jg)%cells%end_block(rl_end) + +!$OMP DO PRIVATE(jk,jc,jb,i_startidx,i_endidx,ic,tdiff,trefdiff), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_c(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + ic = 0 + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = nlev-1, nlev +!DIR$ IVDEP + DO jc = i_startidx, i_endidx + ! Perturbation potential temperature difference between local point and average of the three neighbors + tdiff = theta_v(jc,jk,jb) - & + (theta_v(icidx(jc,jb,1),jk,icblk(jc,jb,1)) + & + theta_v(icidx(jc,jb,2),jk,icblk(jc,jb,2)) + & + theta_v(icidx(jc,jb,3),jk,icblk(jc,jb,3)) ) / 3._wp + trefdiff = diff_inst(jg)%theta_ref_mc(jc,jk,jb) - & + (diff_inst(jg)%theta_ref_mc(icidx(jc,jb,1),jk,icblk(jc,jb,1)) + & + diff_inst(jg)%theta_ref_mc(icidx(jc,jb,2),jk,icblk(jc,jb,2)) + & + diff_inst(jg)%theta_ref_mc(icidx(jc,jb,3),jk,icblk(jc,jb,3)) ) / 3._wp + + ! Enahnced horizontal diffusion is applied if the theta perturbation is either + ! - at least 5 K colder than the average of the neighbor points on valley points (determined by trefdiff < 0.) or + ! - at least 7.5 K colder than the average of the neighbor points otherwise + IF (tdiff-trefdiff < thresh_tdiff .AND. trefdiff < 0._wp .OR. tdiff-trefdiff < 1.5_wp*thresh_tdiff) THEN +#ifndef _OPENACC + ic = ic+1 + iclist(ic,jb) = jc + iklist(ic,jb) = jk + tdlist(ic,jb) = thresh_tdiff - tdiff + trefdiff +#else + ! Enhance Smagorinsky coefficients at the three edges of the cells included in the list +! Attention: this operation is neither vectorizable nor OpenMP-parallelizable (race conditions!) + enh_diffu_3d(jc,jk,jb) = (thresh_tdiff - tdiff + trefdiff)*5.e-4_vp + ELSE + enh_diffu_3d(jc,jk,jb) = -HUGE(0._vp) ! In order that this is never taken as the MAX +#endif + ENDIF + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + icount(jb) = ic + + ENDDO +!$OMP END DO + + ! Enhance Smagorinsky coefficients at the three edges of the cells included in the list + ! Attention: this operation is neither vectorizable nor OpenMP-parallelizable (race conditions!) + +#ifndef _OPENACC +!$OMP MASTER + DO jb = i_startblk,i_endblk + + IF (icount(jb) > 0) THEN + DO ic = 1, icount(jb) + jc = iclist(ic,jb) + jk = iklist(ic,jb) + enh_diffu = tdlist(ic,jb)*5.e-4_vp + kh_smag_e(ieidx(jc,jb,1),jk,ieblk(jc,jb,1)) = MAX(enh_diffu,kh_smag_e(ieidx(jc,jb,1),jk,ieblk(jc,jb,1))) + kh_smag_e(ieidx(jc,jb,2),jk,ieblk(jc,jb,2)) = MAX(enh_diffu,kh_smag_e(ieidx(jc,jb,2),jk,ieblk(jc,jb,2))) + kh_smag_e(ieidx(jc,jb,3),jk,ieblk(jc,jb,3)) = MAX(enh_diffu,kh_smag_e(ieidx(jc,jb,3),jk,ieblk(jc,jb,3))) + ENDDO + ENDIF + + ENDDO + +!$OMP END MASTER +!$OMP BARRIER + +#else + + rl_start = grf_bdywidth_e+1 + rl_end = min_rledge_int + + i_startblk = p_patch_diff(jg)%edges%start_block(rl_start) + i_endblk = p_patch_diff(jg)%edges%end_block(rl_end) + + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = nlev-1, nlev + DO je = i_startidx, i_endidx + kh_smag_e(je,jk,jb) = MAX(kh_smag_e(je,jk,jb), enh_diffu_3d(iecidx(je,jb,1),jk,iecblk(je,jb,1)), & + enh_diffu_3d(iecidx(je,jb,2),jk,iecblk(je,jb,2)) ) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +#endif + + IF (diff_inst(jg)%discr_t == 1) THEN ! use discretization K*nabla(theta) + + rl_start = grf_bdywidth_c+1 + rl_end = min_rlcell_int + + i_startblk = p_patch_diff(jg)%cells%start_block(rl_start) + i_endblk = p_patch_diff(jg)%cells%end_block(rl_end) + +!$OMP DO PRIVATE(jk,jc,jb,i_startidx,i_endidx), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_c(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + ! interpolated diffusion coefficient times nabla2(theta) + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO jc = i_startidx, i_endidx +!DIR$ IVDEP +#ifdef _CRAYFTN +!DIR$ PREFERVECTOR +#endif + DO jk = 1, nlev +#else + DO jk = 1, nlev + DO jc = i_startidx, i_endidx +#endif + z_temp(jc,jk,jb) = & + (kh_smag_e(ieidx(jc,jb,1),jk,ieblk(jc,jb,1))*diff_inst(jg)%e_bln_c_s(jc,1,jb) + & + kh_smag_e(ieidx(jc,jb,2),jk,ieblk(jc,jb,2))*diff_inst(jg)%e_bln_c_s(jc,2,jb) + & + kh_smag_e(ieidx(jc,jb,3),jk,ieblk(jc,jb,3))*diff_inst(jg)%e_bln_c_s(jc,3,jb)) * & + (theta_v(jc,jk,jb) *diff_inst(jg)%geofac_n2s(jc,1,jb) + & + theta_v(icidx(jc,jb,1),jk,icblk(jc,jb,1))*diff_inst(jg)%geofac_n2s(jc,2,jb) + & + theta_v(icidx(jc,jb,2),jk,icblk(jc,jb,2))*diff_inst(jg)%geofac_n2s(jc,3,jb) + & + theta_v(icidx(jc,jb,3),jk,icblk(jc,jb,3))*diff_inst(jg)%geofac_n2s(jc,4,jb)) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + ELSE IF (diff_inst(jg)%discr_t == 2) THEN ! use conservative discretization div(k*grad(theta)) + + rl_start = grf_bdywidth_e + rl_end = min_rledge_int - 1 + + i_startblk = p_patch_diff(jg)%edges%start_block(rl_start) + i_endblk = p_patch_diff(jg)%edges%end_block(rl_end) + +!$OMP DO PRIVATE(jk,je,jb,i_startidx,i_endidx), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + ! compute kh_smag_e * grad(theta) (stored in z_nabla2_e for memory efficiency) + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO je = i_startidx, i_endidx +!DIR$ IVDEP +#ifdef _CRAYFTN +!DIR$ PREFERVECTOR +#endif + DO jk = 1, nlev +#else + DO jk = 1, nlev + DO je = i_startidx, i_endidx +#endif + z_nabla2_e(je,jk,jb) = kh_smag_e(je,jk,jb) * & + p_patch_diff(jg)%edges%inv_dual_edge_length(je,jb)* & + (theta_v(iecidx(je,jb,2),jk,iecblk(je,jb,2)) - & + theta_v(iecidx(je,jb,1),jk,iecblk(je,jb,1))) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + rl_start = grf_bdywidth_c+1 + rl_end = min_rlcell_int + + i_startblk = p_patch_diff(jg)%cells%start_block(rl_start) + i_endblk = p_patch_diff(jg)%cells%end_block(rl_end) + +!$OMP DO PRIVATE(jk,jc,jb,i_startidx,i_endidx), ICON_OMP_RUNTIME_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_c(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + ! now compute the divergence of the quantity above + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) +#ifdef __LOOP_EXCHANGE + DO jc = i_startidx, i_endidx + DO jk = 1, nlev +#else + DO jk = 1, nlev + DO jc = i_startidx, i_endidx +#endif + z_temp(jc,jk,jb) = & + z_nabla2_e(ieidx(jc,jb,1),jk,ieblk(jc,jb,1))*diff_inst(jg)%geofac_div(jc,1,jb) + & + z_nabla2_e(ieidx(jc,jb,2),jk,ieblk(jc,jb,2))*diff_inst(jg)%geofac_div(jc,2,jb) + & + z_nabla2_e(ieidx(jc,jb,3),jk,ieblk(jc,jb,3))*diff_inst(jg)%geofac_div(jc,3,jb) + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + ENDIF + + IF (diff_inst(jg)%l_zdiffu_t) THEN ! Compute temperature diffusion truly horizontally over steep slopes + ! A conservative discretization is not possible here +!$OMP DO PRIVATE(jb,jc,ic,nlen_zdiffu,ishift) ICON_OMP_DEFAULT_SCHEDULE + DO jb = 1, nblks_zdiffu + IF (jb == nblks_zdiffu) THEN + nlen_zdiffu = npromz_zdiffu + ELSE + nlen_zdiffu = nproma_zdiffu + ENDIF + ishift = (jb-1)*nproma_zdiffu + !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRESENT(icell, ilev, iblk, vcoef, zd_geofac) & + !$ACC GANG VECTOR ASYNC(1) IF(i_am_accel_node) +!$NEC ivdep +!DIR$ IVDEP + DO jc = 1, nlen_zdiffu + ic = ishift+jc + z_temp(icell(1,ic),ilev(1,ic),iblk(1,ic)) = & + z_temp(icell(1,ic),ilev(1,ic),iblk(1,ic)) + diff_inst(jg)%zd_diffcoef(ic)* & +! MAX(diff_inst(jg)%zd_diffcoef(ic), & +! kh_smag_e(iedge(1,ic),ilev(1,ic),iedblk(1,ic))* blcoef(1,ic) + & +! kh_smag_e(iedge(2,ic),ilev(1,ic),iedblk(2,ic))* blcoef(2,ic) + & +! kh_smag_e(iedge(3,ic),ilev(1,ic),iedblk(3,ic))* blcoef(3,ic) ) * & + (zd_geofac(1,ic)*theta_v(icell(1,ic),ilev(1,ic),iblk(1,ic)) + & + zd_geofac(2,ic)*(vcoef(1,ic)*theta_v(icell(2,ic),ilev(2,ic),iblk(2,ic))+& + (1._wp-vcoef(1,ic))* theta_v(icell(2,ic),ilev(2,ic)+1,iblk(2,ic))) + & + zd_geofac(3,ic)*(vcoef(2,ic)*theta_v(icell(3,ic),ilev(3,ic),iblk(3,ic))+& + (1._wp-vcoef(2,ic))*theta_v(icell(3,ic),ilev(3,ic)+1,iblk(3,ic))) + & + zd_geofac(4,ic)*(vcoef(3,ic)*theta_v(icell(4,ic),ilev(4,ic),iblk(4,ic))+& + (1._wp-vcoef(3,ic))* theta_v(icell(4,ic),ilev(4,ic)+1,iblk(4,ic))) ) + ENDDO + !$ACC END PARALLEL LOOP + ENDDO +!$OMP END DO + + ENDIF + +!$OMP DO PRIVATE(jk,jc,jb,i_startidx,i_endidx,z_theta) ICON_OMP_DEFAULT_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_c(p_patch_diff(jg), jb, i_startblk, i_endblk, & + i_startidx, i_endidx, rl_start, rl_end) + + !$ACC PARALLEL LOOP DEFAULT(PRESENT) GANG VECTOR COLLAPSE(2) ASYNC(1) IF(i_am_accel_node) + DO jk = 1, nlev +!DIR$ IVDEP + DO jc = i_startidx, i_endidx + z_theta = theta_v(jc,jk,jb) + + theta_v(jc,jk,jb) = theta_v(jc,jk,jb) + & + p_patch_diff(jg)%cells%area(jc,jb)*z_temp(jc,jk,jb) + + exner(jc,jk,jb) = exner(jc,jk,jb) * & + (1._wp+rd_o_cvd*(theta_v(jc,jk,jb)/z_theta-1._wp)) + + ENDDO + ENDDO + !$ACC END PARALLEL LOOP + + ENDDO +!$OMP END DO NOWAIT +!$OMP END PARALLEL + + ! This could be further optimized, but applications without physics are quite rare; + IF ( .NOT. diff_inst(jg)%lhdiff_rcf .OR. linit .OR. .NOT. diff_inst(jg)%lphys ) THEN + CALL sync_patch_array_mult(SYNC_C,p_patch(jg),2,theta_v,exner, & + opt_varname="diffusion: theta and exner") + ENDIF + + ENDIF ! temperature diffusion + + IF ( .NOT. diff_inst(jg)%lhdiff_rcf .OR. linit .OR. .NOT. diff_inst(jg)%lphys ) THEN + IF (diff_inst(jg)%lhdiff_w) THEN + CALL sync_patch_array(SYNC_C,p_patch(jg),w,"diffusion: w") + END IF + ENDIF + + IF (diff_inst(jg)%ltimer) CALL timer_stop(timer_nh_hdiffusion) + + !$ACC END DATA + + !$ACC WAIT + + END SUBROUTINE diffusion_run + + !> + !! finalize_diffusion + !! + !! Prepares the horizontal diffusion of velocity and temperature + !! + !! @par Revision History + !! Initial release by William Sawyer, CSCS (2022-11-25) + !! + SUBROUTINE diffusion_finalize(jg) + INTEGER, INTENT(IN) :: jg + ! Currently nothing to do here + END SUBROUTINE diffusion_finalize + +END MODULE mo_nh_diffusion_new diff --git a/testutils/src/icon4py/testutils/fortran/expected_diffusion_granule_savepoint.f90 b/testutils/src/icon4py/testutils/fortran/expected_diffusion_granule_savepoint.f90 new file mode 100644 index 0000000000..c7b8173854 --- /dev/null +++ b/testutils/src/icon4py/testutils/fortran/expected_diffusion_granule_savepoint.f90 @@ -0,0 +1,664 @@ + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_primal_normal_vert_v1(:,:,:) + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_primal_normal_vert_v2(:,:,:) + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_dual_normal_vert_v1(:,:,:) + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_dual_normal_vert_v2(:,:,:) + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_primal_normal_cell_v1(:,:,:) + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_primal_normal_cell_v2(:,:,:) + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_dual_normal_cell_v1(:,:,:) + + !$ser verbatim real, dimension(:,:,:), allocatable :: edges_dual_normal_cell_v2(:,:,:) + + !$ser init directory="." prefix="f2ser" + + !$ser savepoint diffusion_init_in + + !$ser verbatim allocate(edges_primal_normal_vert_v1(size(edges_primal_normal_vert, 1),size(edges_primal_normal_vert, 2),size(edges_primal_normal_vert, 3))) + !$ser verbatim edges_primal_normal_vert_v1 = edges_primal_normal_vert%v1 + !$ser data edges_primal_normal_vert_v1=edges_primal_normal_vert_v1 + !$ser verbatim deallocate(edges_primal_normal_vert_v1) + + !$ser verbatim allocate(edges_primal_normal_vert_v2(size(edges_primal_normal_vert, 1),size(edges_primal_normal_vert, 2),size(edges_primal_normal_vert, 3))) + !$ser verbatim edges_primal_normal_vert_v2 = edges_primal_normal_vert%v2 + !$ser data edges_primal_normal_vert_v2=edges_primal_normal_vert_v2 + !$ser verbatim deallocate(edges_primal_normal_vert_v2) + + !$ser verbatim allocate(edges_dual_normal_vert_v1(size(edges_dual_normal_vert, 1),size(edges_dual_normal_vert, 2),size(edges_dual_normal_vert, 3))) + !$ser verbatim edges_dual_normal_vert_v1 = edges_dual_normal_vert%v1 + !$ser data edges_dual_normal_vert_v1=edges_dual_normal_vert_v1 + !$ser verbatim deallocate(edges_dual_normal_vert_v1) + + !$ser verbatim allocate(edges_dual_normal_vert_v2(size(edges_dual_normal_vert, 1),size(edges_dual_normal_vert, 2),size(edges_dual_normal_vert, 3))) + !$ser verbatim edges_dual_normal_vert_v2 = edges_dual_normal_vert%v2 + !$ser data edges_dual_normal_vert_v2=edges_dual_normal_vert_v2 + !$ser verbatim deallocate(edges_dual_normal_vert_v2) + + !$ser verbatim allocate(edges_primal_normal_cell_v1(size(edges_primal_normal_cell, 1),size(edges_primal_normal_cell, 2),size(edges_primal_normal_cell, 3))) + !$ser verbatim edges_primal_normal_cell_v1 = edges_primal_normal_cell%v1 + !$ser data edges_primal_normal_cell_v1=edges_primal_normal_cell_v1 + !$ser verbatim deallocate(edges_primal_normal_cell_v1) + + !$ser verbatim allocate(edges_primal_normal_cell_v2(size(edges_primal_normal_cell, 1),size(edges_primal_normal_cell, 2),size(edges_primal_normal_cell, 3))) + !$ser verbatim edges_primal_normal_cell_v2 = edges_primal_normal_cell%v2 + !$ser data edges_primal_normal_cell_v2=edges_primal_normal_cell_v2 + !$ser verbatim deallocate(edges_primal_normal_cell_v2) + + !$ser verbatim allocate(edges_dual_normal_cell_v1(size(edges_dual_normal_cell, 1),size(edges_dual_normal_cell, 2),size(edges_dual_normal_cell, 3))) + !$ser verbatim edges_dual_normal_cell_v1 = edges_dual_normal_cell%v1 + !$ser data edges_dual_normal_cell_v1=edges_dual_normal_cell_v1 + !$ser verbatim deallocate(edges_dual_normal_cell_v1) + + !$ser verbatim allocate(edges_dual_normal_cell_v2(size(edges_dual_normal_cell, 1),size(edges_dual_normal_cell, 2),size(edges_dual_normal_cell, 3))) + !$ser verbatim edges_dual_normal_cell_v2 = edges_dual_normal_cell%v2 + !$ser data edges_dual_normal_cell_v2=edges_dual_normal_cell_v2 + !$ser verbatim deallocate(edges_dual_normal_cell_v2) + + PRINT *, 'Serializing cvd_o_rd=cvd_o_rd' + + !$ser data cvd_o_rd=cvd_o_rd + + PRINT *, 'Serializing grav=grav' + + !$ser data grav=grav + + PRINT *, 'Serializing jg=jg' + + !$ser data jg=jg + + PRINT *, 'Serializing nproma=nproma' + + !$ser data nproma=nproma + + PRINT *, 'Serializing nlev=nlev' + + !$ser data nlev=nlev + + PRINT *, 'Serializing nblks_e=nblks_e' + + !$ser data nblks_e=nblks_e + + PRINT *, 'Serializing nblks_v=nblks_v' + + !$ser data nblks_v=nblks_v + + PRINT *, 'Serializing nblks_c=nblks_c' + + !$ser data nblks_c=nblks_c + + PRINT *, 'Serializing nshift=nshift' + + !$ser data nshift=nshift + + PRINT *, 'Serializing nshift_total=nshift_total' + + !$ser data nshift_total=nshift_total + + PRINT *, 'Serializing nrdmax=nrdmax' + + !$ser data nrdmax=nrdmax + + PRINT *, 'Serializing ndyn_substeps=ndyn_substeps' + + !$ser data ndyn_substeps=ndyn_substeps + + PRINT *, 'Serializing hdiff_order=hdiff_order' + + !$ser data hdiff_order=hdiff_order + + PRINT *, 'Serializing itype_comm=itype_comm' + + !$ser data itype_comm=itype_comm + + PRINT *, 'Serializing itype_sher=itype_sher' + + !$ser data itype_sher=itype_sher + + PRINT *, 'Serializing itype_vn_diffu=itype_vn_diffu' + + !$ser data itype_vn_diffu=itype_vn_diffu + + PRINT *, 'Serializing itype_t_diffu=itype_t_diffu' + + !$ser data itype_t_diffu=itype_t_diffu + + PRINT *, 'Serializing hdiff_smag_z=hdiff_smag_z' + + !$ser data hdiff_smag_z=hdiff_smag_z + + PRINT *, 'Serializing hdiff_smag_z2=hdiff_smag_z2' + + !$ser data hdiff_smag_z2=hdiff_smag_z2 + + PRINT *, 'Serializing hdiff_smag_z3=hdiff_smag_z3' + + !$ser data hdiff_smag_z3=hdiff_smag_z3 + + PRINT *, 'Serializing hdiff_smag_z4=hdiff_smag_z4' + + !$ser data hdiff_smag_z4=hdiff_smag_z4 + + PRINT *, 'Serializing hdiff_smag_fac=hdiff_smag_fac' + + !$ser data hdiff_smag_fac=hdiff_smag_fac + + PRINT *, 'Serializing hdiff_smag_fac2=hdiff_smag_fac2' + + !$ser data hdiff_smag_fac2=hdiff_smag_fac2 + + PRINT *, 'Serializing hdiff_smag_fac3=hdiff_smag_fac3' + + !$ser data hdiff_smag_fac3=hdiff_smag_fac3 + + PRINT *, 'Serializing hdiff_smag_fac4=hdiff_smag_fac4' + + !$ser data hdiff_smag_fac4=hdiff_smag_fac4 + + PRINT *, 'Serializing hdiff_efdt_ratio=hdiff_efdt_ratio' + + !$ser data hdiff_efdt_ratio=hdiff_efdt_ratio + + PRINT *, 'Serializing k4=k4' + + !$ser data k4=k4 + + PRINT *, 'Serializing k4w=k4w' + + !$ser data k4w=k4w + + PRINT *, 'Serializing nudge_max_coeff=nudge_max_coeff' + + !$ser data nudge_max_coeff=nudge_max_coeff + + PRINT *, 'Serializing denom_diffu_v=denom_diffu_v' + + !$ser data denom_diffu_v=denom_diffu_v + + PRINT *, 'Serializing p_test_run=p_test_run' + + !$ser data p_test_run=p_test_run + + PRINT *, 'Serializing lphys=lphys' + + !$ser data lphys=lphys + + PRINT *, 'Serializing lhdiff_rcf=lhdiff_rcf' + + !$ser data lhdiff_rcf=lhdiff_rcf + + PRINT *, 'Serializing lhdiff_w=lhdiff_w' + + !$ser data lhdiff_w=lhdiff_w + + PRINT *, 'Serializing lhdiff_temp=lhdiff_temp' + + !$ser data lhdiff_temp=lhdiff_temp + + PRINT *, 'Serializing l_zdiffu_t=l_zdiffu_t' + + !$ser data l_zdiffu_t=l_zdiffu_t + + PRINT *, 'Serializing l_limited_area=l_limited_area' + + !$ser data l_limited_area=l_limited_area + + PRINT *, 'Serializing lfeedback=lfeedback' + + !$ser data lfeedback=lfeedback + + PRINT *, 'Serializing ltkeshs=ltkeshs' + + !$ser data ltkeshs=ltkeshs + + PRINT *, 'Serializing lsmag_3d=lsmag_3d' + + !$ser data lsmag_3d=lsmag_3d + + PRINT *, 'Serializing lvert_nest=lvert_nest' + + !$ser data lvert_nest=lvert_nest + + PRINT *, 'Serializing ltimer=ltimer' + + !$ser data ltimer=ltimer + + PRINT *, 'Serializing ddt_vn_hdf_is_associated=ddt_vn_hdf_is_associated' + + !$ser data ddt_vn_hdf_is_associated=ddt_vn_hdf_is_associated + + PRINT *, 'Serializing ddt_vn_dyn_is_associated=ddt_vn_dyn_is_associated' + + !$ser data ddt_vn_dyn_is_associated=ddt_vn_dyn_is_associated + + PRINT *, 'Serializing vct_a=vct_a(:)' + + IF (SIZE(vct_a) > 0) THEN + !$ser data vct_a=vct_a(:) + ELSE + PRINT *, 'Warning: Array vct_a has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing c_lin_e=c_lin_e(:,:,:)' + + IF (SIZE(c_lin_e) > 0) THEN + !$ser data c_lin_e=c_lin_e(:,:,:) + ELSE + PRINT *, 'Warning: Array c_lin_e has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing e_bln_c_s=e_bln_c_s(:,:,:)' + + IF (SIZE(e_bln_c_s) > 0) THEN + !$ser data e_bln_c_s=e_bln_c_s(:,:,:) + ELSE + PRINT *, 'Warning: Array e_bln_c_s has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing e_bln_c_u=e_bln_c_u(:,:,:)' + + IF (SIZE(e_bln_c_u) > 0) THEN + !$ser data e_bln_c_u=e_bln_c_u(:,:,:) + ELSE + PRINT *, 'Warning: Array e_bln_c_u has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing e_bln_c_v=e_bln_c_v(:,:,:)' + + IF (SIZE(e_bln_c_v) > 0) THEN + !$ser data e_bln_c_v=e_bln_c_v(:,:,:) + ELSE + PRINT *, 'Warning: Array e_bln_c_v has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_aw_verts=cells_aw_verts(:,:,:)' + + IF (SIZE(cells_aw_verts) > 0) THEN + !$ser data cells_aw_verts=cells_aw_verts(:,:,:) + ELSE + PRINT *, 'Warning: Array cells_aw_verts has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing geofac_div=geofac_div(:,:,:)' + + IF (SIZE(geofac_div) > 0) THEN + !$ser data geofac_div=geofac_div(:,:,:) + ELSE + PRINT *, 'Warning: Array geofac_div has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing geofac_rot=geofac_rot(:,:,:)' + + IF (SIZE(geofac_rot) > 0) THEN + !$ser data geofac_rot=geofac_rot(:,:,:) + ELSE + PRINT *, 'Warning: Array geofac_rot has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing geofac_n2s=geofac_n2s(:,:,:)' + + IF (SIZE(geofac_n2s) > 0) THEN + !$ser data geofac_n2s=geofac_n2s(:,:,:) + ELSE + PRINT *, 'Warning: Array geofac_n2s has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing geofac_grg=geofac_grg(:,:,:,:)' + + IF (SIZE(geofac_grg) > 0) THEN + !$ser data geofac_grg=geofac_grg(:,:,:,:) + ELSE + PRINT *, 'Warning: Array geofac_grg has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing nudgecoeff_e=nudgecoeff_e(:,:)' + + IF (SIZE(nudgecoeff_e) > 0) THEN + !$ser data nudgecoeff_e=nudgecoeff_e(:,:) + ELSE + PRINT *, 'Warning: Array nudgecoeff_e has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing rbf_vec_idx_v=rbf_vec_idx_v(:,:,:)' + + IF (SIZE(rbf_vec_idx_v) > 0) THEN + !$ser data rbf_vec_idx_v=rbf_vec_idx_v(:,:,:) + ELSE + PRINT *, 'Warning: Array rbf_vec_idx_v has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing rbf_vec_blk_v=rbf_vec_blk_v(:,:,:)' + + IF (SIZE(rbf_vec_blk_v) > 0) THEN + !$ser data rbf_vec_blk_v=rbf_vec_blk_v(:,:,:) + ELSE + PRINT *, 'Warning: Array rbf_vec_blk_v has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing rbf_vec_coeff_v=rbf_vec_coeff_v(:,:,:,:)' + + IF (SIZE(rbf_vec_coeff_v) > 0) THEN + !$ser data rbf_vec_coeff_v=rbf_vec_coeff_v(:,:,:,:) + ELSE + PRINT *, 'Warning: Array rbf_vec_coeff_v has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing enhfac_diffu=enhfac_diffu(:)' + + IF (SIZE(enhfac_diffu) > 0) THEN + !$ser data enhfac_diffu=enhfac_diffu(:) + ELSE + PRINT *, 'Warning: Array enhfac_diffu has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing zd_intcoef=zd_intcoef(:,:)' + + IF (SIZE(zd_intcoef) > 0) THEN + !$ser data zd_intcoef=zd_intcoef(:,:) + ELSE + PRINT *, 'Warning: Array zd_intcoef has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing zd_geofac=zd_geofac(:,:)' + + IF (SIZE(zd_geofac) > 0) THEN + !$ser data zd_geofac=zd_geofac(:,:) + ELSE + PRINT *, 'Warning: Array zd_geofac has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing zd_diffcoef=zd_diffcoef(:)' + + IF (SIZE(zd_diffcoef) > 0) THEN + !$ser data zd_diffcoef=zd_diffcoef(:) + ELSE + PRINT *, 'Warning: Array zd_diffcoef has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing wgtfac_c=wgtfac_c(:,:,:)' + + IF (SIZE(wgtfac_c) > 0) THEN + !$ser data wgtfac_c=wgtfac_c(:,:,:) + ELSE + PRINT *, 'Warning: Array wgtfac_c has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing wgtfac_e=wgtfac_e(:,:,:)' + + IF (SIZE(wgtfac_e) > 0) THEN + !$ser data wgtfac_e=wgtfac_e(:,:,:) + ELSE + PRINT *, 'Warning: Array wgtfac_e has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing wgtfacq_e=wgtfacq_e(:,:,:)' + + IF (SIZE(wgtfacq_e) > 0) THEN + !$ser data wgtfacq_e=wgtfacq_e(:,:,:) + ELSE + PRINT *, 'Warning: Array wgtfacq_e has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing wgtfacq1_e=wgtfacq1_e(:,:,:)' + + IF (SIZE(wgtfacq1_e) > 0) THEN + !$ser data wgtfacq1_e=wgtfacq1_e(:,:,:) + ELSE + PRINT *, 'Warning: Array wgtfacq1_e has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing ddqz_z_full_e=ddqz_z_full_e(:,:,:)' + + IF (SIZE(ddqz_z_full_e) > 0) THEN + !$ser data ddqz_z_full_e=ddqz_z_full_e(:,:,:) + ELSE + PRINT *, 'Warning: Array ddqz_z_full_e has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing theta_ref_mc=theta_ref_mc(:,:,:)' + + IF (SIZE(theta_ref_mc) > 0) THEN + !$ser data theta_ref_mc=theta_ref_mc(:,:,:) + ELSE + PRINT *, 'Warning: Array theta_ref_mc has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing zd_indlist=zd_indlist(:,:)' + + IF (SIZE(zd_indlist) > 0) THEN + !$ser data zd_indlist=zd_indlist(:,:) + ELSE + PRINT *, 'Warning: Array zd_indlist has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing zd_blklist=zd_blklist(:,:)' + + IF (SIZE(zd_blklist) > 0) THEN + !$ser data zd_blklist=zd_blklist(:,:) + ELSE + PRINT *, 'Warning: Array zd_blklist has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing zd_vertidx=zd_vertidx(:,:)' + + IF (SIZE(zd_vertidx) > 0) THEN + !$ser data zd_vertidx=zd_vertidx(:,:) + ELSE + PRINT *, 'Warning: Array zd_vertidx has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing zd_listdim=zd_listdim' + + !$ser data zd_listdim=zd_listdim + + PRINT *, 'Serializing edges_start_block=edges_start_block(min_rledge:)' + + IF (SIZE(edges_start_block) > 0) THEN + !$ser data edges_start_block=edges_start_block(min_rledge:) + ELSE + PRINT *, 'Warning: Array edges_start_block has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_end_block=edges_end_block(min_rledge:)' + + IF (SIZE(edges_end_block) > 0) THEN + !$ser data edges_end_block=edges_end_block(min_rledge:) + ELSE + PRINT *, 'Warning: Array edges_end_block has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_start_index=edges_start_index(min_rledge:)' + + IF (SIZE(edges_start_index) > 0) THEN + !$ser data edges_start_index=edges_start_index(min_rledge:) + ELSE + PRINT *, 'Warning: Array edges_start_index has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_end_index=edges_end_index(min_rledge:)' + + IF (SIZE(edges_end_index) > 0) THEN + !$ser data edges_end_index=edges_end_index(min_rledge:) + ELSE + PRINT *, 'Warning: Array edges_end_index has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_vertex_idx=edges_vertex_idx(:,:,:)' + + IF (SIZE(edges_vertex_idx) > 0) THEN + !$ser data edges_vertex_idx=edges_vertex_idx(:,:,:) + ELSE + PRINT *, 'Warning: Array edges_vertex_idx has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_vertex_blk=edges_vertex_blk(:,:,:)' + + IF (SIZE(edges_vertex_blk) > 0) THEN + !$ser data edges_vertex_blk=edges_vertex_blk(:,:,:) + ELSE + PRINT *, 'Warning: Array edges_vertex_blk has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_cell_idx=edges_cell_idx(:,:,:)' + + IF (SIZE(edges_cell_idx) > 0) THEN + !$ser data edges_cell_idx=edges_cell_idx(:,:,:) + ELSE + PRINT *, 'Warning: Array edges_cell_idx has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_cell_blk=edges_cell_blk(:,:,:)' + + IF (SIZE(edges_cell_blk) > 0) THEN + !$ser data edges_cell_blk=edges_cell_blk(:,:,:) + ELSE + PRINT *, 'Warning: Array edges_cell_blk has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_tangent_orientation=edges_tangent_orientation(:,:)' + + IF (SIZE(edges_tangent_orientation) > 0) THEN + !$ser data edges_tangent_orientation=edges_tangent_orientation(:,:) + ELSE + PRINT *, 'Warning: Array edges_tangent_orientation has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_inv_vert_vert_length=edges_inv_vert_vert_length(:,:)' + + IF (SIZE(edges_inv_vert_vert_length) > 0) THEN + !$ser data edges_inv_vert_vert_length=edges_inv_vert_vert_length(:,:) + ELSE + PRINT *, 'Warning: Array edges_inv_vert_vert_length has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_inv_primal_edge_length=edges_inv_primal_edge_length(:,:)' + + IF (SIZE(edges_inv_primal_edge_length) > 0) THEN + !$ser data edges_inv_primal_edge_length=edges_inv_primal_edge_length(:,:) + ELSE + PRINT *, 'Warning: Array edges_inv_primal_edge_length has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_inv_dual_edge_length=edges_inv_dual_edge_length(:,:)' + + IF (SIZE(edges_inv_dual_edge_length) > 0) THEN + !$ser data edges_inv_dual_edge_length=edges_inv_dual_edge_length(:,:) + ELSE + PRINT *, 'Warning: Array edges_inv_dual_edge_length has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing edges_area_edge=edges_area_edge(:,:)' + + IF (SIZE(edges_area_edge) > 0) THEN + !$ser data edges_area_edge=edges_area_edge(:,:) + ELSE + PRINT *, 'Warning: Array edges_area_edge has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_start_block=cells_start_block(min_rlcell:)' + + IF (SIZE(cells_start_block) > 0) THEN + !$ser data cells_start_block=cells_start_block(min_rlcell:) + ELSE + PRINT *, 'Warning: Array cells_start_block has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_end_block=cells_end_block(min_rlcell:)' + + IF (SIZE(cells_end_block) > 0) THEN + !$ser data cells_end_block=cells_end_block(min_rlcell:) + ELSE + PRINT *, 'Warning: Array cells_end_block has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_start_index=cells_start_index(min_rlcell:)' + + IF (SIZE(cells_start_index) > 0) THEN + !$ser data cells_start_index=cells_start_index(min_rlcell:) + ELSE + PRINT *, 'Warning: Array cells_start_index has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_end_index=cells_end_index(min_rlcell:)' + + IF (SIZE(cells_end_index) > 0) THEN + !$ser data cells_end_index=cells_end_index(min_rlcell:) + ELSE + PRINT *, 'Warning: Array cells_end_index has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_neighbor_idx=cells_neighbor_idx(:,:,:)' + + IF (SIZE(cells_neighbor_idx) > 0) THEN + !$ser data cells_neighbor_idx=cells_neighbor_idx(:,:,:) + ELSE + PRINT *, 'Warning: Array cells_neighbor_idx has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_neighbor_blk=cells_neighbor_blk(:,:,:)' + + IF (SIZE(cells_neighbor_blk) > 0) THEN + !$ser data cells_neighbor_blk=cells_neighbor_blk(:,:,:) + ELSE + PRINT *, 'Warning: Array cells_neighbor_blk has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_edge_idx=cells_edge_idx(:,:,:)' + + IF (SIZE(cells_edge_idx) > 0) THEN + !$ser data cells_edge_idx=cells_edge_idx(:,:,:) + ELSE + PRINT *, 'Warning: Array cells_edge_idx has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_edge_blk=cells_edge_blk(:,:,:)' + + IF (SIZE(cells_edge_blk) > 0) THEN + !$ser data cells_edge_blk=cells_edge_blk(:,:,:) + ELSE + PRINT *, 'Warning: Array cells_edge_blk has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing cells_area=cells_area(:,:)' + + IF (SIZE(cells_area) > 0) THEN + !$ser data cells_area=cells_area(:,:) + ELSE + PRINT *, 'Warning: Array cells_area has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing verts_start_block=verts_start_block(min_rlvert:)' + + IF (SIZE(verts_start_block) > 0) THEN + !$ser data verts_start_block=verts_start_block(min_rlvert:) + ELSE + PRINT *, 'Warning: Array verts_start_block has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing verts_end_block=verts_end_block(min_rlvert:)' + + IF (SIZE(verts_end_block) > 0) THEN + !$ser data verts_end_block=verts_end_block(min_rlvert:) + ELSE + PRINT *, 'Warning: Array verts_end_block has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing verts_start_index=verts_start_index(min_rlvert:)' + + IF (SIZE(verts_start_index) > 0) THEN + !$ser data verts_start_index=verts_start_index(min_rlvert:) + ELSE + PRINT *, 'Warning: Array verts_start_index has size 0. Not serializing array.' + END IF + + PRINT *, 'Serializing verts_end_index=verts_end_index(min_rlvert:)' + + IF (SIZE(verts_end_index) > 0) THEN + !$ser data verts_end_index=verts_end_index(min_rlvert:) + ELSE + PRINT *, 'Warning: Array verts_end_index has size 0. Not serializing array.' + END IF diff --git a/testutils/src/icon4py/testutils/fortran/no_deps_subroutine_example.f90 b/testutils/src/icon4py/testutils/fortran/no_deps_subroutine_example.f90 new file mode 100644 index 0000000000..021cb4aa76 --- /dev/null +++ b/testutils/src/icon4py/testutils/fortran/no_deps_subroutine_example.f90 @@ -0,0 +1,25 @@ +MODULE no_deps_example_subroutines + IMPLICIT NONE + + PUBLIC :: no_deps_init, no_deps_run + PRIVATE + + CONTAINS + + SUBROUTINE no_deps_init(a, b, c) + REAL, INTENT(in) :: a + REAL, INTENT(inout) :: b + REAL, INTENT(out) :: c + c = a + b + b = 2.0 * b + END SUBROUTINE no_deps_init + + SUBROUTINE no_deps_run(a, b, c) + REAL, INTENT(in) :: a + REAL, INTENT(inout) :: b + REAL, INTENT(out) :: c + c = a + b + b = 2.0 * b + END SUBROUTINE no_deps_run + +END MODULE no_deps_example_subroutines diff --git a/testutils/src/icon4py/testutils/fortran/subroutine_example.f90 b/testutils/src/icon4py/testutils/fortran/subroutine_example.f90 new file mode 100644 index 0000000000..bb22e8c041 --- /dev/null +++ b/testutils/src/icon4py/testutils/fortran/subroutine_example.f90 @@ -0,0 +1,29 @@ +MODULE example_subroutines + USE ISO_C_BINDING, ONLY: C_DOUBLE + IMPLICIT NONE + + PUBLIC :: mysubroutine, foo_type + PRIVATE + + TYPE, BIND(C) :: foo_type + REAL(C_DOUBLE) :: p1 + REAL(C_DOUBLE) :: p2 + END TYPE foo_type + + CONTAINS + + SUBROUTINE mysubroutine(a, b, c) + REAL, INTENT(in) :: a + REAL, INTENT(inout) :: b + REAL, INTENT(out) :: c + c = a + b + b = 2.0 * b + end SUBROUTINE mysubroutine + +END MODULE example_subroutines + + + + + + diff --git a/liskov/tests/samples/fortran_samples.py b/testutils/src/icon4py/testutils/liskov_fortran_samples.py similarity index 74% rename from liskov/tests/samples/fortran_samples.py rename to testutils/src/icon4py/testutils/liskov_fortran_samples.py index 858b3a2089..8231dbe038 100644 --- a/liskov/tests/samples/fortran_samples.py +++ b/testutils/src/icon4py/testutils/liskov_fortran_samples.py @@ -33,11 +33,11 @@ !$DSL DECLARE(vn=nproma,p_patch%nlev,p_patch%nblks_e; suffix=dsl) - !$DSL DECLARE(vn=nproma,p_patch%nlev,p_patch%nblks_e; a=nproma,p_patch%nlev,p_patch%nblks_e; & + !$DSL DECLARE(vn= nproma,p_patch%nlev,p_patch%nblks_e; a=nproma,p_patch%nlev,p_patch%nblks_e; & !$DSL b=nproma,p_patch%nlev,p_patch%nblks_e; type=REAL(vp)) !$DSL START STENCIL(name=apply_nabla2_to_vn_in_lateral_boundary; & - !$DSL z_nabla2_e=z_nabla2_e(:,:,1); area_edge=p_patch%edges%area_edge(:,1); & + !$DSL z_nabla2_e=z_nabla2_e(:, :, 1); area_edge=p_patch%edges%area_edge(:,1); & !$DSL fac_bdydiff_v=fac_bdydiff_v; vn=p_nh_prog%vn(:,:,1); & !$DSL vertical_lower=1; vertical_upper=nlev; & !$DSL horizontal_lower=i_startidx; horizontal_upper=i_endidx; & @@ -242,3 +242,78 @@ !$DSL END CREATE() """ + + +REPEATED_STENCILS = """\ + !$DSL IMPORTS() + + !$DSL START CREATE() + + !$DSL DECLARE(vn=nproma,p_patch%nlev,p_patch%nblks_e; suffix=dsl) + + !$DSL DECLARE(vn= nproma,p_patch%nlev,p_patch%nblks_e; a=nproma,p_patch%nlev,p_patch%nblks_e; & + !$DSL b=nproma,p_patch%nlev,p_patch%nblks_e; type=REAL(vp)) + + !$DSL START STENCIL(name=apply_nabla2_to_vn_in_lateral_boundary; & + !$DSL z_nabla2_e=z_nabla2_e(:, :, 1); area_edge=p_patch%edges%area_edge(:,1); & + !$DSL fac_bdydiff_v=fac_bdydiff_v; vn=p_nh_prog%vn(:,:,1); & + !$DSL vertical_lower=1; vertical_upper=nlev; & + !$DSL horizontal_lower=i_startidx; horizontal_upper=i_endidx; & + !$DSL accpresent=True) + !$OMP DO PRIVATE(je,jk,jb,i_startidx,i_endidx) ICON_OMP_DEFAULT_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch, jb, i_startblk, i_endblk, & + i_startidx, i_endidx, start_bdydiff_e, grf_bdywidth_e) + + !$ACC PARALLEL IF( i_am_accel_node .AND. acc_on ) DEFAULT(NONE) ASYNC(1) + vn_before(:,:,:) = p_nh_prog%vn(:,:,:) + !$ACC END PARALLEL + + !$ACC PARALLEL LOOP DEFAULT(NONE) GANG VECTOR COLLAPSE(2) ASYNC(1) IF( i_am_accel_node .AND. acc_on ) + DO jk = 1, nlev + !DIR$ IVDEP + DO je = i_startidx, i_endidx + p_nh_prog%vn(je,jk,jb) = & + p_nh_prog%vn(je,jk,jb) + & + z_nabla2_e(je,jk,jb) * & + p_patch%edges%area_edge(je,jb)*fac_bdydiff_v + ENDDO + ENDDO + !$DSL START PROFILE(name=apply_nabla2_to_vn_in_lateral_boundary) + !$ACC END PARALLEL LOOP + !$DSL END PROFILE() + !$DSL END STENCIL(name=apply_nabla2_to_vn_in_lateral_boundary; noprofile=True) + + !$DSL START STENCIL(name=apply_nabla2_to_vn_in_lateral_boundary; & + !$DSL z_nabla2_e=z_nabla2_e(:, :, 1); area_edge=p_patch%edges%area_edge(:,1); & + !$DSL fac_bdydiff_v=fac_bdydiff_v; vn=p_nh_prog%vn(:,:,1); & + !$DSL vertical_lower=1; vertical_upper=nlev; & + !$DSL horizontal_lower=i_startidx; horizontal_upper=i_endidx; & + !$DSL accpresent=True) + !$OMP DO PRIVATE(je,jk,jb,i_startidx,i_endidx) ICON_OMP_DEFAULT_SCHEDULE + DO jb = i_startblk,i_endblk + + CALL get_indices_e(p_patch, jb, i_startblk, i_endblk, & + i_startidx, i_endidx, start_bdydiff_e, grf_bdywidth_e) + + !$ACC PARALLEL IF( i_am_accel_node .AND. acc_on ) DEFAULT(NONE) ASYNC(1) + vn_before(:,:,:) = p_nh_prog%vn(:,:,:) + !$ACC END PARALLEL + + !$ACC PARALLEL LOOP DEFAULT(NONE) GANG VECTOR COLLAPSE(2) ASYNC(1) IF( i_am_accel_node .AND. acc_on ) + DO jk = 1, nlev + !DIR$ IVDEP + DO je = i_startidx, i_endidx + p_nh_prog%vn(je,jk,jb) = & + p_nh_prog%vn(je,jk,jb) + & + z_nabla2_e(je,jk,jb) * & + p_patch%edges%area_edge(je,jb)*fac_bdydiff_v + ENDDO + ENDDO + !$DSL START PROFILE(name=apply_nabla2_to_vn_in_lateral_boundary) + !$ACC END PARALLEL LOOP + !$DSL END PROFILE() + !$DSL END STENCIL(name=apply_nabla2_to_vn_in_lateral_boundary; noprofile=True) + !$DSL END CREATE() + """ diff --git a/testutils/src/icon4py/testutils/liskov_test_utils.py b/testutils/src/icon4py/testutils/liskov_test_utils.py new file mode 100644 index 0000000000..526dcd39cf --- /dev/null +++ b/testutils/src/icon4py/testutils/liskov_test_utils.py @@ -0,0 +1,29 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from pathlib import Path + +import icon4py.liskov.parsing.types as ts +from icon4py.liskov.parsing.scan import DirectivesScanner + + +def scan_for_directives(fpath: Path) -> list[ts.RawDirective]: + collector = DirectivesScanner(fpath) + return collector() + + +def insert_new_lines(fname: Path, lines: list[str]) -> None: + """Append new lines into file.""" + with open(fname, "a") as f: + for ln in lines: + f.write(f"{ln}\n")