diff --git a/.github/windows_arm64_steps/action.yml b/.github/windows_arm64_steps/action.yml new file mode 100644 index 00000000..65f6530a --- /dev/null +++ b/.github/windows_arm64_steps/action.yml @@ -0,0 +1,28 @@ +name: Build Dependencies(Win-ARM64) +description: "Common setup steps for Win-ARM64 CI" +runs: + using: "composite" + steps: + - name: Install LLVM + shell: pwsh + run: | + Invoke-WebRequest https://github.com/llvm/llvm-project/releases/download/llvmorg-20.1.8/LLVM-20.1.8-woa64.exe -UseBasicParsing -OutFile LLVM-woa64.exe + $expectedHash = "7c4ac97eb2ae6b960ca5f9caf3ff6124c8d2a18cc07a7840a4d2ea15537bad8e" + $fileHash = (Get-FileHash -Path "LLVM-woa64.exe" -Algorithm SHA256).Hash + if ($fileHash -ne $expectedHash) { + Write-Error "Checksum verification failed. The downloaded file may be corrupted or tampered with." + exit 1 + } + Start-Process -FilePath ".\LLVM-woa64.exe" -ArgumentList "/S" -Wait + echo "C:\Program Files\LLVM\bin" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append + + - name: Install pkgconf via vcpkg + shell: pwsh + run: | + & "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" arm64 + $env:VCPKG_ROOT = "C:\vcpkg" + Set-Location $env:VCPKG_ROOT + ./vcpkg install pkgconf:arm64-windows + $pkgconfPath = "$env:VCPKG_ROOT\installed\arm64-windows\tools\pkgconf" + Copy-Item "$pkgconfPath\pkgconf.exe" "$pkgconfPath\pkg-config.exe" -Force + echo "$pkgconfPath" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append \ No newline at end of file diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index 85850f9f..7232b1bd 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -17,6 +17,8 @@ jobs: name: ${{ matrix.python[0] }}-${{ matrix.buildplat[1] }} ${{ matrix.buildplat[2] }} runs-on: ${{ matrix.buildplat[0] }} strategy: + #DEBUGGING: turning off the fail-fast can be helpful + # fail-fast: false # This ensures all matrix jobs run to completion matrix: buildplat: # Not supported for now @@ -29,12 +31,20 @@ jobs: - [windows-2025, win, AMD64] - [ubuntu-22.04, manylinux, x86_64] - [macos-14, macosx, arm64, openblas, "12.3"] - - [macos-15-intel, macosx, x86_64, openblas, "10.14"] + - [windows-11-arm, win, ARM64, "", ""] + # removed + # - [macos-15-intel, macosx, x86_64, openblas, "10.14"] - python: [["cp39", "3.9"], ["cp310", "3.10"], ["cp311", "3.11"], ["cp312", "3.12"], ["cp313", "3.13"]] + # DEBUGGING: switch to just one python for debugging # python: [["cp312", "3.12"]] - + + python: [["cp39", "3.9"], ["cp310", "3.10"], ["cp311", "3.11"], ["cp312", "3.12"], ["cp313", "3.13"], ["cp314", "3.14"]] + exclude: + - buildplat: [windows-11-arm, win, ARM64, "", ""] + python: ["cp310", "3.10"] + - buildplat: [windows-11-arm, win, ARM64, "", ""] + python: ["cp39", "3.9"] steps: @@ -55,6 +65,19 @@ jobs: if: ${{ runner.os == 'Windows' }} + - name: Set environment variables for ARM64 + if: matrix.buildplat[1] == 'win' && matrix.buildplat[2] == 'ARM64' + run: | + echo "CC=clang-cl" >> $env:GITHUB_ENV + echo "CXX=clang-cl" >> $env:GITHUB_ENV + echo "FC=flang" >> $env:GITHUB_ENV + echo "TARGET_ARCH=${{ matrix.buildplat[2] }}" >> $env:GITHUB_ENV + + - name: Set up Flang and pkgconf for ARM64 + if: matrix.buildplat[1] == 'win' && matrix.buildplat[2] == 'ARM64' + uses: ./.github/windows_arm64_steps + + - name: Setup macOS if: startsWith( matrix.buildplat[0], 'macos-' ) run: | diff --git a/.gitignore b/.gitignore index 01894e0f..b19204cd 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,11 @@ dev_make_install/optvl tests/test_om_wrapper*_out examples/run_opt_om*_out +examples/*.avl +examples/*.ps +examples/*.png +examples/*.eig +examples/opt_output_sweep/ # ============================================================================== # Github python gitignore template # ============================================================================== @@ -52,7 +57,6 @@ MANIFEST src/ad_src/forward_tmp/ src/ad_src/reverse_tmp/ src/ad_src/preprocessed_files/ -src/includes/AVL.INC examples/*.html examples/*.db examples/reports diff --git a/README.md b/README.md index b8808744..287cb8d9 100644 --- a/README.md +++ b/README.md @@ -57,7 +57,7 @@ print( f'CL:{force_data["CL"]:10.6f} CD:{force_data["CD"]:10.6f} Cm:{force_data["Cm"]:10.6f}' ) -# lets look at the cp countours +# let's look at the cp contours ovl.plot_cp() ``` The plotting calls in this example use matplotlib to produce visualizations like this diff --git a/meson.build b/meson.build index 89d64109..8259b5dc 100644 --- a/meson.build +++ b/meson.build @@ -2,13 +2,12 @@ project( 'optvl', 'c', # don't forget to change the value in pyproject.toml - version: '2.4.0', + version: '2.5.0', license: 'GPL-3.0', meson_version: '>= 0.64.0', default_options: [ 'buildtype=debugoptimized', 'c_std=c99', - 'fortran_std=legacy', ], ) @@ -86,6 +85,14 @@ inc_np = include_directories(incdir_numpy) incdir_f2py = incdir_numpy / '..' / '..' / 'f2py' / 'src' inc_f2py = include_directories(incdir_f2py) +# set up statically allocated variable sizes +if host_machine.system() == 'windows' and host_machine.cpu_family() == 'aarch64' + run_command( + py3, '-c', + 'import shutil; shutil.copy(r"' + (meson.project_source_root() / 'src' / 'includes' / 'ADIMEN.INC.winArm64') + '", r"' + (meson.project_source_root() / 'src' / 'includes' / 'ADIMEN.INC') + '")', + check: true, + ) +endif # Fortran warning flags @@ -111,13 +118,28 @@ ff_args = ff.get_supported_arguments( ff_args += [ '-ffixed-line-length-80', - '-std=legacy', '-fdefault-real-8', '-fdefault-double-8', - '-fPIC', '-O2' ] +if host_machine.system() != 'windows' + ff_args += ['-fPIC'] +endif + +# --- add flags for fortran standard --- +if ff.get_id() == 'gcc' + # -std=legacy is not supported by all Fortran compilers, but very useful with + # gfortran since it avoids a ton of warnings that we don't care about. + # Needs fixing in Meson, see https://github.com/mesonbuild/meson/issues/11633. + ff_args += ['-std=legacy'] +elif ff.get_id() == 'llvm-flang' + + ff_args += ['-fno-init-global-zero', '-ffixed-form', '-w'] + link_args += ['-fno-init-global-zero'] + +endif + numpy_nodepr_api = '-DNPY_NO_DEPRECATED_API=NPY_1_9_API_VERSION' cc_args = [numpy_nodepr_api] @@ -179,4 +201,5 @@ py3.extension_module('libavl', ) + install_subdir('optvl', install_dir: py3.get_install_dir()) diff --git a/optvl/MExt.py b/optvl/MExt.py index a79a5d97..6b8a8310 100644 --- a/optvl/MExt.py +++ b/optvl/MExt.py @@ -67,11 +67,13 @@ def __init__(self, libName, packageName, pip_name, lib_so_file=None, debug=False if not os.path.exists(target_path) and os.path.exists(source_path): # print("Creating symlink from {} to {}".format(source_path, target_path)) os.symlink(source_path, target_path) + elif platform.system() == "Windows": + libs_dir_name = f"{pip_name}.libs" + source_path = os.path.join(spec.submodule_search_locations[0], "..", libs_dir_name) + if os.path.exists(source_path): + os.add_dll_directory(os.path.abspath(source_path)) else: - # raise NotImplementedError("platform not supported") - print(tmpdir) - srcpath = os.path.join(spec.submodule_search_locations[0], "libavl.cp39-win_amd64.dll.a") - pass + raise RuntimeError("Platform not recognized") # add the directory containing the new package to the search path sys.path.insert(0, tmpdir) @@ -90,15 +92,20 @@ def __init__(self, libName, packageName, pip_name, lib_so_file=None, debug=False self.__dict__.update(self._module.__dict__) def __del__(self): - # remove module if not in debug mode if not self.debug: # if the module was imported, remove it from sys.modules if hasattr(self, "_pkg"): del sys.modules[self._module.__name__] del sys.modules[self._pkg.__name__] - # now try to delete the files and directory - shutil.rmtree(self._pkgdir) + try: + # now try to delete the files and directory + shutil.rmtree(self._pkgdir) + except PermissionError: + # On Windows, loaded DLLs are locked and cannot be deleted. + # Clean up as much as possible, ignore what's locked. + pass + # make sure the original module is loaded - # otherwise python crashes on exit # if MExt objects have not been explicitly 'del'd, diff --git a/optvl/om_wrapper.py b/optvl/om_wrapper.py index 9a433951..cb32757e 100644 --- a/optvl/om_wrapper.py +++ b/optvl/om_wrapper.py @@ -90,7 +90,7 @@ def setup(self): ) -AIRFOIL_GEOM_VARS = ["xasec", "casec", "tasec", "xuasec", "xlasec", "zuasec", "zlasec"] +AIRFOIL_GEOM_VARS = ["xasec", "casec", "tasec"] # helper functions used by the AVL components diff --git a/optvl/optvl_class.py b/optvl/optvl_class.py index 8ab2a536..24d55eb2 100644 --- a/optvl/optvl_class.py +++ b/optvl/optvl_class.py @@ -221,22 +221,6 @@ class OVLSolver(object): ad_suffix = "_DIFF" - # Primary array limits: These also need to updated in the Fortran layer if changed - NVMAX = 5000 # number of horseshoe vortices - NSMAX = 500 # number of chord strips - NSECMAX = 301 # nuber of geometry sections - NFMAX = 100 # number of surfaces - NLMAX = 502 # number of source/doublet line nodes - NBMAX = 20 # number of bodies - NUMAX = 6 # number of freestream parameters (V,Omega) - NDMAX = 30 # number of control deflection parameters - NGMAX = 21 # number of design variables - NRMAX = 25 # number of stored run cases - NTMAX = 503 # number of stored time levels - NOBMAX=1 # max number of off body points - ICONX = 20 # - IBX = 200 # max number of airfoil coordinates - def __init__( self, @@ -275,7 +259,6 @@ def __init__( if platform.system() == "Windows": # HACK avl_lib_so_file = glob.glob(os.path.join(module_dir, "libavl*.pyd"))[0] - print("avl_lib_so_file", avl_lib_so_file) if len(avl_lib_so_file) == 0: # print the contents of the module dir print("module_dir", module_dir) @@ -286,6 +269,7 @@ def __init__( # # get just the file name avl_lib_so_file = os.path.basename(avl_lib_so_file) + self.avl = MExt.MExt("libavl", module_name, "optvl", lib_so_file=avl_lib_so_file, debug=debug)._module @@ -296,6 +280,8 @@ def __init__( if timing: self.set_avl_fort_arr("CASE_L", "LTIMING", True) + + self.__set_avl_size_info() if geo_file is not None: try: @@ -418,6 +404,27 @@ def __init__( if timing: print(f"AVL init took {time.time() - start_time} seconds") + def __set_avl_size_info(self): + # Primary array limits: These also need to updated in the Fortran layer if changed + # its ugly, but it works + + output = self.avl.get_avl_constants() + self.NVMAX = output[0] # number of horseshoe vortices + self.NSMAX = output[1] # number of chord strips + self.NSECMAX = output[2] # nuber of geometry sections + self.NFMAX = output[3] # number of surfaces + self.NLMAX = output[4] # number of source/doublet line nodes + self.NBMAX = output[5] # number of bodies + self.NUMAX = output[6] # number of freestream parameters (V,Omega) + self.NDMAX = output[7] # number of control deflection parameters + self.NGMAX = output[8] # number of design variables + self.NRMAX = output[9] # number of stored run cases + self.NTMAX = output[10] # number of stored time levels + self.NOBMAX = output[11] # max number of off body points + self.ICONX = output[12] # + self.IBX = output[13] # max number of airfoil coordinates in a file + self.NASMAX = output[14] # max airfoil point for representation + def _init_map_data(self): """Used in the __init__ method to allocate the slice data for the surfaces""" self.surf_geom_to_fort_var = {} @@ -612,10 +619,6 @@ def _setup_section_maps(self, surf_name: str, idx_surf: int, num_sec: int, nasec "sasec": ["SURF_GEOM_R", "SASEC", sasec_slices], "casec": ["SURF_GEOM_R", "CASEC", casec_slices], "tasec": ["SURF_GEOM_R", "TASEC", tasec_slices], - "xuasec": ["SURF_GEOM_R", "XUASEC", xuasec_slices], - "xlasec": ["SURF_GEOM_R", "XLASEC", xlasec_slices], - "zuasec": ["SURF_GEOM_R", "ZUASEC", zuasec_slices], - "zlasec": ["SURF_GEOM_R", "ZLASEC", zlasec_slices], } def load_input_dict(self, input_dict: dict, pre_check: bool = True, post_check: bool = False): @@ -913,7 +916,7 @@ def check_type(key, avl_vars, given_val, cast_type=True): ) xfminmax_arr = surf_dict.get("xfminmax", np.tile([0.0, 1.0], (num_secs, 1))) - num_pts = min(50, self.IBX) + num_pts = min(self.NASMAX, self.IBX) # setup for manually specifying coordinates if "xasec" in surf_dict.keys(): @@ -937,9 +940,6 @@ def check_type(key, avl_vars, given_val, cast_type=True): if "xasec" in surf_dict.keys(): self.set_avl_fort_arr("SURF_GEOM_I", "NASEC", nasec_list[j], slicer=(idx_surf, j)) - # TODO-JLA: a user should not have to spcify XLASEC, XUASEC, ZLASEC, ZUASEC - # since these can be calculated from the other inputs - for key in self.surf_section_geom_to_fort_var[surf_name]: avl_vars_secs = self.surf_section_geom_to_fort_var[surf_name][key] avl_vars = (avl_vars_secs[0], avl_vars_secs[1], avl_vars_secs[2][j]) @@ -990,11 +990,6 @@ def check_type(key, avl_vars, given_val, cast_type=True): self.set_avl_fort_arr("SURF_GEOM_R", "XASEC", np.array([0.0, 1.0]), slicer=slicer_airfoil_flat) self.set_avl_fort_arr("SURF_GEOM_R", "SASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) self.set_avl_fort_arr("SURF_GEOM_R", "TASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) - - self.set_avl_fort_arr("SURF_GEOM_R", "XLASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) - self.set_avl_fort_arr("SURF_GEOM_R", "XUASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) - self.set_avl_fort_arr("SURF_GEOM_R", "ZLASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) - self.set_avl_fort_arr("SURF_GEOM_R", "ZUASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) self.set_avl_fort_arr("SURF_GEOM_R", "CASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) # --- setup control variables for each section --- @@ -1295,18 +1290,6 @@ def set_section_naca(self, isec: int, isurf: int, nasec: int, naca: str, xfminma self.set_avl_fort_arr("SURF_GEOM_R", "SASEC", slopes, slicer=(isurf, isec, slice(0, nasec))) self.set_avl_fort_arr("SURF_GEOM_R", "TASEC", thickness, slicer=(isurf, isec, slice(0, nasec))) - self.set_avl_fort_arr( - "SURF_GEOM_R", "XLASEC", xf + 0.5 * thickness * np.sin(theta), slicer=(isurf, isec, slice(0, nasec)) - ) - self.set_avl_fort_arr( - "SURF_GEOM_R", "XUASEC", xf - 0.5 * thickness * np.sin(theta), slicer=(isurf, isec, slice(0, nasec)) - ) - self.set_avl_fort_arr( - "SURF_GEOM_R", "ZLASEC", zf - 0.5 * thickness * np.cos(theta), slicer=(isurf, isec, slice(0, nasec)) - ) - self.set_avl_fort_arr( - "SURF_GEOM_R", "ZUASEC", zf + 0.5 * thickness * np.cos(theta), slicer=(isurf, isec, slice(0, nasec)) - ) self.set_avl_fort_arr("SURF_GEOM_R", "CASEC", zf, slicer=(isurf, isec, slice(0, nasec))) def set_section_coordinates( @@ -3753,7 +3736,7 @@ def clear_ad_seeds_fast(self): num_sec_max = self.get_avl_fort_arr("SURF_GEOM_R", "AINCS").shape[-1] num_sec = np.max(self.get_avl_fort_arr("SURF_GEOM_I", "NSEC")) - num_airfoil_pts_max = self.get_avl_fort_arr("SURF_GEOM_R", "XLASEC").shape[-1] + num_airfoil_pts_max = self.get_avl_fort_arr("SURF_GEOM_R", "CASEC").shape[-1] num_airfoil_pts = np.max(self.get_avl_fort_arr("SURF_GEOM_I", "NASEC")) mesh_size_max = 4*num_vor_max diff --git a/pyproject.toml b/pyproject.toml index 9876685e..18be2b20 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,12 +28,12 @@ dependencies = [ "numpy>=1.19", ] readme = "README.md" -version = "2.4.0" # this automatically updates __init__.py and setup_deprecated.py +version = "2.5.0" # this automatically updates __init__.py and setup_deprecated.py # don't forget to change the value in meson.build! [tool.cibuildwheel] -skip = "cp36-* cp37-* pp* *_ppc64le *_i686 *_s390x" +skip = "cp38-* pp* *_ppc64le *_i686 *_s390x" build-verbosity = "3" test-command = "bash {project}/tools/wheels/cibw_test_command.sh {project}" @@ -62,6 +62,11 @@ select = "*-win_amd64" # environment = { CMAKE_PREFIX_PATH="c:/opt/64" } environment = { PKG_CONFIG_PATH="c:/opt/64/lib/pkgconfig" } +# pyproject.toml +[tool.cibuildwheel.config-settings] +setup-args = ["-Dbuildtype=debug"] # optional, but helps +compile-args = ["-v"] # this is the key one — makes ninja verbose + [project.urls] homepage = "https://github.com/joanibal/OptVL" # documentation = diff --git a/src/ad_src/Makefile_tapenade b/src/ad_src/Makefile_tapenade index 908cfead..0f03c82c 100644 --- a/src/ad_src/Makefile_tapenade +++ b/src/ad_src/Makefile_tapenade @@ -26,12 +26,12 @@ PP_FILES = $(addprefix $(PP_DIR)/,$(notdir $(ALL_RES_FILES))) # you also need to add any new files to `src/build/fileList` # ====================== Full List of Routines ================== fullRoutines = "\ - update_surfaces(XYZSCAL,XYZTRAN,ADDINC,XYZLES,CHORDS,MSHBLK,AINCS,XASEC,SASEC,TASEC,CLCDSEC,CLAF)>(ENC,ENV, DXV, CHORDV,CHORD, CHORD1, CHORD2, RLE,RLE1,RLE2, WSTRIP, RV1,RV2,RV,RC,RS,RL, ENSY, ENSZ, ESS, XSREF,YSREF,ZSREF, ENC_D)\ + update_surfaces(XYZSCAL,XYZTRAN,ADDINC,XYZLES,CHORDS,MSHBLK,AINCS,XASEC,SASEC,TASEC,CLCDSEC,CLAF)>(ENC,ENV, DXV, CHORDV,CHORD, CHORD1, CHORD2, RLE,RLE1,RLE2, WSTRIP, RV1,RV2,RV,RC,RS,RL, ENSY, ENSZ, ESS, XSREF,YSREF,ZSREF, ENC_D, CLCD)\ \ get_res(GAM, GAM_D, GAM_U, CONVAL, PARVAL, YSYM, ZSYM, ENC, ENV, DXV, CHORDV, RV1, RV2, RV, RC, RS, RL, ENC_D, XYZREF)>(RES, RES_D, RES_U, GAM, GAM_D, GAM_U, VINF, WROT, VINF_A, VINF_B, ALFA, BETA, DELCON, RV1,RV2, RV, RC,DXV, XYZREF, WV_GAM, CDREF, MACH, SRC, SRC_U)\ \ VElSUM(GAM, GAM_U, GAM_D, VINF, WROT, WV_GAM)>(VV, VV_U, VV_D, WV, WV_U, WV_D, GAM, GAM_U, GAM_D, VINF, WROT, RV, RV1, RV2) \ - AERO(GAM, GAM_U, GAM_D, VV, VV_U, VV_D, WV, WV_U, WV_D, RC, RV, RV1, RV2, RLE, XYZREF, CHORD, CHORD1, CHORD2, RLE1, RLE2, WSTRIP, ALFA, MACH, VINF, WROT, VINF_A, VINF_B, ENSZ, ENSY, ESS, SREF, CREF, BREF, XSREF,YSREF,ZSREF,DXV, CDREF, SRC, SRC_U)>(\ + AERO(GAM, GAM_U, GAM_D, VV, VV_U, VV_D, WV, WV_U, WV_D, RC, RV, RV1, RV2, RLE, XYZREF, CHORD, CHORD1, CHORD2, RLE1, RLE2, WSTRIP, ALFA, MACH, VINF, WROT, VINF_A, VINF_B, ENSZ, ENSY, ESS, SREF, CREF, BREF, XSREF,YSREF,ZSREF,DXV, CLCD, CDREF, SRC, SRC_U)>(\ CLTOT, CYTOT, CDTOT, CDVTOT, CDITOT\ CLFF, CYFF, CDFF,\ CFTOT, \ diff --git a/src/ad_src/ad_utils/edit_ad_src.py b/src/ad_src/ad_utils/edit_ad_src.py index 284a689f..2e628649 100644 --- a/src/ad_src/ad_utils/edit_ad_src.py +++ b/src/ad_src/ad_utils/edit_ad_src.py @@ -37,10 +37,12 @@ # fmt: off hand_made_include = " INCLUDE 'AVL.INC'\n"\ " INCLUDE 'AVL_ad_seeds.inc'\n" +hand_made_surf_include = " INCLUDE 'AVL_surf.INC'\n"\ + " INCLUDE 'AVL_surf_ad_seeds.inc'\n" hand_made_heap_mod = " use avl_heap_diff_inc\n" heap_mod = " use avl_heap_inc\n" fake_include = " INCLUDE 'AVL_tapenade_fake.inc'\n" -fake_include_diff = " INCLUDE 'AVL_tapenade_fake_diff.inc'\n" +fake_include_diff = " INCLUDE 'AVL_tapenade_fake_diff.inc'\n" # fmt: on @@ -55,10 +57,12 @@ # Specify file extension file_ext = "_b.f" tapenade_include = "INCLUDE 'AVL_b.inc'" + tapenade_surf_include = "INCLUDE 'AVL_surf_b.inc'" fake_diff_include = "INCLUDE 'AVL_tapenade_fake_b.inc'" elif args.mode == "forward": # Specify file extension tapenade_include = "INCLUDE 'AVL_d.inc'" + tapenade_surf_include = "INCLUDE 'AVL_surf_d.inc'" fake_diff_include = "INCLUDE 'AVL_tapenade_fake_d.inc'" file_ext = "_d.f" @@ -85,8 +89,8 @@ def read_nvmax(path): # only loop over nvor vortices 'ii1=1,nvmax': 'ii1=1,nvor', 'ii2=1,nvmax': 'ii2=1,nvor', - f'ii1=1,{nvmax}': 'ii1=1,nvor', - f'ii2=1,{nvmax}': 'ii2=1,nvor', + f'ii1=1,{nvmax}': 'ii1=1,nvor', + f'ii2=1,{nvmax}': 'ii2=1,nvor', 'ii1=1,ncdim': 'ii1=1,nc', 'ii2=1,ncdim': 'ii2=1,nc', @@ -102,9 +106,6 @@ def read_nvmax(path): 'ii1=1,nsmax': 'ii1=1,NSTRIP', 'ii2=1,nsmax': 'ii2=1,NSTRIP', - - - } ad_commented_lines = { 'aoper_b.f': { @@ -157,13 +158,17 @@ def read_nvmax(path): # check to see if the next line is the fake include if fake_diff_include.strip() in lines[i+1]: fid_mod.write(heap_mod) - fid_mod.write(hand_made_heap_mod) + fid_mod.write(hand_made_heap_mod) fid_mod.write(hand_made_include) i += 2 continue else: # Insert the hand-written include file fid_mod.write(hand_made_include) + + elif tapenade_surf_include in line: + # check to see if the next line is the fake include + fid_mod.write(hand_made_surf_include) elif "REAL*(" in line: # This syntax is not supported by our compiler diff --git a/src/ad_src/forward_ad_src/aero_d.f b/src/ad_src/forward_ad_src/aero_d.f index 718de70a..77b1de44 100644 --- a/src/ad_src/forward_ad_src/aero_d.f +++ b/src/ad_src/forward_ad_src/aero_d.f @@ -17,9 +17,9 @@ C rr C with respect to varying inputs: alfa vinf vinf_a vinf_b wrot C sref cref bref xyzref mach cdref rle chord rle1 -C chord1 rle2 chord2 wstrip ensy ensz xsref ysref -C zsref rv1 rv2 rv rc gam gam_u gam_d src src_u -C vv vv_u vv_d wv wv_u wv_d +C chord1 rle2 chord2 wstrip clcd ensy ensz xsref +C ysref zsref rv1 rv2 rv rc gam gam_u gam_d src +C src_u vv vv_u vv_d wv wv_u wv_d C RW status of diff variables: alfa:in vinf:in vinf_a:in vinf_b:in C wrot:in sref:in cref:in bref:in xyzref:in mach:in C cdref:in clff:out cyff:out cdff:out spanef:out @@ -40,9 +40,10 @@ C cltot_rz:out cytot_rz:out crtot_rz:out cmtot_rz:out C cntot_rz:out xnp:out sm:out bb:out rr:out rle:in C chord:in rle1:in chord1:in rle2:in chord2:in wstrip:in -C ensy:in ensz:in xsref:in ysref:in zsref:in rv1:in -C rv2:in rv:in rc:in gam:in gam_u:in gam_d:in src:in -C src_u:in vv:in vv_u:in vv_d:in wv:in wv_u:in wv_d:in +C clcd:in ensy:in ensz:in xsref:in ysref:in zsref:in +C rv1:in rv2:in rv:in rc:in gam:in gam_u:in gam_d:in +C src:in src_u:in vv:in vv_u:in vv_d:in wv:in wv_u:in +C wv_d:in C*********************************************************************** C Module: aero.f C @@ -335,8 +336,8 @@ SUBROUTINE AERO_D() C cftot cftot_u cftot_d cmtot cmtot_u cmtot_d cdvtot C with respect to varying inputs: alfa vinf wrot sref cref bref C xyzref rle chord rle1 chord1 rle2 chord2 wstrip -C ensy ensz xsref ysref zsref rv1 rv2 rv gam gam_u -C gam_d vv vv_u vv_d wv wv_u wv_d +C clcd ensy ensz xsref ysref zsref rv1 rv2 rv gam +C gam_u gam_d vv vv_u vv_d wv wv_u wv_d C AERO C C @@ -1710,8 +1711,8 @@ SUBROUTINE SFFORC_D() ENDDO C C--- Get CD from CLCD function using strip CL as parameter - CALL CDCL_D(clcd(1, j), clv, clv_diff, cdv, cdv_diff, cdv_clv - + , cdv_clv_diff) + CALL CDCL_D(clcd(1, j), clcd_diff(1, j), clv, clv_diff, cdv, + + cdv_diff, cdv_clv, cdv_clv_diff) C C--- Strip viscous force contribution (per unit strip area) dcvfx_diff = veffmag*cdv*veff_diff(1) + veff(1)*(cdv* diff --git a/src/ad_src/forward_ad_src/amake_d.f b/src/ad_src/forward_ad_src/amake_d.f index 17bf2b22..2a894ff7 100644 --- a/src/ad_src/forward_ad_src/amake_d.f +++ b/src/ad_src/forward_ad_src/amake_d.f @@ -3,31 +3,33 @@ C C Differentiation of update_surfaces in forward (tangent) mode (with options i4 dr8 r8): C variations of useful results: rle chord rle1 chord1 rle2 -C chord2 wstrip ess ensy ensz xsref ysref zsref +C chord2 wstrip clcd ess ensy ensz xsref ysref zsref C rv1 rv2 rv rc rs dxv chordv enc env enc_d C with respect to varying inputs: xyzscal xyztran addinc xyzles -C chords aincs xasec sasec claf mshblk -C RW status of diff variables: xyzscal:in xyztran:in addinc:in -C xyzles:in chords:in aincs:in xasec:in sasec:in -C claf:in mshblk:in rle:out chord:out rle1:out chord1:out -C rle2:out chord2:out wstrip:out ess:out ensy:out -C ensz:out xsref:out ysref:out zsref:out rv1:out -C rv2:out rv:out rc:out rs:out dxv:out chordv:out -C enc:out env:out enc_d:out +C chords aincs xasec sasec clcdsec claf mshblk +C RW status of diff variables: rle:out chord:out rle1:out chord1:out +C rle2:out chord2:out wstrip:out clcd:out ess:out +C ensy:out ensz:out xsref:out ysref:out zsref:out +C rv1:out rv2:out rv:out rc:out rs:out dxv:out chordv:out +C enc:out env:out enc_d:out xyzscal:in xyztran:in +C addinc:in xyzles:in chords:in aincs:in xasec:in +C sasec:in clcdsec:in claf:in mshblk:in SUBROUTINE UPDATE_SURFACES_D() use avl_heap_inc use avl_heap_diff_inc INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' INTEGER ii INTEGER isurf EXTERNAL AVLHEAP_CLEAN EXTERNAL AVLHEAP_DIFF_CLEAN EXTERNAL AVLHEAP_INIT EXTERNAL AVLHEAP_DIFF_INIT - INTEGER ii3 INTEGER ii2 INTEGER ii1 + INTEGER ii3 nstrip = 0 nvor = 0 isurf = 1 @@ -35,26 +37,6 @@ SUBROUTINE UPDATE_SURFACES_D() DO ii=1,nsurf IF (ldupl(ii)) nsurfdupl = nsurfdupl + 1 ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv1msh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv2msh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rvmsh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rcmsh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO DO ii1=1,NSTRIP DO ii2=1,3 rle_diff(ii2, ii1) = 0.D0 @@ -82,6 +64,11 @@ SUBROUTINE UPDATE_SURFACES_D() DO ii1=1,NSTRIP wstrip_diff(ii1) = 0.D0 ENDDO + DO ii1=1,NSTRIP + DO ii2=1,6 + clcd_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO DO ii1=1,NSTRIP ainc_diff(ii1) = 0.D0 ENDDO @@ -139,6 +126,26 @@ SUBROUTINE UPDATE_SURFACES_D() ENDDO ENDDO ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rv1msh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rv2msh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rvmsh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rcmsh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO C the iterations of this loop are not independent because we count C up the size information as we make each surface DO ii=1,nsurf-nsurfdupl @@ -173,12 +180,13 @@ SUBROUTINE UPDATE_SURFACES_D() C Differentiation of makesurf in forward (tangent) mode (with options i4 dr8 r8): C variations of useful results: rle chord rle1 chord1 rle2 -C chord2 wstrip ainc ainc_g rv1 rv2 rv rc rs dxv -C chordv slopev slopec dcontrol vhinge -C with respect to varying inputs: xyzscal xyztran addinc xyzles -C chords aincs xasec sasec claf rle chord rle1 chord1 -C rle2 chord2 wstrip ainc ainc_g rv1 rv2 rv rc rs +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs C dxv chordv slopev slopec dcontrol vhinge +C with respect to varying inputs: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs +C dxv chordv slopev slopec dcontrol vhinge xyzscal +C xyztran addinc xyzles chords aincs xasec sasec +C clcdsec claf C*********************************************************************** C Module: amake.f C @@ -202,6 +210,9 @@ SUBROUTINE UPDATE_SURFACES_D() SUBROUTINE MAKESURF_D(isurf) INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' + INCLUDE 'AVL_oml.INC' C C REAL xyzlel(3), xyzler(3) @@ -231,6 +242,10 @@ SUBROUTINE MAKESURF_D(isurf) REAL xled(ndmax), xted(ndmax), gainda(ndmax) REAL xled_diff(ndmax), xted_diff(ndmax), gainda_diff(ndmax) INTEGER idx_vor, idx_strip + REAL xlasecl(nasmax), xuasecl(nasmax) + REAL zlasecl(nasmax), zuasecl(nasmax) + REAL xlasecr(nasmax), xuasecr(nasmax) + REAL zlasecr(nasmax), zuasecr(nasmax) INTEGER isec REAL dy REAL dy_diff @@ -362,7 +377,6 @@ SUBROUTINE MAKESURF_D(isurf) REAL result1 REAL result1_diff REAL temp - INTEGER ii1 REAL temp0 INTEGER isurf C @@ -409,9 +423,7 @@ SUBROUTINE MAKESURF_D(isurf) C----------------------------------------------------------------- C---- Arc length positions of sections in wing trace in y-z plane yzlen(1) = 0. - DO ii1=1,ksmax - yzlen_diff(ii1) = 0.D0 - ENDDO + yzlen_diff = 0.D0 DO isec=2,nsec(isurf) dy_diff = xyzles_diff(2, isec, isurf) - xyzles_diff(2, isec-1 + , isurf) @@ -447,15 +459,11 @@ SUBROUTINE MAKESURF_D(isurf) ELSE C nvs(isurf) = 0 - DO ii1=1,ksmax - ypt_diff(ii1) = 0.D0 - ENDDO + ypt_diff = 0.D0 ypt_diff(1) = yzlen_diff(1) ypt(1) = yzlen(1) iptloc(1) = 1 - DO ii1=1,ksmax - ycp_diff(ii1) = 0.D0 - ENDDO + ycp_diff = 0.D0 C DO isec=1,nsec(isurf)-1 dyzlen_diff = yzlen_diff(isec+1) - yzlen_diff(isec) @@ -509,14 +517,10 @@ SUBROUTINE MAKESURF_D(isurf) ELSE CALL SPACER(nspace, sspace(isurf), fspace) C - DO ii1=1,ksmax - ypt_diff(ii1) = 0.D0 - ENDDO + ypt_diff = 0.D0 ypt_diff(1) = yzlen_diff(1) ypt(1) = yzlen(1) - DO ii1=1,ksmax - ycp_diff(ii1) = 0.D0 - ENDDO + ycp_diff = 0.D0 DO ivs=1,nvs(isurf) ycp_diff(ivs) = yzlen_diff(1) + fspace(2*ivs)*(yzlen_diff( + nsec(isurf))-yzlen_diff(1)) @@ -647,36 +651,16 @@ SUBROUTINE MAKESURF_D(isurf) + ndesign STOP ELSE - DO ii1=1,3 - xyzler_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,kcmax - xcp_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chcosl_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chsinr_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - xted_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - xled_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chsinl_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - gainda_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chcosr_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,3 - xyzlel_diff(ii1) = 0.D0 - ENDDO + xyzler_diff = 0.D0 + xcp_diff = 0.D0 + chcosl_g_diff = 0.D0 + chsinr_g_diff = 0.D0 + xted_diff = 0.D0 + xled_diff = 0.D0 + chsinl_g_diff = 0.D0 + gainda_diff = 0.D0 + chcosr_g_diff = 0.D0 + xyzlel_diff = 0.D0 C C---- go over section intervals DO isec=1,nsec(isurf)-1 @@ -1070,6 +1054,9 @@ SUBROUTINE MAKESURF_D(isurf) C C--- Interpolate CD-CL polar defining data from input sections to strips DO l=1,6 + clcd_diff(l, idx_strip) = (1.0-fc)*clcdsec_diff(l, isec + + , isurf) - (clcdsec(l, isec, isurf)-clcdsec(l, isec+1 + + , isurf))*fc_diff + fc*clcdsec_diff(l, isec+1, isurf) clcd(l, idx_strip) = (1.0-fc)*clcdsec(l, isec, isurf) + + fc*clcdsec(l, isec+1, isurf) ENDDO @@ -1173,16 +1160,12 @@ SUBROUTINE MAKESURF_D(isurf) + fc*chordr*(sloper_diff-temp*chordc_diff)/chordc slopec(idx_vor) = temp0*(chordl*slopel) + fc*chordr*temp C - DO ii1=1,kcmax - xvr_diff(ii1) = 0.D0 - ENDDO + xvr_diff = 0.D0 CALL AKIMA_D(xasec(1, isec, isurf), xasec_diff(1, isec, + isurf), sasec(1, isec, isurf), sasec_diff(1 + , isec, isurf), nsl, xvr(ivc), xvr_diff(ivc + ), slopel, slopel_diff, dsdx) - DO ii1=1,kcmax - xvr_diff(ii1) = 0.D0 - ENDDO + xvr_diff = 0.D0 CALL AKIMA_D(xasec(1, isec+1, isurf), xasec_diff(1, isec + +1, isurf), sasec(1, isec+1, isurf), + sasec_diff(1, isec+1, isurf), nsr, xvr(ivc) @@ -1258,15 +1241,29 @@ SUBROUTINE MAKESURF_D(isurf) C... nodal grid associated with vortex strip (aft-panel nodes) C... NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral; C... retained in (x,z) plane at this point - CALL AKIMA(xlasec(1, isec, isurf), zlasec(1, isec, isurf - + ), nsl, xpt(ivc+1), zl_l, dsdx) - CALL AKIMA(xuasec(1, isec, isurf), zuasec(1, isec, isurf - + ), nsl, xpt(ivc+1), zu_l, dsdx) -C - CALL AKIMA(xlasec(1, isec+1, isurf), zlasec(1, isec+1, - + isurf), nsr, xpt(ivc+1), zl_r, dsdx) - CALL AKIMA(xuasec(1, isec+1, isurf), zuasec(1, isec+1, - + isurf), nsr, xpt(ivc+1), zu_r, dsdx) + xlasecl = xasec(:, isec, isurf) + xuasecl = xasec(:, isec, isurf) + zlasecl = casec(:, isec, isurf) - 0.5*tasec(:, isec, + + isurf) + zuasecl = casec(:, isec, isurf) + 0.5*tasec(:, isec, + + isurf) + xlasecr = xasec(:, isec, isurf) + xuasecr = xasec(:, isec, isurf) + zlasecr = casec(:, isec, isurf) - 0.5*tasec(:, isec, + + isurf) + zuasecr = casec(:, isec, isurf) + 0.5*tasec(:, isec, + + isurf) +C +C + CALL AKIMA(xlasecl, zlasecl, nsl, xpt(ivc+1), zl_l, dsdx + + ) + CALL AKIMA(xuasecl, zuasecl, nsl, xpt(ivc+1), zu_l, dsdx + + ) +C + CALL AKIMA(xlasecr, zlasecr, nsr, xpt(ivc+1), zl_r, dsdx + + ) + CALL AKIMA(xuasecr, zuasecr, nsr, xpt(ivc+1), zu_r, dsdx + + ) C xyn1(1, idx_vor) = rle1(1, idx_strip) + xpt(ivc+1)* + chord1(idx_strip) @@ -1329,19 +1326,22 @@ SUBROUTINE MAKESURF_D(isurf) END C Differentiation of makesurf_mesh in forward (tangent) mode (with options i4 dr8 r8): -C variations of useful results: rv1msh rv2msh rvmsh rcmsh rle -C chord rle1 chord1 rle2 chord2 wstrip ainc ainc_g -C rv1 rv2 rv rc rs dxv chordv slopev slopec dcontrol -C vhinge -C with respect to varying inputs: xyzscal xyztran addinc aincs -C xasec sasec claf mshblk rv1msh rv2msh rvmsh rcmsh -C rle chord rle1 chord1 rle2 chord2 wstrip ainc -C ainc_g rv1 rv2 rv rc rs dxv chordv slopev slopec -C dcontrol vhinge +C variations of useful results: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs +C dxv chordv slopev slopec dcontrol vhinge rv1msh +C rv2msh rvmsh rcmsh +C with respect to varying inputs: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs +C dxv chordv slopev slopec dcontrol vhinge xyzscal +C xyztran addinc aincs xasec sasec clcdsec claf +C mshblk rv1msh rv2msh rvmsh rcmsh C SUBROUTINE MAKESURF_MESH_D(isurf) INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' + INCLUDE 'AVL_oml.INC' C working variables (AVL original) INTEGER isurf INTEGER kcmax @@ -1358,8 +1358,12 @@ SUBROUTINE MAKESURF_MESH_D(isurf) REAL chsinl_g_diff(ngmax), chcosl_g_diff(ngmax), chsinr_g_diff( + ngmax), chcosr_g_diff(ngmax), xled_diff(ndmax), xted_diff( + ndmax) -C working variables (OptVL additions) INTEGER isconl(ndmax), isconr(ndmax) + REAL xlasecl(nasmax), xuasecl(nasmax) + REAL zlasecl(nasmax), zuasecl(nasmax) + REAL xlasecr(nasmax), xuasecr(nasmax) +C working variables (OptVL additions) + REAL zlasecr(nasmax), zuasecr(nasmax) REAL m1, m2, m3, f1, f2, fc, dc1, dc2, dc, a1, a2, a3, xptxind1, + xptxind2 REAL m2_diff, m3_diff, dc1_diff, dc2_diff, a1_diff, a2_diff, @@ -1982,6 +1986,9 @@ SUBROUTINE MAKESURF_MESH_D(isurf) C Interpolate CD-CL polar defining data from input to strips C DO idx_coef=1,6 + clcd_diff(idx_coef, idx_strip) = (1.0-fc)*clcdsec_diff( + + idx_coef, iptl, isurf) + fc*clcdsec_diff(idx_coef, iptr, + + isurf) clcd(idx_coef, idx_strip) = (1.0-fc)*clcdsec(idx_coef, iptl + , isurf) + fc*clcdsec(idx_coef, iptr, isurf) ENDDO @@ -2535,21 +2542,28 @@ SUBROUTINE MAKESURF_MESH_D(isurf) C C xptxind2 = (mesh_surf(1,idx_node_yp1+1) C & - RLE2(1,idx_strip))/CHORD2(idx_strip) -C Interpolate cross section on left side xptxind1 = ((mesh_surf(1, idx_node+1)+mesh_surf(1, + idx_node_yp1+1))/2-rle(1, idx_strip))/chord(idx_strip) C C - CALL AKIMA(xlasec(1, iptl, isurf), zlasec(1, iptl, isurf), - + nsl, xptxind1, zl_l, dsdx) - CALL AKIMA(xuasec(1, iptl, isurf), zuasec(1, iptl, isurf), - + nsl, xptxind1, zu_l, dsdx) +C ! Interpolate cross section on left side +C + xlasecl = xasec(:, iptl, isurf) + xuasecl = xasec(:, iptl, isurf) + zlasecl = casec(:, iptl, isurf) - 0.5*tasec(:, iptl, isurf) + zuasecl = casec(:, iptl, isurf) + 0.5*tasec(:, iptl, isurf) + xlasecr = xasec(:, iptr, isurf) + xuasecr = xasec(:, iptr, isurf) + zlasecr = casec(:, iptr, isurf) - 0.5*tasec(:, iptr, isurf) + zuasecr = casec(:, iptr, isurf) + 0.5*tasec(:, iptr, isurf) +C +C + CALL AKIMA(xlasecl, zlasecl, nsl, xptxind1, zl_l, dsdx) + CALL AKIMA(xuasecl, zuasecl, nsl, xptxind1, zu_l, dsdx) C Interpolate cross section on right side C - CALL AKIMA(xlasec(1, iptr, isurf), zlasec(1, iptr, isurf), - + nsr, xptxind1, zl_r, dsdx) - CALL AKIMA(xuasec(1, iptr, isurf), zuasec(1, iptr, isurf), - + nsr, xptxind1, zu_r, dsdx) + CALL AKIMA(xlasecr, zlasecr, nsr, xptxind1, zl_r, dsdx) + CALL AKIMA(xuasecr, zuasecr, nsr, xptxind1, zu_r, dsdx) C Compute the left aft node of panel C X-point C @@ -2614,18 +2628,21 @@ SUBROUTINE MAKESURF_MESH_D(isurf) END C Differentiation of sdupl in forward (tangent) mode (with options i4 dr8 r8): -C variations of useful results: rv1msh rv2msh rvmsh rcmsh rle -C chord rle1 chord1 rle2 chord2 wstrip ainc ainc_g -C rv1 rv2 rv rc dxv chordv slopev slopec dcontrol -C vhinge -C with respect to varying inputs: rv1msh rv2msh rvmsh rcmsh rle -C chord rle1 chord1 rle2 chord2 wstrip ainc ainc_g -C rv1 rv2 rv rc dxv chordv slopev slopec dcontrol -C vhinge +C variations of useful results: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc dxv +C chordv slopev slopec dcontrol vhinge rv1msh rv2msh +C rvmsh rcmsh +C with respect to varying inputs: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc dxv +C chordv slopev slopec dcontrol vhinge rv1msh rv2msh +C rvmsh rcmsh C SUBROUTINE SDUPL_D(nn, ypt, msg) INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' + INCLUDE 'AVL_oml.INC' CHARACTER*(*) msg INTEGER idx_vor INTEGER nni @@ -2789,6 +2806,7 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) C nvstrp(jji) = nvc(nni) DO l=1,6 + clcd_diff(l, jji) = clcd_diff(l, jj) clcd(l, jji) = clcd(l, jj) ENDDO lviscstrp(jji) = lviscstrp(jj) @@ -2912,8 +2930,8 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) C Differentiation of encalc in forward (tangent) mode (with options i4 dr8 r8): C variations of useful results: ess ensy ensz xsref ysref zsref C enc env enc_d -C with respect to varying inputs: rv1msh rv2msh rvmsh rcmsh ainc -C ainc_g rv1 rv2 rv slopev slopec dcontrol vhinge +C with respect to varying inputs: ainc ainc_g rv1 rv2 rv slopev +C slopec dcontrol vhinge rv1msh rv2msh rvmsh rcmsh C BDUPL C C @@ -2924,6 +2942,8 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) SUBROUTINE ENCALC_D() INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' C REAL ep(3), eq(3), es(3), eb(3), ec(3), ecxb(3) REAL ep_diff(3), eq_diff(3), es_diff(3), eb_diff(3), ec_diff(3), diff --git a/src/ad_src/forward_ad_src/asetup_d.f b/src/ad_src/forward_ad_src/asetup_d.f index 94ef1ec9..4d6a6b7e 100644 --- a/src/ad_src/forward_ad_src/asetup_d.f +++ b/src/ad_src/forward_ad_src/asetup_d.f @@ -169,49 +169,46 @@ SUBROUTINE VELSUM_D() C C DO i=1,nvor +C VC(K,I) = 0.0 DO k=1,3 - vc(k, i) = 0.0 vv_diff(k, i) = 0.D0 vv(k, i) = 0.0 C--- h.v. velocity at control points and vortex midpoints +C VC(K,I) = VC(K,I) + WC_GAM(K,I,J)*GAM(J) DO j=1,nvor - vc(k, i) = vc(k, i) + wc_gam(k, i, j)*gam(j) vv_diff(k, i) = vv_diff(k, i) + gam(j)*wv_gam_diff(k, i, j) + + wv_gam(k, i, j)*gam_diff(j) vv(k, i) = vv(k, i) + wv_gam(k, i, j)*gam(j) ENDDO +C VC_U(K,I,N) = 0.0 DO n=1,numax - vc_u(k, i, n) = 0.0 vv_u_diff(k, i, n) = 0.D0 vv_u(k, i, n) = 0.0 +C VC_U(K,I,N) = VC_U(K,I,N) + WC_GAM(K,I,J)*GAM_U(J,N) DO j=1,nvor - vc_u(k, i, n) = vc_u(k, i, n) + wc_gam(k, i, j)*gam_u(j, n - + ) vv_u_diff(k, i, n) = vv_u_diff(k, i, n) + gam_u(j, n)* + wv_gam_diff(k, i, j) + wv_gam(k, i, j)*gam_u_diff(j, n) vv_u(k, i, n) = vv_u(k, i, n) + wv_gam(k, i, j)*gam_u(j, n + ) ENDDO ENDDO +C VC_D(K,I,N) = 0.0 DO n=1,ncontrol - vc_d(k, i, n) = 0.0 vv_d_diff(k, i, n) = 0.D0 vv_d(k, i, n) = 0.0 +C VC_D(K,I,N) = VC_D(K,I,N) + WC_GAM(K,I,J)*GAM_D(J,N) DO j=1,nvor - vc_d(k, i, n) = vc_d(k, i, n) + wc_gam(k, i, j)*gam_d(j, n - + ) vv_d_diff(k, i, n) = vv_d_diff(k, i, n) + gam_d(j, n)* + wv_gam_diff(k, i, j) + wv_gam(k, i, j)*gam_d_diff(j, n) vv_d(k, i, n) = vv_d(k, i, n) + wv_gam(k, i, j)*gam_d(j, n + ) ENDDO ENDDO +C VC_G(K,I,N) = 0.0 DO n=1,ndesign - vc_g(k, i, n) = 0.0 vv_g(k, i, n) = 0.0 +C VC_G(K,I,N) = VC_G(K,I,N) + WC_GAM(K,I,J)*GAM_G(J,N) DO j=1,nvor - vc_g(k, i, n) = vc_g(k, i, n) + wc_gam(k, i, j)*gam_g(j, n - + ) vv_g(k, i, n) = vv_g(k, i, n) + wv_gam(k, i, j)*gam_g(j, n + ) ENDDO @@ -220,6 +217,7 @@ SUBROUTINE VELSUM_D() wcsrd(k, i) = wcsrd_u(k, i, 1)*vinf(1) + wcsrd_u(k, i, 2)*vinf + (2) + wcsrd_u(k, i, 3)*vinf(3) + wcsrd_u(k, i, 4)*wrot(1) + + wcsrd_u(k, i, 5)*wrot(2) + wcsrd_u(k, i, 6)*wrot(3) +C WC(K,I) = VC(K,I) + WCSRD(K,I) wvsrd_diff(k, i) = wvsrd_u(k, i, 1)*vinf_diff(1) + wvsrd_u(k, + i, 2)*vinf_diff(2) + wvsrd_u(k, i, 3)*vinf_diff(3) + wvsrd_u + (k, i, 4)*wrot_diff(1) + wvsrd_u(k, i, 5)*wrot_diff(2) + @@ -228,21 +226,20 @@ SUBROUTINE VELSUM_D() + (2) + wvsrd_u(k, i, 3)*vinf(3) + wvsrd_u(k, i, 4)*wrot(1) + + wvsrd_u(k, i, 5)*wrot(2) + wvsrd_u(k, i, 6)*wrot(3) C--- total velocity at control points and vortex midpoints - wc(k, i) = vc(k, i) + wcsrd(k, i) wv_diff(k, i) = vv_diff(k, i) + wvsrd_diff(k, i) wv(k, i) = vv(k, i) + wvsrd(k, i) +C WC_U(K,I,N) = VC_U(K,I,N) + WCSRD_U(K,I,N) DO n=1,numax - wc_u(k, i, n) = vc_u(k, i, n) + wcsrd_u(k, i, n) wv_u_diff(k, i, n) = vv_u_diff(k, i, n) wv_u(k, i, n) = vv_u(k, i, n) + wvsrd_u(k, i, n) ENDDO +C WC_D(K,I,N) = VC_D(K,I,N) DO n=1,ndmax - wc_d(k, i, n) = vc_d(k, i, n) wv_d_diff(k, i, n) = vv_d_diff(k, i, n) wv_d(k, i, n) = vv_d(k, i, n) ENDDO +C WC_G(K,I,N) = VC_G(K,I,N) DO n=1,ngmax - wc_g(k, i, n) = vc_g(k, i, n) wv_g(k, i, n) = vv_g(k, i, n) ENDDO ENDDO diff --git a/src/ad_src/forward_ad_src/cdcl_d.f b/src/ad_src/forward_ad_src/cdcl_d.f index cf7d1b63..1e3aa10e 100644 --- a/src/ad_src/forward_ad_src/cdcl_d.f +++ b/src/ad_src/forward_ad_src/cdcl_d.f @@ -3,7 +3,7 @@ C C Differentiation of cdcl in forward (tangent) mode (with options i4 dr8 r8): C variations of useful results: cd_cl cd -C with respect to varying inputs: cl +C with respect to varying inputs: cl cdclpol C*********************************************************************** C Module: cdcl.f C @@ -24,8 +24,8 @@ C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. C*********************************************************************** C - SUBROUTINE CDCL_D(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, - + cd_cl_diff) + SUBROUTINE CDCL_D(cdclpol, cdclpol_diff, cl, cl_diff, cd, cd_diff + + , cd_cl, cd_cl_diff) C------------------------------------------------------------------------- C Returns cd as a function of cl, or drag polar. C @@ -59,16 +59,26 @@ SUBROUTINE CDCL_D(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, C cdclpol(6) cd (cdpos) at clpos C------------------------------------------------------------------------- REAL cdclpol(6) + REAL cdclpol_diff(6) REAL clmin + REAL clmin_diff REAL cdmin + REAL cdmin_diff REAL cl0 + REAL cl0_diff REAL cd0 + REAL cd0_diff REAL clmax + REAL clmax_diff REAL cdmax + REAL cdmax_diff REAL cdx1 + REAL cdx1_diff REAL cdx2 + REAL cdx2_diff REAL clfac REAL temp + REAL temp0 REAL clinc REAL cd_cl REAL cd_cl_diff @@ -87,11 +97,17 @@ SUBROUTINE CDCL_D(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, cd_cl = 0.0 C C---- unpack the cd,cl parameters for the polar + clmin_diff = cdclpol_diff(1) clmin = cdclpol(1) + cdmin_diff = cdclpol_diff(2) cdmin = cdclpol(2) + cl0_diff = cdclpol_diff(3) cl0 = cdclpol(3) + cd0_diff = cdclpol_diff(4) cd0 = cdclpol(4) + clmax_diff = cdclpol_diff(5) clmax = cdclpol(5) + cdmax_diff = cdclpol_diff(6) cdmax = cdclpol(6) IF (clmax .LE. cl0 .OR. cl0 .LE. clmin) THEN WRITE(*, *) '* CDCL: input CL data out of order' @@ -101,45 +117,73 @@ SUBROUTINE CDCL_D(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, ELSE C C---- matching parameters that make the slopes smooth at the stall joins - cdx1 = 2.0*(cdmin-cd0)*(clmin-cl0)/(clmin-cl0)**2 - cdx2 = 2.0*(cdmax-cd0)*(clmax-cl0)/(clmax-cl0)**2 + temp = (cdmin-cd0)/(clmin-cl0) + cdx1_diff = 2.0*(cdmin_diff-cd0_diff-temp*(clmin_diff-cl0_diff)) + + /(clmin-cl0) + cdx1 = 2.0*temp + temp = (cdmax-cd0)/(clmax-cl0) + cdx2_diff = 2.0*(cdmax_diff-cd0_diff-temp*(clmax_diff-cl0_diff)) + + /(clmax-cl0) + cdx2 = 2.0*temp clfac = 1.0/clinc C IF (cl .LT. clmin) THEN C------ negative stall region C- slope matches lower side, quadratic drag rise temp = cdinc*(clfac*clfac) - cd_diff = (temp*2*(cl-clmin)-cdx1/(clmin-cl0))*cl_diff - cd = cdmin + temp*((cl-clmin)*(cl-clmin)) + cdx1*(1.0-(cl-cl0) - + /(clmin-cl0)) - cd_cl_diff = clfac**2*cdinc*2.0*cl_diff - cd_cl = cdinc*clfac**2*(cl-clmin)*2.0 - cdx1/(clmin-cl0) + temp0 = (cl-cl0)/(clmin-cl0) + cd_diff = cdmin_diff + temp*2*(cl-clmin)*(cl_diff-clmin_diff) + + + (1.0-temp0)*cdx1_diff - cdx1*(cl_diff-cl0_diff-temp0*( + + clmin_diff-cl0_diff))/(clmin-cl0) + cd = cdmin + temp*((cl-clmin)*(cl-clmin)) + cdx1*(1.0-temp0) + temp0 = cdx1/(clmin-cl0) + cd_cl_diff = clfac**2*cdinc*2.0*(cl_diff-clmin_diff) - ( + + cdx1_diff-temp0*(clmin_diff-cl0_diff))/(clmin-cl0) + cd_cl = cdinc*2.0*(cl-clmin)*(clfac*clfac) - temp0 C ELSE IF (cl .LT. cl0) THEN C------ lower quadratic C- slope matches negative stall, and has zero slope at min drag point - cd_diff = (cdmin-cd0)*2*(cl-cl0)*cl_diff/(clmin-cl0)**2 - cd = cd0 + (cdmin-cd0)*(cl-cl0)**2/(clmin-cl0)**2 - cd_cl_diff = (cdmin-cd0)*2.0*cl_diff/(clmin-cl0)**2 - cd_cl = (cdmin-cd0)*(cl-cl0)*2.0/(clmin-cl0)**2 + temp0 = (cdmin-cd0)*((cl-cl0)*(cl-cl0))/((clmin-cl0)*(clmin- + + cl0)) + cd_diff = cd0_diff + ((cl-cl0)**2*(cdmin_diff-cd0_diff)+(cdmin + + -cd0)*2*(cl-cl0)*(cl_diff-cl0_diff)-temp0*2*(clmin-cl0)*( + + clmin_diff-cl0_diff))/(clmin-cl0)**2 + cd = cd0 + temp0 + temp0 = (cdmin-cd0)*(cl-cl0)/((clmin-cl0)*(clmin-cl0)) + cd_cl_diff = 2.0*((cl-cl0)*(cdmin_diff-cd0_diff)+(cdmin-cd0)*( + + cl_diff-cl0_diff)-temp0*2*(clmin-cl0)*(clmin_diff-cl0_diff)) + + /(clmin-cl0)**2 + cd_cl = 2.0*temp0 C ELSE IF (cl .LT. clmax) THEN C------ upper quadratic C- slope matches positive stall, and has zero slope at min drag point - cd_diff = (cdmax-cd0)*2*(cl-cl0)*cl_diff/(clmax-cl0)**2 - cd = cd0 + (cdmax-cd0)*(cl-cl0)**2/(clmax-cl0)**2 - cd_cl_diff = (cdmax-cd0)*2.0*cl_diff/(clmax-cl0)**2 - cd_cl = (cdmax-cd0)*(cl-cl0)*2.0/(clmax-cl0)**2 + temp0 = (cdmax-cd0)*((cl-cl0)*(cl-cl0))/((clmax-cl0)*(clmax- + + cl0)) + cd_diff = cd0_diff + ((cl-cl0)**2*(cdmax_diff-cd0_diff)+(cdmax + + -cd0)*2*(cl-cl0)*(cl_diff-cl0_diff)-temp0*2*(clmax-cl0)*( + + clmax_diff-cl0_diff))/(clmax-cl0)**2 + cd = cd0 + temp0 + temp0 = (cdmax-cd0)*(cl-cl0)/((clmax-cl0)*(clmax-cl0)) + cd_cl_diff = 2.0*((cl-cl0)*(cdmax_diff-cd0_diff)+(cdmax-cd0)*( + + cl_diff-cl0_diff)-temp0*2*(clmax-cl0)*(clmax_diff-cl0_diff)) + + /(clmax-cl0)**2 + cd_cl = 2.0*temp0 C ELSE C------ positive stall region C- slope matches upper side, quadratic drag rise - temp = cdinc*(clfac*clfac) - cd_diff = (temp*2*(cl-clmax)+cdx2/(clmax-cl0))*cl_diff - cd = cdmax + temp*((cl-clmax)*(cl-clmax)) - cdx2*(1.0-(cl-cl0) - + /(clmax-cl0)) - cd_cl_diff = clfac**2*cdinc*2.0*cl_diff - cd_cl = cdinc*clfac**2*(cl-clmax)*2.0 + cdx2/(clmax-cl0) + temp0 = cdinc*(clfac*clfac) + temp = (cl-cl0)/(clmax-cl0) + cd_diff = cdmax_diff + temp0*2*(cl-clmax)*(cl_diff-clmax_diff) + + - (1.0-temp)*cdx2_diff + cdx2*(cl_diff-cl0_diff-temp*( + + clmax_diff-cl0_diff))/(clmax-cl0) + cd = cdmax + temp0*((cl-clmax)*(cl-clmax)) - cdx2*(1.0-temp) + temp0 = cdx2/(clmax-cl0) + cd_cl_diff = clfac**2*cdinc*2.0*(cl_diff-clmax_diff) + ( + + cdx2_diff-temp0*(clmax_diff-cl0_diff))/(clmax-cl0) + cd_cl = cdinc*2.0*(cl-clmax)*(clfac*clfac) + temp0 END IF C RETURN diff --git a/src/ad_src/reverse_ad_src/aero_b.f b/src/ad_src/reverse_ad_src/aero_b.f index d9f76a25..af9aa7b9 100644 --- a/src/ad_src/reverse_ad_src/aero_b.f +++ b/src/ad_src/reverse_ad_src/aero_b.f @@ -28,9 +28,9 @@ C cytot_rx crtot_rx cmtot_rx cntot_rx cdtot_ry cltot_ry C cytot_ry crtot_ry cmtot_ry cntot_ry cdtot_rz cltot_rz C cytot_rz crtot_rz cmtot_rz cntot_rz xnp sm bb -C rr rle chord rle1 chord1 rle2 chord2 wstrip ensy -C ensz xsref ysref zsref rv1 rv2 rv rc gam gam_u -C gam_d src src_u vv vv_u vv_d wv wv_u wv_d +C rr rle chord rle1 chord1 rle2 chord2 wstrip clcd +C ensy ensz xsref ysref zsref rv1 rv2 rv rc gam +C gam_u gam_d src src_u vv vv_u vv_d wv wv_u wv_d C RW status of diff variables: alfa:out vinf:out vinf_a:out vinf_b:out C wrot:out sref:out cref:out bref:out xyzref:out C mach:out cdref:out clff:in-zero cyff:in-zero cdff:in-zero @@ -56,11 +56,11 @@ C cytot_rz:in-zero crtot_rz:in-zero cmtot_rz:in-zero C cntot_rz:in-zero xnp:in-out sm:in-out bb:in-out C rr:in-out rle:out chord:out rle1:out chord1:out -C rle2:out chord2:out wstrip:out ensy:out ensz:out -C xsref:out ysref:out zsref:out rv1:out rv2:out -C rv:out rc:out gam:out gam_u:out gam_d:out src:out -C src_u:out vv:out vv_u:out vv_d:out wv:out wv_u:out -C wv_d:out +C rle2:out chord2:out wstrip:out clcd:out ensy:out +C ensz:out xsref:out ysref:out zsref:out rv1:out +C rv2:out rv:out rc:out gam:out gam_u:out gam_d:out +C src:out src_u:out vv:out vv_u:out vv_d:out wv:out +C wv_u:out wv_d:out C*********************************************************************** C Module: aero.f C @@ -435,8 +435,8 @@ SUBROUTINE AERO_B() C with respect to varying inputs: alfa vinf wrot sref cref bref C xyzref cdtot_d cytot_d cltot_d cftot cftot_d cmtot C cmtot_d rle chord rle1 chord1 rle2 chord2 wstrip -C ensy ensz xsref ysref zsref rv1 rv2 rv gam gam_u -C gam_d vv vv_u vv_d wv wv_u wv_d +C clcd ensy ensz xsref ysref zsref rv1 rv2 rv gam +C gam_u gam_d vv vv_u vv_d wv wv_u wv_d C AERO C C @@ -2004,6 +2004,11 @@ SUBROUTINE SFFORC_B() DO ii1=1,NSTRIP chord2_diff(ii1) = 0.D0 ENDDO + DO ii1=1,NSTRIP + DO ii2=1,6 + clcd_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO DO ii1=1,NSTRIP ensy_diff(ii1) = 0.D0 ENDDO @@ -3296,8 +3301,8 @@ SUBROUTINE SFFORC_B() + veff(2)*dcvfy_diff + veffmag*veff(1)*dcvfx_diff veff_diff(2) = veff_diff(2) + veffmag*cdv*dcvfy_diff veff_diff(1) = veff_diff(1) + veffmag*cdv*dcvfx_diff - CALL CDCL_B(clcd(1, j), clv, clv_diff, cdv, cdv_diff, cdv_clv - + , cdv_clv_diff) + CALL CDCL_B(clcd(1, j), clcd_diff(1, j), clv, clv_diff, cdv, + + cdv_diff, cdv_clv, cdv_clv_diff) CALL POPREAL8ARRAY(cfy_d, ndmax) CALL POPREAL8ARRAY(cfz_d, ndmax) CALL POPREAL8ARRAY(cfx_d, ndmax) diff --git a/src/ad_src/reverse_ad_src/amake_b.f b/src/ad_src/reverse_ad_src/amake_b.f index abb90cbe..12ff736c 100644 --- a/src/ad_src/reverse_ad_src/amake_b.f +++ b/src/ad_src/reverse_ad_src/amake_b.f @@ -3,26 +3,29 @@ C C Differentiation of update_surfaces in reverse (adjoint) mode (with options i4 dr8 r8): C gradient of useful results: rle chord rle1 chord1 rle2 -C chord2 wstrip ess ensy ensz xsref ysref zsref +C chord2 wstrip clcd ess ensy ensz xsref ysref zsref C rv1 rv2 rv rc rs dxv chordv enc env enc_d -C with respect to varying inputs: xyzscal xyztran addinc xyzles -C chords aincs xasec sasec claf mshblk rle chord -C rle1 chord1 rle2 chord2 wstrip ess ensy ensz xsref -C ysref zsref rv1 rv2 rv rc rs dxv chordv enc env -C enc_d -C RW status of diff variables: xyzscal:out xyztran:out addinc:out -C xyzles:out chords:out aincs:out xasec:out sasec:out -C claf:out mshblk:out rle:in-out chord:in-out rle1:in-out +C with respect to varying inputs: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ess ensy ensz xsref ysref zsref +C rv1 rv2 rv rc rs dxv chordv enc env enc_d xyzscal +C xyztran addinc xyzles chords aincs xasec sasec +C clcdsec claf mshblk +C RW status of diff variables: rle:in-out chord:in-out rle1:in-out C chord1:in-out rle2:in-out chord2:in-out wstrip:in-out -C ess:in-out ensy:in-out ensz:in-out xsref:in-out -C ysref:in-out zsref:in-out rv1:in-out rv2:in-out -C rv:in-out rc:in-out rs:in-out dxv:in-out chordv:in-out -C enc:in-out env:in-out enc_d:in-out +C clcd:in-out ess:in-out ensy:in-out ensz:in-out +C xsref:in-out ysref:in-out zsref:in-out rv1:in-out +C rv2:in-out rv:in-out rc:in-out rs:in-out dxv:in-out +C chordv:in-out enc:in-out env:in-out enc_d:in-out +C xyzscal:out xyztran:out addinc:out xyzles:out +C chords:out aincs:out xasec:out sasec:out clcdsec:out +C claf:out mshblk:out SUBROUTINE UPDATE_SURFACES_B() use avl_heap_inc use avl_heap_diff_inc INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' INTEGER ii INTEGER isurf EXTERNAL AVLHEAP_CLEAN @@ -50,6 +53,8 @@ SUBROUTINE UPDATE_SURFACES_B() DO ii=1,nsurf-nsurfdupl IF (lverbose) WRITE(*, *) 'Updating surface ', isurf IF (lsurfmsh(isurf)) THEN + CALL PUSHINTEGER4ARRAY(nvs, nfmax) + CALL PUSHINTEGER4ARRAY(nvc, nfmax) CALL PUSHREAL8ARRAY(dxv, nvmax) CALL PUSHREAL8ARRAY(rc, 3*nvmax) CALL PUSHREAL8ARRAY(rv, 3*nvmax) @@ -61,18 +66,16 @@ SUBROUTINE UPDATE_SURFACES_B() CALL PUSHREAL8ARRAY(rle, 3*nsmax) CALL PUSHINTEGER4ARRAY(nvstrp, nsmax) CALL PUSHINTEGER4ARRAY(ijfrst, nsmax) - CALL PUSHINTEGER4ARRAY(nvs, nfmax) - CALL PUSHINTEGER4ARRAY(nvc, nfmax) CALL PUSHINTEGER4ARRAY(jfrst, nfmax) CALL PUSHINTEGER4ARRAY(nj, nfmax) CALL MAKESURF_MESH(isurf) CALL PUSHCONTROL1B(0) ELSE + CALL PUSHINTEGER4ARRAY(nvs, nfmax) + CALL PUSHINTEGER4ARRAY(nvc, nfmax) CALL PUSHREAL8ARRAY(chord, nsmax) CALL PUSHINTEGER4ARRAY(nvstrp, nsmax) CALL PUSHINTEGER4ARRAY(ijfrst, nsmax) - CALL PUSHINTEGER4ARRAY(nvs, nfmax) - CALL PUSHINTEGER4ARRAY(nvc, nfmax) CALL PUSHINTEGER4ARRAY(jfrst, nfmax) CALL PUSHINTEGER4ARRAY(nj, nfmax) CALL MAKESURF(isurf) @@ -80,14 +83,14 @@ SUBROUTINE UPDATE_SURFACES_B() END IF IF (ldupl(isurf)) THEN IF (lverbose) WRITE(*, *) ' reduplicating ', isurf - CALL PUSHREAL8ARRAY(vrefl, nsmax*ndmax) - CALL PUSHINTEGER4ARRAY(nvstrp, nsmax) - CALL PUSHINTEGER4ARRAY(ijfrst, nsmax) CALL PUSHBOOLEANARRAY(lmeshflat, nfmax) CALL PUSHBOOLEANARRAY(lsurfmsh, nfmax) CALL PUSHINTEGER4ARRAY(nvs, nfmax) CALL PUSHINTEGER4ARRAY(nvc, nfmax) CALL PUSHBOOLEANARRAY(lsurfspacing, nfmax) + CALL PUSHREAL8ARRAY(vrefl, nsmax*ndmax) + CALL PUSHINTEGER4ARRAY(nvstrp, nsmax) + CALL PUSHINTEGER4ARRAY(ijfrst, nsmax) CALL PUSHINTEGER4ARRAY(jfrst, nfmax) CALL PUSHINTEGER4ARRAY(nk, nfmax) CALL PUSHINTEGER4ARRAY(nj, nfmax) @@ -137,18 +140,25 @@ SUBROUTINE UPDATE_SURFACES_B() ENDDO DO ii1=1,NSURF DO ii2=1,NSEC(ii1) - DO ii3=1,ibx + DO ii3=1,nasmax xasec_diff(ii3, ii2, ii1) = 0.D0 ENDDO ENDDO ENDDO DO ii1=1,NSURF DO ii2=1,NSEC(ii1) - DO ii3=1,ibx + DO ii3=1,nasmax sasec_diff(ii3, ii2, ii1) = 0.D0 ENDDO ENDDO ENDDO + DO ii1=1,NSURF + DO ii2=1,NSEC(ii1) + DO ii3=1,6 + clcdsec_diff(ii3, ii2, ii1) = 0.D0 + ENDDO + ENDDO + ENDDO DO ii1=1,NSURF DO ii2=1,NSEC(ii1) claf_diff(ii2, ii1) = 0.D0 @@ -170,22 +180,20 @@ SUBROUTINE UPDATE_SURFACES_B() CALL POPINTEGER4ARRAY(nj, nfmax) CALL POPINTEGER4ARRAY(nk, nfmax) CALL POPINTEGER4ARRAY(jfrst, nfmax) + CALL POPINTEGER4ARRAY(ijfrst, nsmax) + CALL POPINTEGER4ARRAY(nvstrp, nsmax) + CALL POPREAL8ARRAY(vrefl, nsmax*ndmax) CALL POPBOOLEANARRAY(lsurfspacing, nfmax) CALL POPINTEGER4ARRAY(nvc, nfmax) CALL POPINTEGER4ARRAY(nvs, nfmax) CALL POPBOOLEANARRAY(lsurfmsh, nfmax) CALL POPBOOLEANARRAY(lmeshflat, nfmax) - CALL POPINTEGER4ARRAY(ijfrst, nsmax) - CALL POPINTEGER4ARRAY(nvstrp, nsmax) - CALL POPREAL8ARRAY(vrefl, nsmax*ndmax) CALL SDUPL_B(isurf, ydupl(isurf), 'ydup') END IF CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN CALL POPINTEGER4ARRAY(nj, nfmax) CALL POPINTEGER4ARRAY(jfrst, nfmax) - CALL POPINTEGER4ARRAY(nvc, nfmax) - CALL POPINTEGER4ARRAY(nvs, nfmax) CALL POPINTEGER4ARRAY(ijfrst, nsmax) CALL POPINTEGER4ARRAY(nvstrp, nsmax) CALL POPREAL8ARRAY(rle, 3*nsmax) @@ -197,15 +205,17 @@ SUBROUTINE UPDATE_SURFACES_B() CALL POPREAL8ARRAY(rv, 3*nvmax) CALL POPREAL8ARRAY(rc, 3*nvmax) CALL POPREAL8ARRAY(dxv, nvmax) + CALL POPINTEGER4ARRAY(nvc, nfmax) + CALL POPINTEGER4ARRAY(nvs, nfmax) CALL MAKESURF_MESH_B(isurf) ELSE CALL POPINTEGER4ARRAY(nj, nfmax) CALL POPINTEGER4ARRAY(jfrst, nfmax) - CALL POPINTEGER4ARRAY(nvc, nfmax) - CALL POPINTEGER4ARRAY(nvs, nfmax) CALL POPINTEGER4ARRAY(ijfrst, nsmax) CALL POPINTEGER4ARRAY(nvstrp, nsmax) CALL POPREAL8ARRAY(chord, nsmax) + CALL POPINTEGER4ARRAY(nvc, nfmax) + CALL POPINTEGER4ARRAY(nvs, nfmax) CALL MAKESURF_B(isurf) END IF ENDDO @@ -215,14 +225,16 @@ SUBROUTINE UPDATE_SURFACES_B() END C Differentiation of makesurf in reverse (adjoint) mode (with options i4 dr8 r8): -C gradient of useful results: xyzscal xyztran addinc xyzles -C chords aincs xasec sasec claf rle chord rle1 chord1 -C rle2 chord2 wstrip ainc ainc_g rv1 rv2 rv rc rs -C dxv chordv slopev slopec dcontrol vhinge -C with respect to varying inputs: xyzscal xyztran addinc xyzles -C chords aincs xasec sasec claf rle chord rle1 chord1 -C rle2 chord2 wstrip ainc ainc_g rv1 rv2 rv rc rs -C dxv chordv slopev slopec dcontrol vhinge +C gradient of useful results: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs +C dxv chordv slopev slopec dcontrol vhinge xyzscal +C xyztran addinc xyzles chords aincs xasec sasec +C clcdsec claf +C with respect to varying inputs: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs +C dxv chordv slopev slopec dcontrol vhinge xyzscal +C xyztran addinc xyzles chords aincs xasec sasec +C clcdsec claf C*********************************************************************** C Module: amake.f C @@ -246,6 +258,9 @@ SUBROUTINE UPDATE_SURFACES_B() SUBROUTINE MAKESURF_B(isurf) INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' + INCLUDE 'AVL_oml.INC' C C REAL xyzlel(3), xyzler(3) @@ -275,6 +290,10 @@ SUBROUTINE MAKESURF_B(isurf) REAL xled(ndmax), xted(ndmax), gainda(ndmax) REAL xled_diff(ndmax), xted_diff(ndmax), gainda_diff(ndmax) INTEGER idx_vor, idx_strip + REAL xlasecl(nasmax), xuasecl(nasmax) + REAL zlasecl(nasmax), zuasecl(nasmax) + REAL xlasecr(nasmax), xuasecr(nasmax) + REAL zlasecr(nasmax), zuasecr(nasmax) INTEGER isec REAL dy REAL dy_diff @@ -408,14 +427,11 @@ SUBROUTINE MAKESURF_B(isurf) REAL temp_diff0 REAL temp0 REAL temp_diff1 - INTEGER ii1 INTEGER ad_to INTEGER ad_to0 INTEGER ad_count INTEGER i INTEGER branch - INTEGER ii3 - INTEGER ii2 INTEGER ad_to1 INTEGER ad_to2 INTEGER ad_to3 @@ -1014,6 +1030,8 @@ SUBROUTINE MAKESURF_B(isurf) C C C +C +C Cc#endif CALL PUSHINTEGER4(idx_vor) idx_vor = idx_vor + 1 @@ -1024,42 +1042,18 @@ SUBROUTINE MAKESURF_B(isurf) ENDDO CALL PUSHINTEGER4(ispan - 1) ENDDO - DO ii1=1,3 - xyzler_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,kcmax - xcp_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chcosl_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chsinr_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - xted_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - xled_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chsinl_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - gainda_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - ycp_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - ypt_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - chcosr_g_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,3 - xyzlel_diff(ii1) = 0.D0 - ENDDO + xyzler_diff = 0.D0 + xcp_diff = 0.D0 + chcosl_g_diff = 0.D0 + chsinr_g_diff = 0.D0 + xted_diff = 0.D0 + xled_diff = 0.D0 + chsinl_g_diff = 0.D0 + gainda_diff = 0.D0 + ycp_diff = 0.D0 + ypt_diff = 0.D0 + chcosr_g_diff = 0.D0 + xyzlel_diff = 0.D0 DO isec=nsec(isurf)-1,1,-1 clafr = claf(isec+1, isurf) chordl = xyzscal(1, isurf)*chords(isec, isurf) @@ -1142,17 +1136,13 @@ SUBROUTINE MAKESURF_B(isurf) slopev_diff(idx_vor) = 0.D0 sloper_diff = temp_diff0 CALL POPREAL8(sloper) - DO ii1=1,kcmax - xvr_diff(ii1) = 0.D0 - ENDDO + xvr_diff = 0.D0 CALL AKIMA_B(xasec(1, isec+1, isurf), xasec_diff(1, isec + +1, isurf), sasec(1, isec+1, isurf), + sasec_diff(1, isec+1, isurf), nsr, xvr(ivc) + , xvr_diff(ivc), sloper, sloper_diff, dsdx) CALL POPREAL8(slopel) - DO ii1=1,kcmax - xvr_diff(ii1) = 0.D0 - ENDDO + xvr_diff = 0.D0 CALL AKIMA_B(xasec(1, isec, isurf), xasec_diff(1, isec, + isurf), sasec(1, isec, isurf), sasec_diff(1 + , isec, isurf), nsl, xvr(ivc), xvr_diff(ivc @@ -1246,6 +1236,15 @@ SUBROUTINE MAKESURF_B(isurf) CALL POPREAL8(chordc) chord_diff(idx_strip) = chord_diff(idx_strip) + + chordc_diff + DO l=6,1,-1 + fc_diff = fc_diff + (clcdsec(l, isec+1, isurf)-clcdsec(l + + , isec, isurf))*clcd_diff(l, idx_strip) + clcdsec_diff(l, isec, isurf) = clcdsec_diff(l, isec, + + isurf) + (1.0-fc)*clcd_diff(l, idx_strip) + clcdsec_diff(l, isec+1, isurf) = clcdsec_diff(l, isec+1 + + , isurf) + fc*clcd_diff(l, idx_strip) + clcd_diff(l, idx_strip) = 0.D0 + ENDDO DO n=ncontrol,1,-1 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN @@ -1632,153 +1631,43 @@ SUBROUTINE MAKESURF_B(isurf) IF (i .EQ. 1) THEN CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN - DO ii1=1,ksmax - yzlen_diff(ii1) = 0.D0 - ENDDO + yzlen_diff = 0.D0 GOTO 140 ELSE - DO ii1=1,NSURF - DO ii2=1,3 - xyzscal_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,3 - xyztran_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - addinc_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,3 - xyzles_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - chords_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - aincs_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,ibx - xasec_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,ibx - sasec_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - claf_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle1_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord1_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle2_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord2_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - wstrip_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - ainc_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - DO ii2=1,NSTRIP - ainc_g_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv1_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv2_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rc_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rs_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - dxv_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - chordv_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - slopev_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - slopec_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - DO ii2=1,nvor - dcontrol_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,ndmax - DO ii2=1,NSTRIP - DO ii3=1,3 - vhinge_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,ksmax - ycp_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - ypt_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - yzlen_diff(ii1) = 0.D0 - ENDDO + rle_diff = 0.D0 + chord_diff = 0.D0 + rle1_diff = 0.D0 + chord1_diff = 0.D0 + rle2_diff = 0.D0 + chord2_diff = 0.D0 + wstrip_diff = 0.D0 + clcd_diff = 0.D0 + ainc_diff = 0.D0 + ainc_g_diff = 0.D0 + rv1_diff = 0.D0 + rv2_diff = 0.D0 + rv_diff = 0.D0 + rc_diff = 0.D0 + rs_diff = 0.D0 + dxv_diff = 0.D0 + chordv_diff = 0.D0 + slopev_diff = 0.D0 + slopec_diff = 0.D0 + dcontrol_diff = 0.D0 + vhinge_diff = 0.D0 + xyzscal_diff = 0.D0 + xyztran_diff = 0.D0 + addinc_diff = 0.D0 + xyzles_diff = 0.D0 + chords_diff = 0.D0 + aincs_diff = 0.D0 + xasec_diff = 0.D0 + sasec_diff = 0.D0 + clcdsec_diff = 0.D0 + claf_diff = 0.D0 + ycp_diff = 0.D0 + ypt_diff = 0.D0 + yzlen_diff = 0.D0 dyzlen_diff = 0.D0 END IF ELSE @@ -1808,297 +1697,79 @@ SUBROUTINE MAKESURF_B(isurf) IF (i0 .EQ. 1) THEN CALL POPCONTROL2B(branch) IF (branch .EQ. 0) THEN - DO ii1=1,ksmax - yzlen_diff(ii1) = 0.D0 - ENDDO + yzlen_diff = 0.D0 GOTO 160 ELSE IF (branch .EQ. 1) THEN - DO ii1=1,NSURF - DO ii2=1,3 - xyzscal_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,3 - xyztran_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - addinc_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,3 - xyzles_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - chords_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - aincs_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,ibx - xasec_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,ibx - sasec_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - claf_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle1_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord1_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle2_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord2_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - wstrip_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - ainc_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - DO ii2=1,NSTRIP - ainc_g_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv1_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv2_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rc_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rs_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - dxv_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - chordv_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - slopev_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - slopec_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - DO ii2=1,nvor - dcontrol_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,ndmax - DO ii2=1,NSTRIP - DO ii3=1,3 - vhinge_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,ksmax - ycp_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - ypt_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - yzlen_diff(ii1) = 0.D0 - ENDDO + rle_diff = 0.D0 + chord_diff = 0.D0 + rle1_diff = 0.D0 + chord1_diff = 0.D0 + rle2_diff = 0.D0 + chord2_diff = 0.D0 + wstrip_diff = 0.D0 + clcd_diff = 0.D0 + ainc_diff = 0.D0 + ainc_g_diff = 0.D0 + rv1_diff = 0.D0 + rv2_diff = 0.D0 + rv_diff = 0.D0 + rc_diff = 0.D0 + rs_diff = 0.D0 + dxv_diff = 0.D0 + chordv_diff = 0.D0 + slopev_diff = 0.D0 + slopec_diff = 0.D0 + dcontrol_diff = 0.D0 + vhinge_diff = 0.D0 + xyzscal_diff = 0.D0 + xyztran_diff = 0.D0 + addinc_diff = 0.D0 + xyzles_diff = 0.D0 + chords_diff = 0.D0 + aincs_diff = 0.D0 + xasec_diff = 0.D0 + sasec_diff = 0.D0 + clcdsec_diff = 0.D0 + claf_diff = 0.D0 + ycp_diff = 0.D0 + ypt_diff = 0.D0 + yzlen_diff = 0.D0 GOTO 150 ELSE - DO ii1=1,NSURF - DO ii2=1,3 - xyzscal_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,3 - xyztran_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - addinc_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,3 - xyzles_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - chords_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - aincs_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,ibx - xasec_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - DO ii3=1,ibx - sasec_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,NSURF - DO ii2=1,NSEC(ii1) - claf_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle1_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord1_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - rle2_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,NSTRIP - chord2_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - wstrip_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,NSTRIP - ainc_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ngmax - DO ii2=1,NSTRIP - ainc_g_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv1_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv2_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rc_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rs_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - dxv_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - chordv_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - slopev_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,nvor - slopec_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ndmax - DO ii2=1,nvor - dcontrol_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,ndmax - DO ii2=1,NSTRIP - DO ii3=1,3 - vhinge_diff(ii3, ii2, ii1) = 0.D0 - ENDDO - ENDDO - ENDDO - DO ii1=1,ksmax - ycp_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - ypt_diff(ii1) = 0.D0 - ENDDO - DO ii1=1,ksmax - yzlen_diff(ii1) = 0.D0 - ENDDO + rle_diff = 0.D0 + chord_diff = 0.D0 + rle1_diff = 0.D0 + chord1_diff = 0.D0 + rle2_diff = 0.D0 + chord2_diff = 0.D0 + wstrip_diff = 0.D0 + clcd_diff = 0.D0 + ainc_diff = 0.D0 + ainc_g_diff = 0.D0 + rv1_diff = 0.D0 + rv2_diff = 0.D0 + rv_diff = 0.D0 + rc_diff = 0.D0 + rs_diff = 0.D0 + dxv_diff = 0.D0 + chordv_diff = 0.D0 + slopev_diff = 0.D0 + slopec_diff = 0.D0 + dcontrol_diff = 0.D0 + vhinge_diff = 0.D0 + xyzscal_diff = 0.D0 + xyztran_diff = 0.D0 + addinc_diff = 0.D0 + xyzles_diff = 0.D0 + chords_diff = 0.D0 + aincs_diff = 0.D0 + xasec_diff = 0.D0 + sasec_diff = 0.D0 + clcdsec_diff = 0.D0 + claf_diff = 0.D0 + ycp_diff = 0.D0 + ypt_diff = 0.D0 + yzlen_diff = 0.D0 END IF ELSE yscale_diff = 0.D0 @@ -2224,20 +1895,23 @@ SUBROUTINE MAKESURF_B(isurf) END C Differentiation of makesurf_mesh in reverse (adjoint) mode (with options i4 dr8 r8): -C gradient of useful results: xyzscal xyztran addinc aincs -C xasec sasec claf mshblk rv1msh rv2msh rvmsh rcmsh -C rle chord rle1 chord1 rle2 chord2 wstrip ainc -C ainc_g rv1 rv2 rv rc rs dxv chordv slopev slopec -C dcontrol vhinge -C with respect to varying inputs: xyzscal xyztran addinc aincs -C xasec sasec claf mshblk rv1msh rv2msh rvmsh rcmsh -C rle chord rle1 chord1 rle2 chord2 wstrip ainc -C ainc_g rv1 rv2 rv rc rs dxv chordv slopev slopec -C dcontrol vhinge +C gradient of useful results: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs +C dxv chordv slopev slopec dcontrol vhinge xyzscal +C xyztran addinc aincs xasec sasec clcdsec claf +C mshblk rv1msh rv2msh rvmsh rcmsh +C with respect to varying inputs: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc rs +C dxv chordv slopev slopec dcontrol vhinge xyzscal +C xyztran addinc aincs xasec sasec clcdsec claf +C mshblk rv1msh rv2msh rvmsh rcmsh C SUBROUTINE MAKESURF_MESH_B(isurf) INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' + INCLUDE 'AVL_oml.INC' C working variables (AVL original) INTEGER isurf INTEGER kcmax @@ -2254,8 +1928,12 @@ SUBROUTINE MAKESURF_MESH_B(isurf) REAL chsinl_g_diff(ngmax), chcosl_g_diff(ngmax), chsinr_g_diff( + ngmax), chcosr_g_diff(ngmax), xled_diff(ndmax), xted_diff( + ndmax) -C working variables (OptVL additions) INTEGER isconl(ndmax), isconr(ndmax) + REAL xlasecl(nasmax), xuasecl(nasmax) + REAL zlasecl(nasmax), zuasecl(nasmax) + REAL xlasecr(nasmax), xuasecr(nasmax) +C working variables (OptVL additions) + REAL zlasecr(nasmax), zuasecr(nasmax) REAL m1, m2, m3, f1, f2, fc, dc1, dc2, dc, a1, a2, a3, xptxind1, + xptxind2 REAL m2_diff, m3_diff, dc1_diff, dc2_diff, a1_diff, a2_diff, @@ -2972,7 +2650,10 @@ SUBROUTINE MAKESURF_MESH_B(isurf) C C xptxind2 = (mesh_surf(1,idx_node_yp1+1) C & - RLE2(1,idx_strip))/CHORD2(idx_strip) -C Interpolate cross section on left side +C +C +C ! Interpolate cross section on left side +C C C C Interpolate cross section on right side @@ -3571,6 +3252,13 @@ SUBROUTINE MAKESURF_MESH_B(isurf) + idx_strip) chordl_diff = chordl_diff + clafl*temp_diff2 clafl_diff = chordl*temp_diff2 + DO idx_coef=6,1,-1 + clcdsec_diff(idx_coef, iptl, isurf) = clcdsec_diff(idx_coef + + , iptl, isurf) + (1.0-fc)*clcd_diff(idx_coef, idx_strip) + clcdsec_diff(idx_coef, iptr, isurf) = clcdsec_diff(idx_coef + + , iptr, isurf) + fc*clcd_diff(idx_coef, idx_strip) + clcd_diff(idx_coef, idx_strip) = 0.D0 + ENDDO DO n=ncontrol,1,-1 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN @@ -3938,18 +3626,21 @@ SUBROUTINE MAKESURF_MESH_B(isurf) END C Differentiation of sdupl in reverse (adjoint) mode (with options i4 dr8 r8): -C gradient of useful results: rv1msh rv2msh rvmsh rcmsh rle -C chord rle1 chord1 rle2 chord2 wstrip ainc ainc_g -C rv1 rv2 rv rc dxv chordv slopev slopec dcontrol -C vhinge -C with respect to varying inputs: rv1msh rv2msh rvmsh rcmsh rle -C chord rle1 chord1 rle2 chord2 wstrip ainc ainc_g -C rv1 rv2 rv rc dxv chordv slopev slopec dcontrol -C vhinge +C gradient of useful results: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc dxv +C chordv slopev slopec dcontrol vhinge rv1msh rv2msh +C rvmsh rcmsh +C with respect to varying inputs: rle chord rle1 chord1 rle2 +C chord2 wstrip clcd ainc ainc_g rv1 rv2 rv rc dxv +C chordv slopev slopec dcontrol vhinge rv1msh rv2msh +C rvmsh rcmsh C SUBROUTINE SDUPL_B(nn, ypt, msg) INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' + INCLUDE 'AVL_oml.INC' CHARACTER*(*) msg INTEGER idx_vor INTEGER nni @@ -3999,41 +3690,42 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) REAL(kind=avl_real) tmp15 REAL(kind=avl_real) tmp16 REAL(kind=avl_real) tmp17 - REAL(kind=avl_real) tmp18 REAL(kind=avl_real) tmp_diff9 - REAL(kind=avl_real) tmp19 + REAL(kind=avl_real) tmp18 REAL(kind=avl_real) tmp_diff10 - REAL(kind=avl_real) tmp20 + REAL(kind=avl_real) tmp19 REAL(kind=avl_real) tmp_diff11 - REAL(kind=avl_real) tmp21 + REAL(kind=avl_real) tmp20 REAL(kind=avl_real) tmp_diff12 - REAL(kind=avl_real) tmp22 + REAL(kind=avl_real) tmp21 REAL(kind=avl_real) tmp_diff13 - REAL(kind=avl_real) tmp23 + REAL(kind=avl_real) tmp22 REAL(kind=avl_real) tmp_diff14 - REAL(kind=avl_real) tmp24 + REAL(kind=avl_real) tmp23 REAL(kind=avl_real) tmp_diff15 - REAL(kind=avl_real) tmp25 + REAL(kind=avl_real) tmp24 REAL(kind=avl_real) tmp_diff16 - REAL(kind=avl_real) tmp26 + REAL(kind=avl_real) tmp25 REAL(kind=avl_real) tmp_diff17 + REAL(kind=avl_real) tmp26 + REAL(kind=avl_real) tmp_diff18 REAL(kind=avl_real) tmp27 REAL(kind=avl_real) tmp28 - REAL(kind=avl_real) tmp_diff18 - REAL(kind=avl_real) tmp29 REAL(kind=avl_real) tmp_diff19 - REAL(kind=avl_real) tmp30 + REAL(kind=avl_real) tmp29 REAL(kind=avl_real) tmp_diff20 - REAL(kind=avl_real) tmp31 + REAL(kind=avl_real) tmp30 REAL(kind=avl_real) tmp_diff21 - REAL(kind=avl_real) tmp32 + REAL(kind=avl_real) tmp31 REAL(kind=avl_real) tmp_diff22 - REAL(kind=avl_real) tmp33 + REAL(kind=avl_real) tmp32 REAL(kind=avl_real) tmp_diff23 - REAL(kind=avl_real) tmp34 + REAL(kind=avl_real) tmp33 REAL(kind=avl_real) tmp_diff24 - REAL(kind=avl_real) tmp35 + REAL(kind=avl_real) tmp34 REAL(kind=avl_real) tmp_diff25 + REAL(kind=avl_real) tmp35 + REAL(kind=avl_real) tmp_diff26 INTEGER ad_count INTEGER i INTEGER branch @@ -4042,9 +3734,9 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) INTEGER ad_to1 INTEGER ad_count0 INTEGER i0 - INTEGER ii3 INTEGER ii2 INTEGER ii1 + INTEGER ii3 INTEGER nn REAL ypt C @@ -4132,6 +3824,7 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) ijfrst(jji) = ijfrst(jji-1) + nvstrp(jji-1) C nvstrp(jji) = nvc(nni) + l = 7 C idx_vor = ijfrst(jji) ad_count0 = 1 @@ -4186,26 +3879,6 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) IF (i0 .EQ. 1) THEN CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN - DO ii1=1,nvor - DO ii2=1,3 - rv1msh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv2msh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rvmsh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rcmsh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO DO ii1=1,NSTRIP DO ii2=1,3 rle_diff(ii2, ii1) = 0.D0 @@ -4233,6 +3906,11 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) DO ii1=1,NSTRIP wstrip_diff(ii1) = 0.D0 ENDDO + DO ii1=1,NSTRIP + DO ii2=1,6 + clcd_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO DO ii1=1,NSTRIP ainc_diff(ii1) = 0.D0 ENDDO @@ -4285,36 +3963,56 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) ENDDO ENDDO ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rv1msh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rv2msh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rvmsh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rcmsh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO END IF ELSE CALL POPINTEGER4(ad_to1) DO n=ad_to1,1,-1 - tmp_diff25 = dcontrol_diff(iii, n) + tmp_diff26 = dcontrol_diff(iii, n) dcontrol_diff(iii, n) = 0.D0 dcontrol_diff(ii, n) = dcontrol_diff(ii, n) - rsgn* - + tmp_diff25 + + tmp_diff26 CALL POPREAL8(rsgn) ENDDO CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN - tmp_diff24 = rcmsh_diff(3, iii) + tmp_diff25 = rcmsh_diff(3, iii) rcmsh_diff(3, iii) = 0.D0 - rcmsh_diff(3, ii) = rcmsh_diff(3, ii) + tmp_diff24 - tmp_diff23 = rcmsh_diff(2, iii) + rcmsh_diff(3, ii) = rcmsh_diff(3, ii) + tmp_diff25 + tmp_diff24 = rcmsh_diff(2, iii) rcmsh_diff(2, iii) = 0.D0 - rcmsh_diff(2, ii) = rcmsh_diff(2, ii) - tmp_diff23 - tmp_diff22 = rcmsh_diff(1, iii) + rcmsh_diff(2, ii) = rcmsh_diff(2, ii) - tmp_diff24 + tmp_diff23 = rcmsh_diff(1, iii) rcmsh_diff(1, iii) = 0.D0 - rcmsh_diff(1, ii) = rcmsh_diff(1, ii) + tmp_diff22 - tmp_diff21 = rvmsh_diff(3, iii) + rcmsh_diff(1, ii) = rcmsh_diff(1, ii) + tmp_diff23 + tmp_diff22 = rvmsh_diff(3, iii) rvmsh_diff(3, iii) = 0.D0 - rvmsh_diff(3, ii) = rvmsh_diff(3, ii) + tmp_diff21 - tmp_diff20 = rvmsh_diff(2, iii) + rvmsh_diff(3, ii) = rvmsh_diff(3, ii) + tmp_diff22 + tmp_diff21 = rvmsh_diff(2, iii) rvmsh_diff(2, iii) = 0.D0 - rvmsh_diff(2, ii) = rvmsh_diff(2, ii) - tmp_diff20 - tmp_diff19 = rvmsh_diff(1, iii) + rvmsh_diff(2, ii) = rvmsh_diff(2, ii) - tmp_diff21 + tmp_diff20 = rvmsh_diff(1, iii) rvmsh_diff(1, iii) = 0.D0 - rvmsh_diff(1, ii) = rvmsh_diff(1, ii) + tmp_diff19 + rvmsh_diff(1, ii) = rvmsh_diff(1, ii) + tmp_diff20 rv1msh_diff(3, ii) = rv1msh_diff(3, ii) + rv2msh_diff(3 + , iii) rv2msh_diff(3, iii) = 0.D0 @@ -4334,36 +4032,36 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) + , iii) rv1msh_diff(1, iii) = 0.D0 END IF - tmp_diff18 = chordv_diff(iii) + tmp_diff19 = chordv_diff(iii) chordv_diff(iii) = 0.D0 - chordv_diff(ii) = chordv_diff(ii) + tmp_diff18 - tmp_diff17 = dxv_diff(iii) + chordv_diff(ii) = chordv_diff(ii) + tmp_diff19 + tmp_diff18 = dxv_diff(iii) dxv_diff(iii) = 0.D0 - dxv_diff(ii) = dxv_diff(ii) + tmp_diff17 - tmp_diff16 = slopev_diff(iii) + dxv_diff(ii) = dxv_diff(ii) + tmp_diff18 + tmp_diff17 = slopev_diff(iii) slopev_diff(iii) = 0.D0 - slopev_diff(ii) = slopev_diff(ii) + tmp_diff16 - tmp_diff15 = slopec_diff(iii) + slopev_diff(ii) = slopev_diff(ii) + tmp_diff17 + tmp_diff16 = slopec_diff(iii) slopec_diff(iii) = 0.D0 - slopec_diff(ii) = slopec_diff(ii) + tmp_diff15 - tmp_diff14 = rc_diff(3, iii) + slopec_diff(ii) = slopec_diff(ii) + tmp_diff16 + tmp_diff15 = rc_diff(3, iii) rc_diff(3, iii) = 0.D0 - rc_diff(3, ii) = rc_diff(3, ii) + tmp_diff14 - tmp_diff13 = rc_diff(2, iii) + rc_diff(3, ii) = rc_diff(3, ii) + tmp_diff15 + tmp_diff14 = rc_diff(2, iii) rc_diff(2, iii) = 0.D0 - rc_diff(2, ii) = rc_diff(2, ii) - tmp_diff13 - tmp_diff12 = rc_diff(1, iii) + rc_diff(2, ii) = rc_diff(2, ii) - tmp_diff14 + tmp_diff13 = rc_diff(1, iii) rc_diff(1, iii) = 0.D0 - rc_diff(1, ii) = rc_diff(1, ii) + tmp_diff12 - tmp_diff11 = rv_diff(3, iii) + rc_diff(1, ii) = rc_diff(1, ii) + tmp_diff13 + tmp_diff12 = rv_diff(3, iii) rv_diff(3, iii) = 0.D0 - rv_diff(3, ii) = rv_diff(3, ii) + tmp_diff11 - tmp_diff10 = rv_diff(2, iii) + rv_diff(3, ii) = rv_diff(3, ii) + tmp_diff12 + tmp_diff11 = rv_diff(2, iii) rv_diff(2, iii) = 0.D0 - rv_diff(2, ii) = rv_diff(2, ii) - tmp_diff10 - tmp_diff9 = rv_diff(1, iii) + rv_diff(2, ii) = rv_diff(2, ii) - tmp_diff11 + tmp_diff10 = rv_diff(1, iii) rv_diff(1, iii) = 0.D0 - rv_diff(1, ii) = rv_diff(1, ii) + tmp_diff9 + rv_diff(1, ii) = rv_diff(1, ii) + tmp_diff10 rv1_diff(3, ii) = rv1_diff(3, ii) + rv2_diff(3, iii) rv2_diff(3, iii) = 0.D0 rv1_diff(2, ii) = rv1_diff(2, ii) - rv2_diff(2, iii) @@ -4380,6 +4078,11 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) CALL POPINTEGER4(iii) END IF ENDDO + DO l=6,1,-1 + tmp_diff9 = clcd_diff(l, jji) + clcd_diff(l, jji) = 0.D0 + clcd_diff(l, jj) = clcd_diff(l, jj) + tmp_diff9 + ENDDO CALL POPINTEGER4(ad_to0) DO n=ad_to0,1,-1 tmp_diff8 = vhinge_diff(3, jji, n) @@ -4451,9 +4154,9 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) C Differentiation of encalc in reverse (adjoint) mode (with options i4 dr8 r8): C gradient of useful results: ess ensy ensz xsref ysref zsref C rv1 rv2 rv enc env enc_d -C with respect to varying inputs: rv1msh rv2msh rvmsh rcmsh ess -C ensy ensz xsref ysref zsref ainc ainc_g rv1 rv2 -C rv slopev slopec dcontrol vhinge enc env enc_d +C with respect to varying inputs: ess ensy ensz xsref ysref zsref +C ainc ainc_g rv1 rv2 rv slopev slopec dcontrol +C vhinge enc env enc_d rv1msh rv2msh rvmsh rcmsh C BDUPL C C @@ -4464,6 +4167,8 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) SUBROUTINE ENCALC_B() INCLUDE 'AVL.INC' INCLUDE 'AVL_ad_seeds.inc' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_surf_ad_seeds.inc' C REAL ep(3), eq(3), es(3), eb(3), ec(3), ecxb(3) REAL ep_diff(3), eq_diff(3), es_diff(3), eb_diff(3), ec_diff(3), @@ -4541,8 +4246,8 @@ SUBROUTINE ENCALC_B() INTEGER branch INTEGER ad_to INTEGER ii1 - INTEGER ii3 INTEGER ii2 + INTEGER ii3 C C...Calculate the normal vector at control points and bound vortex midpoints C @@ -4909,26 +4614,6 @@ SUBROUTINE ENCALC_B() ENDDO CALL PUSHINTEGER4(ii - 1) ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv1msh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rv2msh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rvmsh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO - DO ii1=1,nvor - DO ii2=1,3 - rcmsh_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO DO ii1=1,NSTRIP ainc_diff(ii1) = 0.D0 ENDDO @@ -4955,6 +4640,26 @@ SUBROUTINE ENCALC_B() ENDDO ENDDO ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rv1msh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rv2msh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rvmsh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,nvor + DO ii2=1,3 + rcmsh_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO DO ii1=1,3 eb_diff(ii1) = 0.D0 ENDDO diff --git a/src/ad_src/reverse_ad_src/cdcl_b.f b/src/ad_src/reverse_ad_src/cdcl_b.f index 453e6da6..5860d495 100644 --- a/src/ad_src/reverse_ad_src/cdcl_b.f +++ b/src/ad_src/reverse_ad_src/cdcl_b.f @@ -2,8 +2,8 @@ C Tapenade 3.16 (develop) - 15 Jan 2021 14:26 C C Differentiation of cdcl in reverse (adjoint) mode (with options i4 dr8 r8): -C gradient of useful results: cd_cl cd -C with respect to varying inputs: cl +C gradient of useful results: cd_cl cd cdclpol +C with respect to varying inputs: cl cdclpol C*********************************************************************** C Module: cdcl.f C @@ -24,8 +24,8 @@ C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. C*********************************************************************** C - SUBROUTINE CDCL_B(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, - + cd_cl_diff) + SUBROUTINE CDCL_B(cdclpol, cdclpol_diff, cl, cl_diff, cd, cd_diff + + , cd_cl, cd_cl_diff) C------------------------------------------------------------------------- C Returns cd as a function of cl, or drag polar. C @@ -59,15 +59,31 @@ SUBROUTINE CDCL_B(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, C cdclpol(6) cd (cdpos) at clpos C------------------------------------------------------------------------- REAL cdclpol(6) + REAL cdclpol_diff(6) REAL clmin + REAL clmin_diff REAL cdmin + REAL cdmin_diff REAL cl0 + REAL cl0_diff REAL cd0 + REAL cd0_diff REAL clmax + REAL clmax_diff REAL cdmax + REAL cdmax_diff REAL cdx1 + REAL cdx1_diff REAL cdx2 + REAL cdx2_diff REAL clfac + REAL temp + REAL temp_diff + REAL temp0 + REAL temp_diff0 + REAL temp1 + REAL temp_diff1 + REAL temp_diff2 REAL clinc REAL cd_cl REAL cd_cl_diff @@ -92,6 +108,12 @@ SUBROUTINE CDCL_B(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, cdmax = cdclpol(6) IF (clmax .LE. cl0 .OR. cl0 .LE. clmin) THEN cl_diff = 0.D0 + clmin_diff = 0.D0 + cl0_diff = 0.D0 + cdmax_diff = 0.D0 + clmax_diff = 0.D0 + cd0_diff = 0.D0 + cdmin_diff = 0.D0 ELSE C C---- matching parameters that make the slopes smooth at the stall joins @@ -100,19 +122,103 @@ SUBROUTINE CDCL_B(cdclpol, cl, cl_diff, cd, cd_diff, cd_cl, clfac = 1.0/clinc C IF (cl .LT. clmin) THEN - cl_diff = cdinc*2.0*clfac**2*cd_cl_diff + (2*(cl-clmin)*cdinc* - + clfac**2-cdx1/(clmin-cl0))*cd_diff - ELSE IF (cl .LT. cl0) THEN - cl_diff = (cdmin-cd0)*2.0*cd_cl_diff/(clmin-cl0)**2 + 2*(cl- - + cl0)*(cdmin-cd0)*cd_diff/(clmin-cl0)**2 - ELSE IF (cl .LT. clmax) THEN - cl_diff = (cdmax-cd0)*2.0*cd_cl_diff/(clmax-cl0)**2 + 2*(cl- - + cl0)*(cdmax-cd0)*cd_diff/(clmax-cl0)**2 + temp = (cl-cl0)/(clmin-cl0) + temp_diff1 = cdinc*2.0*clfac**2*cd_cl_diff + temp_diff0 = -(cd_cl_diff/(clmin-cl0)) + cdx1_diff = temp_diff0 + (1.0-temp)*cd_diff + temp_diff = -(cdx1*temp_diff0/(clmin-cl0)) + clmin_diff = temp_diff - temp_diff1 + cl0_diff = -temp_diff + cdmin_diff = cd_diff + temp_diff0 = 2*(cl-clmin)*cdinc*clfac**2*cd_diff + temp_diff = -(cdx1*cd_diff/(clmin-cl0)) + cl_diff = temp_diff1 + temp_diff + temp_diff0 + temp_diff1 = -(temp*temp_diff) + cl0_diff = cl0_diff - temp_diff - temp_diff1 + clmin_diff = clmin_diff + temp_diff1 - temp_diff0 + cdmax_diff = 0.D0 + cdx2_diff = 0.D0 + clmax_diff = 0.D0 + cd0_diff = 0.D0 ELSE - cl_diff = cdinc*2.0*clfac**2*cd_cl_diff + (2*(cl-clmax)*cdinc* - + clfac**2+cdx2/(clmax-cl0))*cd_diff + IF (cl .LT. cl0) THEN + temp_diff2 = 2.0*cd_cl_diff/(clmin-cl0)**2 + temp_diff1 = -(2*(cdmin-cd0)*(cl-cl0)*temp_diff2/(clmin-cl0) + + ) + cl0_diff = -((cdmin-cd0)*temp_diff2) - temp_diff1 + clmin_diff = temp_diff1 + temp0 = (cl-cl0)*(cl-cl0) + temp_diff1 = cd_diff/(clmin-cl0)**2 + cdmin_diff = (cl-cl0)*temp_diff2 + temp0*temp_diff1 + cd0_diff = cd_diff - (cl-cl0)*temp_diff2 - temp0*temp_diff1 + temp_diff = 2*(cl-cl0)*(cdmin-cd0)*temp_diff1 + cl_diff = (cdmin-cd0)*temp_diff2 + temp_diff + temp_diff2 = -(2*(cdmin-cd0)*temp0*temp_diff1/(clmin-cl0)) + clmin_diff = clmin_diff + temp_diff2 + cl0_diff = cl0_diff - temp_diff2 - temp_diff + cdmax_diff = 0.D0 + cdx2_diff = 0.D0 + clmax_diff = 0.D0 + ELSE + IF (cl .LT. clmax) THEN + temp_diff2 = 2.0*cd_cl_diff/(clmax-cl0)**2 + cdmax_diff = (cl-cl0)*temp_diff2 + cd0_diff = -((cl-cl0)*temp_diff2) + cl_diff = (cdmax-cd0)*temp_diff2 + cl0_diff = -((cdmax-cd0)*temp_diff2) + temp_diff1 = -(2*(cdmax-cd0)*(cl-cl0)*temp_diff2/(clmax- + + cl0)) + temp1 = (cl-cl0)*(cl-cl0) + temp_diff2 = cd_diff/(clmax-cl0)**2 + cd0_diff = cd0_diff + cd_diff - temp1*temp_diff2 + cdmax_diff = cdmax_diff + temp1*temp_diff2 + temp_diff0 = 2*(cl-cl0)*(cdmax-cd0)*temp_diff2 + temp_diff = -(2*(cdmax-cd0)*temp1*temp_diff2/(clmax-cl0)) + clmax_diff = temp_diff1 + temp_diff + cl0_diff = cl0_diff - temp_diff1 - temp_diff - temp_diff0 + cl_diff = cl_diff + temp_diff0 + cdx2_diff = 0.D0 + ELSE + temp1 = (cl-cl0)/(clmax-cl0) + temp_diff2 = cdinc*2.0*clfac**2*cd_cl_diff + temp_diff1 = cd_cl_diff/(clmax-cl0) + cdx2_diff = temp_diff1 - (1.0-temp1)*cd_diff + temp_diff0 = -(cdx2*temp_diff1/(clmax-cl0)) + clmax_diff = temp_diff0 - temp_diff2 + cl0_diff = -temp_diff0 + cl_diff = temp_diff2 + cdmax_diff = cd_diff + temp_diff2 = 2*(cl-clmax)*cdinc*clfac**2*cd_diff + temp_diff1 = cdx2*cd_diff/(clmax-cl0) + cl_diff = cl_diff + temp_diff1 + temp_diff2 + temp_diff0 = -(temp1*temp_diff1) + cl0_diff = cl0_diff - temp_diff1 - temp_diff0 + clmax_diff = clmax_diff + temp_diff0 - temp_diff2 + cd0_diff = 0.D0 + END IF + clmin_diff = 0.D0 + cdmin_diff = 0.D0 + END IF + cdx1_diff = 0.D0 END IF + temp_diff0 = 2.0*cdx2_diff/(clmax-cl0) + cdmax_diff = cdmax_diff + temp_diff0 + temp_diff = -((cdmax-cd0)*temp_diff0/(clmax-cl0)) + clmax_diff = clmax_diff + temp_diff + cl0_diff = cl0_diff - temp_diff + temp_diff = 2.0*cdx1_diff/(clmin-cl0) + cd0_diff = cd0_diff - temp_diff0 - temp_diff + cdmin_diff = cdmin_diff + temp_diff + temp_diff0 = -((cdmin-cd0)*temp_diff/(clmin-cl0)) + clmin_diff = clmin_diff + temp_diff0 + cl0_diff = cl0_diff - temp_diff0 END IF + cdclpol_diff(6) = cdclpol_diff(6) + cdmax_diff + cdclpol_diff(5) = cdclpol_diff(5) + clmax_diff + cdclpol_diff(4) = cdclpol_diff(4) + cd0_diff + cdclpol_diff(3) = cdclpol_diff(3) + cl0_diff + cdclpol_diff(2) = cdclpol_diff(2) + cdmin_diff + cdclpol_diff(1) = cdclpol_diff(1) + clmin_diff END C cdcl diff --git a/src/ainput.f b/src/ainput.f index 93c42cb2..d11c3fc1 100644 --- a/src/ainput.f +++ b/src/ainput.f @@ -23,6 +23,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C Reads an processes an AVL configuration input file C--------------------------------------------------------- INCLUDE 'AVL.INC' + INCLUDE 'AVL_surf.INC' C CHARACTER*(*) FNAME LOGICAL FERR @@ -604,14 +605,14 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) TASEC(2,ISEC, ISURF) = 0.0 cc#ifdef USE_CPOML C - XLASEC(1,ISEC,ISURF) = 0 - XLASEC(2,ISEC,ISURF) = 0 - XUASEC(1,ISEC,ISURF) = 0 - XUASEC(2,ISEC,ISURF) = 0 - ZLASEC(1,ISEC,ISURF) = 0 - ZLASEC(2,ISEC,ISURF) = 0 - ZUASEC(1,ISEC,ISURF) = 0 - ZUASEC(2,ISEC,ISURF) = 0 + ! XLASEC(1,ISEC,ISURF) = 0 + ! XLASEC(2,ISEC,ISURF) = 0 + ! XUASEC(1,ISEC,ISURF) = 0 + ! XUASEC(2,ISEC,ISURF) = 0 + ! ZLASEC(1,ISEC,ISURF) = 0 + ! ZLASEC(2,ISEC,ISURF) = 0 + ! ZUASEC(1,ISEC,ISURF) = 0 + ! ZUASEC(2,ISEC,ISURF) = 0 CASEC(2,ISEC,ISURF) = 0 CASEC(2,ISEC,ISURF) = 0 cc#endif @@ -708,12 +709,9 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) ELSE ZF = 0 ENDIF - TH = ATAN(SLP) - XLASEC(I,ISEC,ISURF) = XF + 0.5*THK*SIN(TH) - XUASEC(I,ISEC,ISURF) = XF - 0.5*THK*SIN(TH) - ZLASEC(I,ISEC,ISURF) = ZF - 0.5*THK*COS(TH) - ZUASEC(I,ISEC,ISURF) = ZF + 0.5*THK*COS(TH) + ! TH = ATAN(SLP) CASEC(I,ISEC,ISURF) = ZF + cc#endif ENDDO @@ -724,15 +722,15 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) idx_upper = NASEC(ISEC,ISURF) - I +1 idx_lower = I - XSEC(I ,isec, isurf) = - & XUASEC(idx_upper, ISEC, ISURF) - XSEC(I+NASEC(ISEC,ISURF),isec, isurf) = - & XLASEC(idx_lower, ISEC, ISURF) +! XSEC(I ,isec, isurf) = +! & XUASEC(idx_upper, ISEC, ISURF) +! XSEC(I+NASEC(ISEC,ISURF),isec, isurf) = +! & XLASEC(idx_lower, ISEC, ISURF) - YSEC(I ,isec, isurf) = - & ZUASEC(idx_upper, ISEC, ISURF) - YSEC(I+NASEC(ISEC,ISURF),isec, isurf) = - & ZLASEC(idx_lower, ISEC, ISURF) +! YSEC(I ,isec, isurf) = +! & ZUASEC(idx_upper, ISEC, ISURF) +! YSEC(I+NASEC(ISEC,ISURF),isec, isurf) = +! & ZLASEC(idx_lower, ISEC, ISURF) enddo @@ -826,10 +824,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) cc#ifdef USE_CPOML C... lower/upper coordinates (for oml) CALL AKIMA(XIN,YIN,NIN,XASEC(I,ISEC,ISURF),ZC,DUMMY) - XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) - XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) - ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF) - ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF) + ! XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) + ! XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) + ! ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF) + ! ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF) CASEC(I,ISEC,ISURF) = ZC cc#endif END DO @@ -918,10 +916,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) SASEC(I,ISEC,ISURF) = 0.0 TASEC(I,ISEC,ISURF) = 0.0 cc#ifdef USE_CPOML - XLASEC(I,ISEC,ISURF) = 0 - XUASEC(I,ISEC,ISURF) = 0 - ZLASEC(I,ISEC,ISURF) = 0 - ZUASEC(I,ISEC,ISURF) = 0 + ! XLASEC(I,ISEC,ISURF) = 0 + ! XUASEC(I,ISEC,ISURF) = 0 + ! ZLASEC(I,ISEC,ISURF) = 0 + ! ZUASEC(I,ISEC,ISURF) = 0 CASEC(I,ISEC,ISURF) = 0 cc#endif ENDDO @@ -949,11 +947,9 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) cc#ifdef USE_CPOML C... lower/upper coordinates (for oml) CALL AKIMA(XIN,YIN,NIN,XASEC(I,ISEC,ISURF),ZC,DUMMY) - XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) - XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) - ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF) - ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF) CASEC(I,ISEC,ISURF) = ZC + + cc#endif END DO CALL NRMLIZ (NASEC(ISEC,ISURF),XASEC(1,ISEC,ISURF)) diff --git a/src/aio.f b/src/aio.f index 38716cea..0b24d645 100644 --- a/src/aio.f +++ b/src/aio.f @@ -1,6 +1,7 @@ subroutine write_tecplot(file_name, add_time, solution_time) C INCLUDE 'AVL.INC' + INCLUDE 'AVL_oml.INC' C LOGICAL LRANGEALL INTEGER LU, nChords, nStrips, nPts diff --git a/src/amake.f b/src/amake.f index 6571b47e..9cd6dcc3 100644 --- a/src/amake.f +++ b/src/amake.f @@ -24,6 +24,8 @@ SUBROUTINE MAKESURF(ISURF) C using info from configuration input file. C-------------------------------------------------------------- INCLUDE 'AVL.INC' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_oml.INC' C C REAL XYZLEL(3), XYZLER(3) @@ -46,6 +48,11 @@ SUBROUTINE MAKESURF(ISURF) INTEGER ISCONL(NDMAX), ISCONR(NDMAX) REAL XLED(NDMAX), XTED(NDMAX), GAINDA(NDMAX) integer idx_vor, idx_strip + REAL XLASECL(NASMAX), XUASECL(NASMAX) + REAL ZLASECL(NASMAX), ZUASECL(NASMAX) + REAL XLASECR(NASMAX), XUASECR(NASMAX) + REAL ZLASECR(NASMAX), ZUASECR(NASMAX) + C C IF(NSEC(ISURF).LT.2) THEN @@ -577,14 +584,25 @@ SUBROUTINE MAKESURF(ISURF) C... nodal grid associated with vortex strip (aft-panel nodes) C... NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral; C... retained in (x,z) plane at this point - CALL AKIMA( XLASEC(1,ISEC,ISURF), ZLASEC(1,ISEC,ISURF), NSL, + XLASECL = XASEC(:,ISEC,ISURF) + XUASECL = XASEC(:,ISEC,ISURF) + ZLASECL = CASEC(:,ISEC,ISURF) - 0.5*TASEC(:,ISEC,ISURF) + ZUASECL = CASEC(:,ISEC,ISURF) + 0.5*TASEC(:,ISEC,ISURF) + + XLASECR = XASEC(:,ISEC,ISURF) + XUASECR = XASEC(:,ISEC,ISURF) + ZLASECR = CASEC(:,ISEC,ISURF) - 0.5*TASEC(:,ISEC,ISURF) + ZUASECR = CASEC(:,ISEC,ISURF) + 0.5*TASEC(:,ISEC,ISURF) + + + CALL AKIMA( XLASECL, ZLASECL, NSL, & XPT(IVC+1), ZL_L, DSDX ) - CALL AKIMA( XUASEC(1,ISEC,ISURF), ZUASEC(1,ISEC,ISURF), NSL, + CALL AKIMA( XUASECL, ZUASECL, NSL, & XPT(IVC+1), ZU_L, DSDX ) C - CALL AKIMA( XLASEC(1,ISEC+1,ISURF), ZLASEC(1,ISEC+1,ISURF), + CALL AKIMA( XLASECR, ZLASECR, & NSR, XPT(IVC+1), ZL_R, DSDX ) - CALL AKIMA( XUASEC(1,ISEC+1,ISURF), ZUASEC(1,ISEC+1,ISURF), + CALL AKIMA( XUASECR, ZUASECR, & NSR, XPT(IVC+1), ZU_R, DSDX ) C XYN1(1,idx_vor) = RLE1(1,idx_strip) + @@ -642,6 +660,7 @@ SUBROUTINE MAKESURF(ISURF) integer function flatidx(idx_x, idx_y, idx_surf) include 'AVL.INC' + INCLUDE 'AVL_surf.INC' ! store MFRST and NVC in the common block integer idx_x, idx_y, idx_surf flatidx = idx_x + (idx_y - 1) * (NVC(idx_surf)+1) @@ -655,6 +674,8 @@ subroutine makesurf_mesh(isurf) C and the given mesh coordinate array. c-------------------------------------------------------------- INCLUDE 'AVL.INC' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_oml.INC' ! input/output integer isurf @@ -668,6 +689,10 @@ subroutine makesurf_mesh(isurf) & CHSINR_G(NGMAX),CHCOSR_G(NGMAX), & XLED(NDMAX), XTED(NDMAX), GAINDA(NDMAX) INTEGER ISCONL(NDMAX), ISCONR(NDMAX) + REAL XLASECL(NASMAX), XUASECL(NASMAX) + REAL ZLASECL(NASMAX), ZUASECL(NASMAX) + REAL XLASECR(NASMAX), XUASECR(NASMAX) + REAL ZLASECR(NASMAX), ZUASECR(NASMAX) ! working variables (OptVL additions) real m1, m2, m3, f1, f2, fc, dc1, dc2, dc, a1, a2, a3, xptxind1, @@ -1332,18 +1357,30 @@ subroutine makesurf_mesh(isurf) ! xptxind2 = (mesh_surf(1,idx_node_yp1+1) ! & - RLE2(1,idx_strip))/CHORD2(idx_strip) - ! Interpolate cross section on left side - CALL AKIMA( XLASEC(1,iptl,isurf), ZLASEC(1,iptl,isurf), +! ! Interpolate cross section on left side + + XLASECL = XASEC(:,iptl,ISURF) + XUASECL = XASEC(:,iptl,ISURF) + ZLASECL = CASEC(:,iptl,ISURF) - 0.5*TASEC(:,iptl,ISURF) + ZUASECL = CASEC(:,iptl,ISURF) + 0.5*TASEC(:,iptl,ISURF) + + XLASECR = XASEC(:,iptr,ISURF) + XUASECR = XASEC(:,iptr,ISURF) + ZLASECR = CASEC(:,iptr,ISURF) - 0.5*TASEC(:,iptr,ISURF) + ZUASECR = CASEC(:,iptr,ISURF) + 0.5*TASEC(:,iptr,ISURF) + + + CALL AKIMA( XLASECL, ZLASECL, & NSL,xptxind1, ZL_L, DSDX ) - CALL AKIMA( XUASEC(1,iptl,isurf), ZUASEC(1,iptl,isurf), + CALL AKIMA( XUASECL, ZUASECL, & NSL,xptxind1, ZU_L, DSDX ) ! Interpolate cross section on right side - CALL AKIMA(XLASEC(1,iptr,isurf), - & ZLASEC(1,iptr,isurf),NSR, xptxind1, ZL_R, DSDX) + CALL AKIMA(XLASECR, + & ZLASECR,NSR, xptxind1, ZL_R, DSDX) - CALL AKIMA(XUASEC(1,iptr,isurf), - & ZUASEC(1,iptr,isurf),NSR, xptxind1, ZU_R, DSDX) + CALL AKIMA(XUASECR, + & ZUASECR,NSR, xptxind1, ZU_R, DSDX) ! Compute the left aft node of panel @@ -1417,6 +1454,7 @@ subroutine update_surfaces() c-------------------------------------------------------------- use avl_heap_inc include 'AVL.INC' + INCLUDE 'AVL_surf.INC' integer ii NSTRIP = 0 @@ -1476,6 +1514,7 @@ SUBROUTINE MAKEBODY(IBODY) C using info from configuration input file. C-------------------------------------------------------------- INCLUDE 'AVL.INC' + INCLUDE 'AVL_surf.INC' C C REAL XBOD(IBX), YBOD(IBX), TBOD(IBX) C @@ -1574,6 +1613,7 @@ subroutine update_bodies() c Updates all bodies, using the stored data. c-------------------------------------------------------------- include 'AVL.INC' + INCLUDE 'AVL_surf.INC' integer IBODY, NBOD, NBLDS character*120 upname @@ -1615,6 +1655,7 @@ subroutine set_section_coordinates(isec,isurf,x,y,n,nin,xfmin, c Sets the airfoil coodinate data for the given section and surface c-------------------------------------------------------------- include 'AVL.INC' + INCLUDE 'AVL_surf.INC' c input integer isec, isurf, n, nin real x(n), y(n) @@ -1658,10 +1699,10 @@ subroutine set_section_coordinates(isec,isurf,x,y,n,nin,xfmin, & SASEC(i,isec,isurf)) call AKIMA(xin,tin,nin,XASEC(i,isec,isurf), & TASEC(i,isec,isurf),dummy) - XLASEC(i,isec,isurf) = XASEC(i,isec,isurf) - XUASEC(i,isec,isurf) = XASEC(i,isec,isurf) - ZLASEC(i,isec,isurf) = zc - 0.5*TASEC(i,isec,isurf) - ZUASEC(i,isec,isurf) = zc + 0.5*TASEC(i,isec,isurf) + ! XLASEC(i,isec,isurf) = XASEC(i,isec,isurf) + ! XUASEC(i,isec,isurf) = XASEC(i,isec,isurf) + ! ZLASEC(i,isec,isurf) = zc - 0.5*TASEC(i,isec,isurf) + ! ZUASEC(i,isec,isurf) = zc + 0.5*TASEC(i,isec,isurf) CASEC(i,isec,isurf) = zc end do @@ -1674,6 +1715,7 @@ subroutine set_body_coordinates(ibod,xb,yb,nb,nin,storecoords) c Sets the body oml coodinate data for the given section and surface c-------------------------------------------------------------- include 'AVL.INC' + INCLUDE 'AVL_surf.INC' integer ibod, nb, nin real xb(nb), yb(nb) logical storecoords @@ -1702,8 +1744,12 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) C reflected about y=Ypt. C----------------------------------- INCLUDE 'AVL.INC' + INCLUDE 'AVL_surf.INC' + INCLUDE 'AVL_oml.INC' CHARACTER*(*) MSG integer idx_vor + REAL :: Ypt + integer :: NN C C NNI = NN + 1 @@ -1923,6 +1969,7 @@ SUBROUTINE BDUPL(NN,Ypt,MSG) C reflected about y=Ypt. C----------------------------------- INCLUDE 'AVL.INC' + INCLUDE 'AVL_surf.INC' CHARACTER*(*) MSG C NNI = NBODY + 1 @@ -2009,6 +2056,7 @@ SUBROUTINE ENCALC C...COMMENTS C INCLUDE 'AVL.INC' + INCLUDE 'AVL_surf.INC' C REAL EP(3), EQ(3), ES(3), EB(3), EC(3), ECXB(3) REAL EC_G(3,NDMAX), ECXB_G(3) diff --git a/src/aoml.f b/src/aoml.f index 6bfbbd26..8d8dce33 100644 --- a/src/aoml.f +++ b/src/aoml.f @@ -23,6 +23,7 @@ SUBROUTINE CPOML(save_file) C INCLUDE 'AVL.INC' + INCLUDE 'AVL_oml.INC' C logical :: save_file LOGICAL LRANGEALL @@ -59,6 +60,7 @@ SUBROUTINE CPTHK C...COMMENTS C INCLUDE 'AVL.INC' + INCLUDE 'AVL_oml.INC' C PARAMETER (NCMAX=256) DIMENSION AICT(NCMAX,NCMAX) ! AIC for no normal flow BC @@ -186,6 +188,7 @@ SUBROUTINE CPDUMP(save_file) C...COMMENTS C code reader in read_cpoml.c C INCLUDE 'AVL.INC' + INCLUDE 'AVL_oml.INC' logical :: save_file integer idx_mesh, idx_mesh_surf, idx_LE, idx_LE_lo, idx_LE_up diff --git a/src/asetup.f b/src/asetup.f index 90df01ca..8ce5400c 100644 --- a/src/asetup.f +++ b/src/asetup.f @@ -520,34 +520,34 @@ SUBROUTINE VELSUM() C DO I = 1, NVOR DO K = 1, 3 - VC(K,I) = 0.0 + ! VC(K,I) = 0.0 VV(K,I) = 0.0 C--- h.v. velocity at control points and vortex midpoints DO J = 1, NVOR - VC(K,I) = VC(K,I) + WC_GAM(K,I,J)*GAM(J) + ! VC(K,I) = VC(K,I) + WC_GAM(K,I,J)*GAM(J) VV(K,I) = VV(K,I) + WV_GAM(K,I,J)*GAM(J) ENDDO DO N = 1, NUMAX - VC_U(K,I,N) = 0.0 + ! VC_U(K,I,N) = 0.0 VV_U(K,I,N) = 0.0 DO J = 1, NVOR - VC_U(K,I,N) = VC_U(K,I,N) + WC_GAM(K,I,J)*GAM_U(J,N) + ! VC_U(K,I,N) = VC_U(K,I,N) + WC_GAM(K,I,J)*GAM_U(J,N) VV_U(K,I,N) = VV_U(K,I,N) + WV_GAM(K,I,J)*GAM_U(J,N) ENDDO ENDDO DO N = 1, NCONTROL - VC_D(K,I,N) = 0.0 + ! VC_D(K,I,N) = 0.0 VV_D(K,I,N) = 0.0 DO J = 1, NVOR - VC_D(K,I,N) = VC_D(K,I,N) + WC_GAM(K,I,J)*GAM_D(J,N) + ! VC_D(K,I,N) = VC_D(K,I,N) + WC_GAM(K,I,J)*GAM_D(J,N) VV_D(K,I,N) = VV_D(K,I,N) + WV_GAM(K,I,J)*GAM_D(J,N) ENDDO ENDDO DO N = 1, NDESIGN - VC_G(K,I,N) = 0.0 + ! VC_G(K,I,N) = 0.0 VV_G(K,I,N) = 0.0 DO J = 1, NVOR - VC_G(K,I,N) = VC_G(K,I,N) + WC_GAM(K,I,J)*GAM_G(J,N) + ! VC_G(K,I,N) = VC_G(K,I,N) + WC_GAM(K,I,J)*GAM_G(J,N) VV_G(K,I,N) = VV_G(K,I,N) + WV_GAM(K,I,J)*GAM_G(J,N) ENDDO ENDDO @@ -565,18 +565,18 @@ SUBROUTINE VELSUM() & + WVSRD_U(K,I,5)*WROT(2) & + WVSRD_U(K,I,6)*WROT(3) C--- total velocity at control points and vortex midpoints - WC(K,I) = VC(K,I) + WCSRD(K,I) + ! WC(K,I) = VC(K,I) + WCSRD(K,I) WV(K,I) = VV(K,I) + WVSRD(K,I) DO N = 1, NUMAX - WC_U(K,I,N) = VC_U(K,I,N) + WCSRD_U(K,I,N) + ! WC_U(K,I,N) = VC_U(K,I,N) + WCSRD_U(K,I,N) WV_U(K,I,N) = VV_U(K,I,N) + WVSRD_U(K,I,N) ENDDO DO N = 1, NDMAX - WC_D(K,I,N) = VC_D(K,I,N) + ! WC_D(K,I,N) = VC_D(K,I,N) WV_D(K,I,N) = VV_D(K,I,N) ENDDO DO N = 1, NGMAX - WC_G(K,I,N) = VC_G(K,I,N) + ! WC_G(K,I,N) = VC_G(K,I,N) WV_G(K,I,N) = VV_G(K,I,N) ENDDO C diff --git a/src/avl.f b/src/avl.f index 90d55f01..36ae52e1 100644 --- a/src/avl.f +++ b/src/avl.f @@ -794,11 +794,11 @@ SUBROUTINE VARINI WROT(2) = 0. WROT(3) = 0. C - DO N = 1, NCONTROL + DO N = 1, NDMAX DELCON(N) = 0.0 ENDDO C - DO N = 1, NDESIGN + DO N = 1, NGMAX DELDES(N) = 0.0 ENDDO LSOL = .FALSE. @@ -1336,3 +1336,29 @@ SUBROUTINE AOCFIL(FNAME,IFILE) RETURN END + subroutine get_avl_constants( NVMAX_out, NSMAX_out, NSECMAX_out, + & NFMAX_out, NLMAX_out, NBMAX_out, NUMAX_out, + & NDMAX_out, NGMAX_out, NRMAX_out, NTMAX_out, + & NOBMAX_out, ICONX_out,IBX_out, NASMAX_out) + INCLUDE 'AVL.INC' + integer :: NVMAX_out, NSMAX_out, NSECMAX_out + integer :: NFMAX_out, NLMAX_out, NBMAX_out, NUMAX_out + integer :: NDMAX_out, NGMAX_out, NRMAX_out, NTMAX_out + integer :: NOBMAX_out, ICONX_out, IBX_out, NASMAX_out + + NVMAX_out = NVMAX + NSMAX_out = NSMAX + NSECMAX_out = NSECMAX + NFMAX_out = NFMAX + NLMAX_out = NLMAX + NBMAX_out = NBMAX + NUMAX_out = NUMAX + NDMAX_out = NDMAX + NGMAX_out = NGMAX + NRMAX_out = NRMAX + NTMAX_out = NTMAX + NOBMAX_out = NOBMAX + ICONX_out = ICONX + IBX_out = IBX + NASMAX_out = NASMAX + end subroutine get_avl_constants \ No newline at end of file diff --git a/src/avl_heap.f90 b/src/avl_heap.f90 index 18c8e49f..93944b68 100644 --- a/src/avl_heap.f90 +++ b/src/avl_heap.f90 @@ -8,14 +8,17 @@ subroutine avlheap_init(n) integer :: n ! Allocate AIC variable storage +! Use heap_allocated flag instead of allocated() to avoid false positives +! when Fortran allocatable descriptors contain garbage with -fno-init-global-zero - if (.not. allocated(AICN)) then + if (.not. heap_allocated) then allocate(AICN(n,n)) allocate(AICN_LU(n,n)) allocate(WC_GAM(3,n,n)) allocate(WV_GAM(3,n,n)) + heap_allocated = .TRUE. endif - + NAIC = n end subroutine avlheap_init @@ -28,13 +31,14 @@ subroutine avlheap_clean() ! Deallocate heap storage for AIC's - if (allocated(AICN)) then + if (heap_allocated) then deallocate(AICN) deallocate(AICN_LU) deallocate(WC_GAM) deallocate(WV_GAM) + heap_allocated = .FALSE. endif - + NAIC = -1 end subroutine avlheap_clean diff --git a/src/avl_heap_diff.f90 b/src/avl_heap_diff.f90 index 8d44378c..afd7622a 100644 --- a/src/avl_heap_diff.f90 +++ b/src/avl_heap_diff.f90 @@ -9,11 +9,15 @@ subroutine avlheap_diff_init(n) ! Allocate AIC variable storage - if (.not. allocated(AICN_DIFF)) then +! Use heap_diff_allocated flag instead of allocated() to avoid false positives +! when Fortran allocatable descriptors contain garbage with -fno-init-global-zero + + if (.not. heap_diff_allocated) then allocate(AICN_DIFF(n,n)) allocate(AICN_LU_DIFF(n,n)) allocate(WC_GAM_DIFF(3,n,n)) allocate(WV_GAM_DIFF(3,n,n)) + heap_diff_allocated = .TRUE. endif end subroutine avlheap_diff_init @@ -26,11 +30,12 @@ subroutine avlheap_diff_clean() ! Deallocate heap storage for AIC's - if (allocated(AICN_DIFF)) then + if (heap_diff_allocated) then deallocate(AICN_DIFF) deallocate(AICN_LU_DIFF) deallocate(WC_GAM_DIFF) deallocate(WV_GAM_DIFF) + heap_diff_allocated = .FALSE. endif end subroutine avlheap_diff_clean diff --git a/src/avl_heap_diff_inc.f90 b/src/avl_heap_diff_inc.f90 index 51052ab1..df141756 100644 --- a/src/avl_heap_diff_inc.f90 +++ b/src/avl_heap_diff_inc.f90 @@ -20,6 +20,9 @@ module avl_heap_diff_inc REAL(8), DIMENSION(:,:,:), ALLOCATABLE :: WV_GAM_DIFF !! !$omp threadprivate(WV_GAM) + + ! Explicit init ensures .data section placement (reliable with -fno-init-global-zero) + LOGICAL :: heap_diff_allocated = .FALSE. end module avl_heap_diff_inc diff --git a/src/avl_heap_inc.f90 b/src/avl_heap_inc.f90 index e82c4825..aa341c17 100644 --- a/src/avl_heap_inc.f90 +++ b/src/avl_heap_inc.f90 @@ -6,8 +6,9 @@ module avl_heap_inc ! Heap array include file for large AVL AIC arrays INTEGER, PARAMETER :: NVX=5000 ! must match nvmax in ADIMEN.INC) - INTEGER :: NAIC - + INTEGER :: NAIC = -1 + ! Explicit init ensures .data section placement (reliable with -fno-init-global-zero) + LOGICAL :: heap_allocated = .FALSE. ! All non-constant variables are declared as threadprivate for OpenMP ! warning! hardcoding precision! @@ -22,6 +23,6 @@ module avl_heap_inc REAL(8), DIMENSION(:,:,:), ALLOCATABLE :: WV_GAM !! !$omp threadprivate(WV_GAM) - + end module avl_heap_inc diff --git a/src/avl_heap_inc_diff.f90 b/src/avl_heap_inc_diff.f90 deleted file mode 100644 index 8307738e..00000000 --- a/src/avl_heap_inc_diff.f90 +++ /dev/null @@ -1,25 +0,0 @@ -! Declarations for heap memory storage for AVL -! Harold Youngren 6/21/2023 - -module avl_heap_inc_diff - -! Heap array include file for large AVL AIC arrays - - INTEGER, PARAMETER :: NVX=5000 ! must match nvmax in ADIMEN.INC) - -! All non-constant variables are declared as threadprivate for OpenMP - ! warning! hardcoding precision! - REAL(8), DIMENSION(:,:), ALLOCATABLE :: AICN_DIFF -!! !$omp threadprivate(AICN) - - REAL(8), DIMENSION(:,:), ALLOCATABLE :: AICN_LU_DIFF -!! !$omp threadprivate(AICN_LU) - - REAL(8), DIMENSION(:,:,:), ALLOCATABLE :: WC_GAM_DIFF -!! !$omp threadprivate(WC_GAM) - - REAL(8), DIMENSION(:,:,:), ALLOCATABLE :: WV_GAM_DIFF -!! !$omp threadprivate(WV_GAM) - - -end module avl_heap_inc_diff diff --git a/src/f2py/libavl.pyf b/src/f2py/libavl.pyf index c5c91579..ce5919af 100644 --- a/src/f2py/libavl.pyf +++ b/src/f2py/libavl.pyf @@ -4,10 +4,15 @@ python module libavl ! in interface ! in :libavl ! block data ! end block data - ! ! TODO: import from AVL.INC directly to avoid duplication subroutine avl ! in :libavl:avl.f + ! makesure you get that capitalization right! + ! certain platforms/compilers will be case sensitive + ! linux/gfortran for example include '../includes/AVL.INC' + include '../includes/AVL_surf.INC' + include '../includes/AVL_surf_ad_seeds.inc' + include '../includes/AVL_oml.INC' include '../includes/AVL_ad_seeds.inc' ! parameters from @@ -68,13 +73,13 @@ python module libavl ! in subroutine sdupl(nn,ypt,msg) ! in :libavl:amake.f integer :: nn - real :: ypt + real*8 :: ypt character*(*) msg end subroutine sdupl subroutine bdupl(nn,ypt,msg) ! in :libavl:amake.f integer :: nn - real :: ypt + real*8 :: ypt character*(*) msg end subroutine bdupl @@ -152,6 +157,12 @@ python module libavl ! in real*8 :: xb(nb), yb(nb) logical :: storecoords end subroutine set_body_coordinates - + + subroutine get_avl_constants(nvmax_out, nsmax_out, nsecmax_out, nfmax_out, nlmax_out, nbmax_out, numax_out, ndmax_out, ngmax_out, nrmax_out, ntmax_out, nobmax_out, iconx_out,ibx_out, nasmax_out) + integer, intent(out) :: nvmax_out, nsmax_out, nsecmax_out + integer, intent(out) :: nfmax_out, nlmax_out, nbmax_out, numax_out + integer, intent(out) :: ndmax_out, ngmax_out, nrmax_out, ntmax_out + integer, intent(out) :: nobmax_out, iconx_out,ibx_out, nasmax_out + end subroutine get_avl_constants end interface end python module libavl diff --git a/src/includes/ADIMEN.INC b/src/includes/ADIMEN.INC index cf2a3437..09aa1653 100644 --- a/src/includes/ADIMEN.INC +++ b/src/includes/ADIMEN.INC @@ -32,11 +32,8 @@ ! to make it easier to termine the true size in clear_ad_seeds_fast() ! I increased some values to avoid matching - ! PARAMETER (NVMAX=4000) PARAMETER (NVMAX=5000) - ! PARAMETER (NSMAX=400) ! max number of strips - ! PARAMETER (NSECMAX=401) ! max number of geometry sections PARAMETER (NSMAX=500) ! max number of strips PARAMETER (NSECMAX=301) ! max number of geometry sections PARAMETER (NFMAX=100) ! max number of surfaces @@ -54,6 +51,6 @@ ! PARAMETER (NOBMAX=1000) ! max number of off body points PARAMETER (ICONX=20) ! - PARAMETER (IBX=200) ! max number of airfoil coordinates - ! PARAMETER (IBX=300) ! max number of airfoil coordinates - + PARAMETER (IBX=300) ! max number of airfoil coordinates in a file + + PARAMETER (NASMAX=50) ! max number of airfoil data point diff --git a/src/includes/ADIMEN.INC.winArm64 b/src/includes/ADIMEN.INC.winArm64 new file mode 100644 index 00000000..ca091bca --- /dev/null +++ b/src/includes/ADIMEN.INC.winArm64 @@ -0,0 +1,56 @@ +! Parameters for Athena 3-D vortex lattice code +! +! Primary array limits +! NVMAX number of horseshoe vortices +! +! NSMAX number of chord strips +! NFMAX number of surfaces +! +! NLMAX number of source/doublet line nodes +! NBMAX number of bodies +! +! NUMAX number of freestream parameters (V,Omega) +! NDMAX number of control deflection parameters +! NGMAX number of design variables +! +! NRMAX number of stored run cases +! NTMAX number of stored time levels +! +! NOBMAX number of off-body points +! +! ICONX max number of control or design variable declaration lines per section (20) +! IBX max number of airfoil section points + INTEGER NVMAX + INTEGER NSMAX,NFMAX,NLMAX,NBMAX + INTEGER NUMAX,NDMAX,NGMAX + INTEGER NRMAX,NTMAX + INTEGER NOBMAX + INTEGER ICONX + INTEGER IBX + +! WARNING: This values have been tweaked compared to AVL +! to make it easier to termine the true size in clear_ad_seeds_fast() +! I increased some values to avoid matching + + PARAMETER (NVMAX=5000) + + PARAMETER (NSMAX=500) ! max number of strips + PARAMETER (NSECMAX=201) ! max number of geometry sections + PARAMETER (NFMAX=75) ! max number of surfaces + PARAMETER (NLMAX=502) ! max number of body nodes + PARAMETER (NBMAX=20) ! max number of bodies + + PARAMETER (NUMAX=6) ! number of velocity components + PARAMETER (NDMAX=15) ! max number of control surfaces + PARAMETER (NGMAX=11) ! max number of design variables + + PARAMETER (NRMAX=25) ! max number of run cases + PARAMETER (NTMAX=503) ! max number of time values + + PARAMETER (NOBMAX=1) ! max number of off body points + ! PARAMETER (NOBMAX=1) ! max number of off body points + + PARAMETER (ICONX=11) ! max number of control or design variable declaration lines per section + PARAMETER (IBX=300) ! max number of airfoil coordinates + + PARAMETER (NASMAX=50) ! max number of airfoil data point \ No newline at end of file diff --git a/src/includes/AVL.INC b/src/includes/AVL.INC index fc47c9cc..e042fdad 100644 --- a/src/includes/AVL.INC +++ b/src/includes/AVL.INC @@ -44,10 +44,10 @@ C & FEVDEF, ! default eigenvalue save file & TITLE, ! configuration title & STITLE(NFMAX), ! surface title - & AFILES(NSMAX,NFMAX), ! airfoil file names + & AFILES(NSECMAX,NFMAX), ! airfoil file names & BTITLE(NBMAX), ! body title & BFILES(NBMAX), ! body file names - & NACA(NSMAX,NFMAX), ! naca section names + & NACA(NSECMAX,NFMAX), ! naca section names & RTITLE(NRMAX), ! run case title & DNAME(NDMAX), ! control variable name & GNAME(NGMAX), ! design variable name @@ -395,144 +395,6 @@ C & CF_LSRF(3,NFMAX), CM_LSRF(3,NFMAX), & CL_LSRF(NFMAX), CD_LSRF(NFMAX), CMLE_LSRF(NFMAX) -C !start added variables for python geometry minipulation - - LOGICAL LDUPL_B - COMMON /BODY_GEOM_L/ - & LDUPL_B(NBMAX) ! T if body is a duplicated -c - COMMON /BODY_GEOM_I/ - & NSEC_B(NBMAX), ! number of sections in body - & NVB(NFMAX), ! number of body elements - & NBOD(NBMAX) ! number of body oml points - - REAL(kind=avl_real) XYZSCAL_B - REAL(kind=avl_real) XYZTRAN_B - REAL(kind=avl_real) YDUPL_B - REAL(kind=avl_real) XYZLES_B - REAL(kind=avl_real) BSPACE - REAL(kind=avl_real) XBOD - REAL(kind=avl_real) YBOD - REAL(kind=avl_real) TBOD - REAL(kind=avl_real) XBOD_R - REAL(kind=avl_real) YBOD_R - COMMON /BODY_GEOM_R/ - & XYZSCAL_B(3, NBMAX), ! scaling factors for XYZ coordinates - & XYZTRAN_B(3, NBMAX), ! translation factors for XYZ coordinates - & YDUPL_B(NBMAX), ! y duplicate - & XYZLES_B(3, NSMAX, NBMAX), ! leading edge cordinate vector - & BSPACE(NBMAX), ! body spacing - & XBOD(IBX, NBMAX), ! body oml x-coordinates - & YBOD(IBX, NBMAX), ! body oml y-coordinates - & TBOD(IBX, NBMAX), ! body oml thickness - & XBOD_R(IBX, NBMAX), ! raw input oml x-coordinates - & YBOD_R(IBX, NBMAX) ! raw input oml y-coordinates -c - LOGICAL LDUPL - LOGICAL LSURFSPACING - COMMON /SURF_GEOM_L/ - & LDUPL(NFMAX), ! T if surface is a duplicated - & LSURFSPACING(NFMAX) ! surface spacing set under the surface heeading - - COMMON /SURF_GEOM_I/ - & NSEC(NFMAX), ! number of sections in surface - & NVC(NFMAX), ! number of chordwise elements - & NVS(NFMAX), ! number of spanwise elements - & NSPANS(NSECMAX, NFMAX), ! number of spanwise elements vector - & NASEC(NSECMAX, NFMAX), ! number of points used to represent section geometry - & NAPTSSEC(NSECMAX, IBX), ! number of points of input airfoil geometry - & ICONTD(ICONX, NSECMAX, NFMAX), ! control variable index - & NSCON(NSECMAX, NFMAX), ! number of control variables - & IDESTD(ICONX, NSECMAX, NFMAX), ! design variable index - & NSDES(NSECMAX, NFMAX) ! number of design variables - - REAL(kind=avl_real) XYZSCAL - REAL(kind=avl_real) XYZTRAN - REAL(kind=avl_real) YDUPL - REAL(kind=avl_real) ADDINC - REAL(kind=avl_real) CSPACE - REAL(kind=avl_real) SSPACE - REAL(kind=avl_real) SSPACES - REAL(kind=avl_real) XYZLES - REAL(kind=avl_real) XSEC - REAL(kind=avl_real) YSEC - REAL(kind=avl_real) XFMIN_R - REAL(kind=avl_real) XFMAX_R - REAL(kind=avl_real) XLASEC - REAL(kind=avl_real) ZLASEC - REAL(kind=avl_real) XUASEC - REAL(kind=avl_real) ZUASEC - REAL(kind=avl_real) CHORDS - REAL(kind=avl_real) AINCS - REAL(kind=avl_real) XASEC - REAL(kind=avl_real) CASEC - REAL(kind=avl_real) SASEC - REAL(kind=avl_real) TASEC - REAL(kind=avl_real) CLCDSEC - REAL(kind=avl_real) CLCDSRF - REAL(kind=avl_real) CLAF - REAL(kind=avl_real) XHINGED - REAL(kind=avl_real) VHINGED - REAL(kind=avl_real) GAIND - REAL(kind=avl_real) REFLD - REAL(kind=avl_real) GAING - COMMON /SURF_GEOM_R/ - & XYZSCAL(3, NFMAX), ! scaling factors for XYZ coordinates - & XYZTRAN(3, NFMAX), ! translation factors for XYZ coordinates - & YDUPL(NFMAX), ! y position of duplicate plane - & ADDINC(NFMAX), ! additional incidence angle - & CSPACE(NFMAX), ! chordwise spacing - & SSPACE(NFMAX), ! spanwise spacing - & SSPACES(NSECMAX, NFMAX), ! spanwise spacing vector - & XYZLES(3, NSECMAX, NFMAX), ! leading edge cordinate vector - & XSEC(IBX, NSECMAX, NFMAX), ! raw section input x coordinates - & YSEC(IBX, NSECMAX, NFMAX), ! raw section input y coordinates - & XFMIN_R(NSECMAX, NFMAX), ! section input lower x bound - & XFMAX_R(NSECMAX, NFMAX), ! section input upper x bound - & XLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface x-coords - & ZLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface z-coords - & XUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface x-coords - & ZUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface z-coords - & CHORDS(NSECMAX, NFMAX), ! chord length vectorm - & AINCS(NSECMAX, NFMAX), ! incidence angle vector - & XASEC(IBX, NSECMAX, NFMAX), ! the x coordinate aifoil section - & CASEC(IBX, NSECMAX, NFMAX), ! camber line at xasec - & SASEC(IBX, NSECMAX, NFMAX), ! slope of camber line at xasec - & TASEC(IBX, NSECMAX, NFMAX), ! thickness at xasec - & CLCDSEC(6, NSECMAX, NFMAX), - & CLCDSRF(6, NFMAX), - & CLAF(NSECMAX, NFMAX), ! CL alpha (dCL/da) scaling factor - & XHINGED(ICONX, NSECMAX, NFMAX), ! hinge location - & VHINGED(3, ICONX, NSECMAX, NFMAX), ! hinge vector - & GAIND(ICONX, NSECMAX, NFMAX), ! control surface gain - & REFLD(ICONX, NSECMAX, NFMAX), ! control surface reflection - & GAING(ICONX, NSECMAX, NFMAX) ! desgin variable gain - - COMMON /SURF_MESH_I/ - & MFRST(NFMAX), ! stores the index in the MSHBLK where each surface's mesh begins - & IPTSEC(NSMAX,NFMAX) ! stores the iptloc vector for each surface - - REAL(kind=avl_real) MSHBLK - REAL(kind=avl_real) RV1MSH - REAL(kind=avl_real) RV2MSH - REAL(kind=avl_real) RVMSH - REAL(kind=avl_real) RCMSH - COMMON /SURF_MESH_R/ - & MSHBLK(3, 4*NVMAX), ! block to store all surface meshes - & RV1MSH(3,NVMAX), ! mesh h.v. vortex left points - & RV2MSH(3,NVMAX), ! mesh h.v. vortex right points - & RVMSH(3,NVMAX), ! mesh h.v. vortex center points - & RCMSH(3,NVMAX) ! mesh h.v. control points - - LOGICAL LSURFMSH - LOGICAL LMESHFLAT - COMMON /SURF_MESH_L/ - & LSURFMSH(NFMAX), ! T if surface uses a mesh cordinates for its geometry - & LMESHFLAT(NFMAX) ! T if the surface's mesh should be flattened when computing vorticies and control points - - -C !!--- end added variables for python geometry manipulation --- - COMMON /STRP_I/ & LSSURF(NSMAX), ! index of surface which contains this strip & IJFRST(NSMAX), ! index of first element in strip @@ -738,22 +600,6 @@ cc#else & RES_U(NVMAX,NUMAX), ! residual for gamma_u lin sys & RES_D(NVMAX,NDMAX) ! residual for gamma_d lin. sys -cc#ifdef USE_CPOML - REAL(kind=avl_real) XYN1 - REAL(kind=avl_real) XYN2 - REAL(kind=avl_real) ZLON1 - REAL(kind=avl_real) ZLON2 - REAL(kind=avl_real) ZUPN1 - REAL(kind=avl_real) ZUPN2 - REAL(kind=avl_real) XYZSURF - REAL(kind=avl_real) CPSURF - COMMON /VRTX_S/ - & XYN1(2,NVMAX), XYN2(2,NVMAX), ! left/right aft-node of element - & ZLON1(NVMAX), ZLON2(NVMAX), ! left/right lower z-coord of aft-node - & ZUPN1(NVMAX), ZUPN2(NVMAX), ! left/right upper z-coord of aft-node - & XYZSURF(3, NVMAX, NFMAX), - & CPSURF(NVMAX,NFMAX) -cc#endif C COMMON /BODY_I/ & NL(NBMAX), ! number of source-line nodes in body @@ -782,13 +628,13 @@ C REAL(kind=avl_real) AMACH REAL(kind=avl_real) VC - REAL(kind=avl_real) VC_U - REAL(kind=avl_real) VC_D - REAL(kind=avl_real) VC_G - REAL(kind=avl_real) WC - REAL(kind=avl_real) WC_U - REAL(kind=avl_real) WC_D - REAL(kind=avl_real) WC_G +c REAL(kind=avl_real) VC_U +c REAL(kind=avl_real) VC_D +c REAL(kind=avl_real) VC_G +c REAL(kind=avl_real) WC +c REAL(kind=avl_real) WC_U +c REAL(kind=avl_real) WC_D +c REAL(kind=avl_real) WC_G REAL(kind=avl_real) WV REAL(kind=avl_real) WV_U REAL(kind=avl_real) WV_D @@ -805,14 +651,14 @@ cc & AICN(NVMAX,NVMAX), ! normalwash AIC matrix (and VL system matrix) cc & WC_GAM(3,NVMAX,NVMAX), ! c.p. velocity/Gamma influence matrix cc & WV_GAM(3,NVMAX,NVMAX), ! h.v. velocity/Gamma influence matrix & VC(3,NVMAX), ! h.v. induced velocity at c.p. - & VC_U(3,NVMAX,NUMAX), - & VC_D(3,NVMAX,NDMAX), - & VC_G(3,NVMAX,NGMAX), +c & VC_U(3,NVMAX,NUMAX), +c & VC_D(3,NVMAX,NDMAX), +c & VC_G(3,NVMAX,NGMAX), C - & WC(3,NVMAX), ! total induced velocity at c.p. (used in residual) - & WC_U(3,NVMAX,NUMAX), - & WC_D(3,NVMAX,NDMAX), - & WC_G(3,NVMAX,NGMAX), +c & WC(3,NVMAX), ! total induced velocity at c.p. (used in residual) +C & WC_U(3,NVMAX,NUMAX), +C & WC_D(3,NVMAX,NDMAX), +C & WC_G(3,NVMAX,NGMAX), C & VV(3,NVMAX), ! h.v. induced velocity at h.v. & VV_U(3,NVMAX,NUMAX), diff --git a/src/includes/AVL_ad_seeds.inc b/src/includes/AVL_ad_seeds.inc index dae2b9fb..4cbe5529 100644 --- a/src/includes/AVL_ad_seeds.inc +++ b/src/includes/AVL_ad_seeds.inc @@ -1,5 +1,5 @@ C automatically generated by gen_ad_inc.py -C must be precended by 'include AVL_ad_seed' +C must be precended by 'AVL.INC' real(kind=avl_real) YSYM_DIFF real(kind=avl_real) ZSYM_DIFF real(kind=avl_real) ALFA_DIFF @@ -11,36 +11,16 @@ C must be precended by 'include AVL_ad_seed' real(kind=avl_real) PARVAL_DIFF real(kind=avl_real) CONVAL_DIFF real(kind=avl_real) DELCON_DIFF - real(kind=avl_real) DELDES_DIFF real(kind=avl_real) SREF_DIFF real(kind=avl_real) CREF_DIFF real(kind=avl_real) BREF_DIFF real(kind=avl_real) XYZREF_DIFF - real(kind=avl_real) XYZREF0_DIFF real(kind=avl_real) MACH_DIFF - real(kind=avl_real) MACH0_DIFF real(kind=avl_real) CDREF_DIFF - real(kind=avl_real) CDREF0_DIFF - real(kind=avl_real) VRCOREC_DIFF - real(kind=avl_real) VRCOREW_DIFF - real(kind=avl_real) SRCORE_DIFF real(kind=avl_real) CLFF_DIFF real(kind=avl_real) CYFF_DIFF real(kind=avl_real) CDFF_DIFF - real(kind=avl_real) CLFF_U_DIFF - real(kind=avl_real) CYFF_U_DIFF - real(kind=avl_real) CDFF_U_DIFF - real(kind=avl_real) CLFF_D_DIFF - real(kind=avl_real) CYFF_D_DIFF - real(kind=avl_real) CDFF_D_DIFF - real(kind=avl_real) CLFF_G_DIFF - real(kind=avl_real) CYFF_G_DIFF - real(kind=avl_real) CDFF_G_DIFF real(kind=avl_real) SPANEF_DIFF - real(kind=avl_real) SPANEF_A_DIFF - real(kind=avl_real) SPANEF_U_DIFF - real(kind=avl_real) SPANEF_D_DIFF - real(kind=avl_real) SPANEF_G_DIFF real(kind=avl_real) CDTOT_DIFF real(kind=avl_real) CYTOT_DIFF real(kind=avl_real) CLTOT_DIFF @@ -52,13 +32,9 @@ C must be precended by 'include AVL_ad_seed' real(kind=avl_real) CDTOT_D_DIFF real(kind=avl_real) CYTOT_D_DIFF real(kind=avl_real) CLTOT_D_DIFF - real(kind=avl_real) CDTOT_G_DIFF - real(kind=avl_real) CYTOT_G_DIFF - real(kind=avl_real) CLTOT_G_DIFF real(kind=avl_real) CFTOT_DIFF real(kind=avl_real) CFTOT_U_DIFF real(kind=avl_real) CFTOT_D_DIFF - real(kind=avl_real) CFTOT_G_DIFF real(kind=avl_real) CXBAX_DIFF real(kind=avl_real) CYBAX_DIFF real(kind=avl_real) CZBAX_DIFF @@ -68,15 +44,10 @@ C must be precended by 'include AVL_ad_seed' real(kind=avl_real) CXTOT_D_BA_DIFF real(kind=avl_real) CYTOT_D_BA_DIFF real(kind=avl_real) CZTOT_D_BA_DIFF - real(kind=avl_real) CXTOT_G_BA_DIFF - real(kind=avl_real) CYTOT_G_BA_DIFF - real(kind=avl_real) CZTOT_G_BA_DIFF real(kind=avl_real) CMTOT_DIFF real(kind=avl_real) CMTOT_U_DIFF real(kind=avl_real) CMTOT_D_DIFF - real(kind=avl_real) CMTOT_G_DIFF real(kind=avl_real) CRSAX_DIFF - real(kind=avl_real) CMSAX_DIFF real(kind=avl_real) CNSAX_DIFF real(kind=avl_real) CRSAX_D_DIFF real(kind=avl_real) CMSAX_D_DIFF @@ -90,9 +61,6 @@ C must be precended by 'include AVL_ad_seed' real(kind=avl_real) CRTOT_D_BA_DIFF real(kind=avl_real) CMTOT_D_BA_DIFF real(kind=avl_real) CNTOT_D_BA_DIFF - real(kind=avl_real) CRTOT_G_BA_DIFF - real(kind=avl_real) CMTOT_G_BA_DIFF - real(kind=avl_real) CNTOT_G_BA_DIFF real(kind=avl_real) CDVTOT_DIFF real(kind=avl_real) CDITOT_DIFF real(kind=avl_real) CDTOT_AL_DIFF @@ -125,14 +93,6 @@ C must be precended by 'include AVL_ad_seed' real(kind=avl_real) CRTOT_RZ_DIFF real(kind=avl_real) CMTOT_RZ_DIFF real(kind=avl_real) CNTOT_RZ_DIFF - real(kind=avl_real) CHINGE_DIFF - real(kind=avl_real) CHINGE_U_DIFF - real(kind=avl_real) CHINGE_D_DIFF - real(kind=avl_real) CHINGE_G_DIFF - real(kind=avl_real) DCL_A0_DIFF - real(kind=avl_real) DCM_A0_DIFF - real(kind=avl_real) DCL_U0_DIFF - real(kind=avl_real) DCM_U0_DIFF real(kind=avl_real) XNP_DIFF real(kind=avl_real) SM_DIFF real(kind=avl_real) BB_DIFF @@ -149,36 +109,16 @@ C must be precended by 'include AVL_ad_seed' & PARVAL_DIFF(IPMAX,NRMAX), & CONVAL_DIFF(ICMAX,NRMAX), & DELCON_DIFF(NDMAX), - & DELDES_DIFF(NGMAX), & SREF_DIFF, & CREF_DIFF, & BREF_DIFF, & XYZREF_DIFF(3), - & XYZREF0_DIFF(3), & MACH_DIFF, - & MACH0_DIFF, & CDREF_DIFF, - & CDREF0_DIFF, - & VRCOREC_DIFF, - & VRCOREW_DIFF, - & SRCORE_DIFF, & CLFF_DIFF, & CYFF_DIFF, & CDFF_DIFF, - & CLFF_U_DIFF(NUMAX), - & CYFF_U_DIFF(NUMAX), - & CDFF_U_DIFF(NUMAX), - & CLFF_D_DIFF(NDMAX), - & CYFF_D_DIFF(NDMAX), - & CDFF_D_DIFF(NDMAX), - & CLFF_G_DIFF(NGMAX), - & CYFF_G_DIFF(NGMAX), - & CDFF_G_DIFF(NGMAX), & SPANEF_DIFF, - & SPANEF_A_DIFF, - & SPANEF_U_DIFF(NUMAX), - & SPANEF_D_DIFF(NDMAX), - & SPANEF_G_DIFF(NGMAX), & CDTOT_DIFF, & CYTOT_DIFF, & CLTOT_DIFF, @@ -190,13 +130,9 @@ C must be precended by 'include AVL_ad_seed' & CDTOT_D_DIFF(NDMAX), & CYTOT_D_DIFF(NDMAX), & CLTOT_D_DIFF(NDMAX), - & CDTOT_G_DIFF(NGMAX), - & CYTOT_G_DIFF(NGMAX), - & CLTOT_G_DIFF(NGMAX), & CFTOT_DIFF(3), & CFTOT_U_DIFF(3,NUMAX), & CFTOT_D_DIFF(3,NDMAX), - & CFTOT_G_DIFF(3,NGMAX), & CXBAX_DIFF, & CYBAX_DIFF, & CZBAX_DIFF, @@ -206,15 +142,10 @@ C must be precended by 'include AVL_ad_seed' & CXTOT_D_BA_DIFF(NDMAX), & CYTOT_D_BA_DIFF(NDMAX), & CZTOT_D_BA_DIFF(NDMAX), - & CXTOT_G_BA_DIFF(NGMAX), - & CYTOT_G_BA_DIFF(NGMAX), - & CZTOT_G_BA_DIFF(NGMAX), & CMTOT_DIFF(3), & CMTOT_U_DIFF(3,NUMAX), & CMTOT_D_DIFF(3,NDMAX), - & CMTOT_G_DIFF(3,NGMAX), & CRSAX_DIFF, - & CMSAX_DIFF, & CNSAX_DIFF, & CRSAX_D_DIFF(NDMAX), & CMSAX_D_DIFF(NDMAX), @@ -228,9 +159,6 @@ C must be precended by 'include AVL_ad_seed' & CRTOT_D_BA_DIFF(NDMAX), & CMTOT_D_BA_DIFF(NDMAX), & CNTOT_D_BA_DIFF(NDMAX), - & CRTOT_G_BA_DIFF(NGMAX), - & CMTOT_G_BA_DIFF(NGMAX), - & CNTOT_G_BA_DIFF(NGMAX), & CDVTOT_DIFF, & CDITOT_DIFF, & CDTOT_AL_DIFF, @@ -263,40 +191,15 @@ C must be precended by 'include AVL_ad_seed' & CRTOT_RZ_DIFF, & CMTOT_RZ_DIFF, & CNTOT_RZ_DIFF, - & CHINGE_DIFF(NDMAX), - & CHINGE_U_DIFF(NDMAX,NUMAX), - & CHINGE_D_DIFF(NDMAX,NDMAX), - & CHINGE_G_DIFF(NDMAX,NGMAX), - & DCL_A0_DIFF, - & DCM_A0_DIFF, - & DCL_U0_DIFF, - & DCM_U0_DIFF, & XNP_DIFF, & SM_DIFF, & BB_DIFF, & RR_DIFF -C - real(kind=avl_real) RHO0_DIFF - real(kind=avl_real) GEE0_DIFF - real(kind=avl_real) XYZMASS0_DIFF - real(kind=avl_real) RMASS0_DIFF - real(kind=avl_real) RINER0_DIFF - real(kind=avl_real) AMASS_DIFF - real(kind=avl_real) AINER_DIFF - COMMON /MASS_R_DIFF/ - & RHO0_DIFF, - & GEE0_DIFF, - & XYZMASS0_DIFF(3), - & RMASS0_DIFF, - & RINER0_DIFF(3,3), - & AMASS_DIFF(3,3), - & AINER_DIFF(3,3) C real(kind=avl_real) CDSURF_DIFF real(kind=avl_real) CYSURF_DIFF real(kind=avl_real) CLSURF_DIFF real(kind=avl_real) CDS_A_DIFF - real(kind=avl_real) CYS_A_DIFF real(kind=avl_real) CLS_A_DIFF real(kind=avl_real) CDS_U_DIFF real(kind=avl_real) CYS_U_DIFF @@ -304,33 +207,18 @@ C real(kind=avl_real) CDS_D_DIFF real(kind=avl_real) CYS_D_DIFF real(kind=avl_real) CLS_D_DIFF - real(kind=avl_real) CDS_G_DIFF - real(kind=avl_real) CYS_G_DIFF - real(kind=avl_real) CLS_G_DIFF real(kind=avl_real) CFSURF_DIFF real(kind=avl_real) CFS_U_DIFF real(kind=avl_real) CFS_D_DIFF - real(kind=avl_real) CFS_G_DIFF real(kind=avl_real) CMSURF_DIFF - real(kind=avl_real) CMSURFBAX_DIFF real(kind=avl_real) CMS_U_DIFF real(kind=avl_real) CMS_D_DIFF - real(kind=avl_real) CMS_G_DIFF real(kind=avl_real) CDVSURF_DIFF - real(kind=avl_real) CDISURF_DIFF - real(kind=avl_real) SSURF_DIFF - real(kind=avl_real) CAVESURF_DIFF - real(kind=avl_real) CF_LSRF_DIFF - real(kind=avl_real) CM_LSRF_DIFF - real(kind=avl_real) CL_LSRF_DIFF - real(kind=avl_real) CD_LSRF_DIFF - real(kind=avl_real) CMLE_LSRF_DIFF COMMON /SURF_R_DIFF/ & CDSURF_DIFF(NFMAX), & CYSURF_DIFF(NFMAX), & CLSURF_DIFF(NFMAX), & CDS_A_DIFF(NFMAX), - & CYS_A_DIFF(NFMAX), & CLS_A_DIFF(NFMAX), & CDS_U_DIFF(NFMAX,NUMAX), & CYS_U_DIFF(NFMAX,NUMAX), @@ -338,123 +226,13 @@ C & CDS_D_DIFF(NFMAX,NDMAX), & CYS_D_DIFF(NFMAX,NDMAX), & CLS_D_DIFF(NFMAX,NDMAX), - & CDS_G_DIFF(NFMAX,NGMAX), - & CYS_G_DIFF(NFMAX,NGMAX), - & CLS_G_DIFF(NFMAX,NGMAX), & CFSURF_DIFF(3,NFMAX), & CFS_U_DIFF(3,NFMAX,NUMAX), & CFS_D_DIFF(3,NFMAX,NDMAX), - & CFS_G_DIFF(3,NFMAX,NGMAX), & CMSURF_DIFF(3,NFMAX), - & CMSURFBAX_DIFF(3,NFMAX), & CMS_U_DIFF(3,NFMAX,NUMAX), & CMS_D_DIFF(3,NFMAX,NDMAX), - & CMS_G_DIFF(3,NFMAX,NGMAX), - & CDVSURF_DIFF(NFMAX), - & CDISURF_DIFF(NFMAX), - & SSURF_DIFF(NFMAX), - & CAVESURF_DIFF(NFMAX), - & CF_LSRF_DIFF(3,NFMAX), - & CM_LSRF_DIFF(3,NFMAX), - & CL_LSRF_DIFF(NFMAX), - & CD_LSRF_DIFF(NFMAX), - & CMLE_LSRF_DIFF(NFMAX) -C - real(kind=avl_real) XYZSCAL_B_DIFF - real(kind=avl_real) XYZTRAN_B_DIFF - real(kind=avl_real) YDUPL_B_DIFF - real(kind=avl_real) XYZLES_B_DIFF - real(kind=avl_real) BSPACE_DIFF - real(kind=avl_real) XBOD_DIFF - real(kind=avl_real) YBOD_DIFF - real(kind=avl_real) TBOD_DIFF - real(kind=avl_real) XBOD_R_DIFF - real(kind=avl_real) YBOD_R_DIFF - COMMON /BODY_GEOM_R_DIFF/ - & XYZSCAL_B_DIFF(3, NBMAX), - & XYZTRAN_B_DIFF(3, NBMAX), - & YDUPL_B_DIFF(NBMAX), - & XYZLES_B_DIFF(3, NSMAX, NBMAX), - & BSPACE_DIFF(NBMAX), - & XBOD_DIFF(IBX, NBMAX), - & YBOD_DIFF(IBX, NBMAX), - & TBOD_DIFF(IBX, NBMAX), - & XBOD_R_DIFF(IBX, NBMAX), - & YBOD_R_DIFF(IBX, NBMAX) -C - real(kind=avl_real) XYZSCAL_DIFF - real(kind=avl_real) XYZTRAN_DIFF - real(kind=avl_real) YDUPL_DIFF - real(kind=avl_real) ADDINC_DIFF - real(kind=avl_real) CSPACE_DIFF - real(kind=avl_real) SSPACE_DIFF - real(kind=avl_real) SSPACES_DIFF - real(kind=avl_real) XYZLES_DIFF - real(kind=avl_real) XSEC_DIFF - real(kind=avl_real) YSEC_DIFF - real(kind=avl_real) XFMIN_R_DIFF - real(kind=avl_real) XFMAX_R_DIFF - real(kind=avl_real) XLASEC_DIFF - real(kind=avl_real) ZLASEC_DIFF - real(kind=avl_real) XUASEC_DIFF - real(kind=avl_real) ZUASEC_DIFF - real(kind=avl_real) CHORDS_DIFF - real(kind=avl_real) AINCS_DIFF - real(kind=avl_real) XASEC_DIFF - real(kind=avl_real) CASEC_DIFF - real(kind=avl_real) SASEC_DIFF - real(kind=avl_real) TASEC_DIFF - real(kind=avl_real) CLCDSEC_DIFF - real(kind=avl_real) CLCDSRF_DIFF - real(kind=avl_real) CLAF_DIFF - real(kind=avl_real) XHINGED_DIFF - real(kind=avl_real) VHINGED_DIFF - real(kind=avl_real) GAIND_DIFF - real(kind=avl_real) REFLD_DIFF - real(kind=avl_real) GAING_DIFF - COMMON /SURF_GEOM_R_DIFF/ - & XYZSCAL_DIFF(3, NFMAX), - & XYZTRAN_DIFF(3, NFMAX), - & YDUPL_DIFF(NFMAX), - & ADDINC_DIFF(NFMAX), - & CSPACE_DIFF(NFMAX), - & SSPACE_DIFF(NFMAX), - & SSPACES_DIFF(NSECMAX, NFMAX), - & XYZLES_DIFF(3, NSECMAX, NFMAX), - & XSEC_DIFF(IBX, NSECMAX, NFMAX), - & YSEC_DIFF(IBX, NSECMAX, NFMAX), - & XFMIN_R_DIFF(NSECMAX, NFMAX), - & XFMAX_R_DIFF(NSECMAX, NFMAX), - & XLASEC_DIFF(IBX, NSECMAX, NFMAX), - & ZLASEC_DIFF(IBX, NSECMAX, NFMAX), - & XUASEC_DIFF(IBX, NSECMAX, NFMAX), - & ZUASEC_DIFF(IBX, NSECMAX, NFMAX), - & CHORDS_DIFF(NSECMAX, NFMAX), - & AINCS_DIFF(NSECMAX, NFMAX), - & XASEC_DIFF(IBX, NSECMAX, NFMAX), - & CASEC_DIFF(IBX, NSECMAX, NFMAX), - & SASEC_DIFF(IBX, NSECMAX, NFMAX), - & TASEC_DIFF(IBX, NSECMAX, NFMAX), - & CLCDSEC_DIFF(6, NSECMAX, NFMAX), - & CLCDSRF_DIFF(6, NFMAX), - & CLAF_DIFF(NSECMAX, NFMAX), - & XHINGED_DIFF(ICONX, NSECMAX, NFMAX), - & VHINGED_DIFF(3, ICONX, NSECMAX, NFMAX), - & GAIND_DIFF(ICONX, NSECMAX, NFMAX), - & REFLD_DIFF(ICONX, NSECMAX, NFMAX), - & GAING_DIFF(ICONX, NSECMAX, NFMAX) -C - real(kind=avl_real) MSHBLK_DIFF - real(kind=avl_real) RV1MSH_DIFF - real(kind=avl_real) RV2MSH_DIFF - real(kind=avl_real) RVMSH_DIFF - real(kind=avl_real) RCMSH_DIFF - COMMON /SURF_MESH_R_DIFF/ - & MSHBLK_DIFF(3, 4*NVMAX), - & RV1MSH_DIFF(3,NVMAX), - & RV2MSH_DIFF(3,NVMAX), - & RVMSH_DIFF(3,NVMAX), - & RCMSH_DIFF(3,NVMAX) + & CDVSURF_DIFF(NFMAX) C real(kind=avl_real) RLE_DIFF real(kind=avl_real) CHORD_DIFF @@ -463,11 +241,7 @@ C real(kind=avl_real) RLE2_DIFF real(kind=avl_real) CHORD2_DIFF real(kind=avl_real) WSTRIP_DIFF - real(kind=avl_real) TANLE_DIFF - real(kind=avl_real) TANTE_DIFF - real(kind=avl_real) GINCSTRIP_DIFF real(kind=avl_real) CLCD_DIFF - real(kind=avl_real) SAXFR_DIFF real(kind=avl_real) ESS_DIFF real(kind=avl_real) ENSY_DIFF real(kind=avl_real) ENSZ_DIFF @@ -480,7 +254,6 @@ C real(kind=avl_real) CYSTRP_DIFF real(kind=avl_real) CLSTRP_DIFF real(kind=avl_real) CDST_A_DIFF - real(kind=avl_real) CYST_A_DIFF real(kind=avl_real) CLST_A_DIFF real(kind=avl_real) CDST_U_DIFF real(kind=avl_real) CYST_U_DIFF @@ -488,33 +261,13 @@ C real(kind=avl_real) CDST_D_DIFF real(kind=avl_real) CYST_D_DIFF real(kind=avl_real) CLST_D_DIFF - real(kind=avl_real) CDST_G_DIFF - real(kind=avl_real) CYST_G_DIFF - real(kind=avl_real) CLST_G_DIFF real(kind=avl_real) CFSTRP_DIFF real(kind=avl_real) CFST_U_DIFF real(kind=avl_real) CFST_D_DIFF - real(kind=avl_real) CFST_G_DIFF real(kind=avl_real) CMSTRP_DIFF real(kind=avl_real) CMST_U_DIFF real(kind=avl_real) CMST_D_DIFF - real(kind=avl_real) CMST_G_DIFF - real(kind=avl_real) CF_LSTRP_DIFF - real(kind=avl_real) CM_LSTRP_DIFF - real(kind=avl_real) CN_LSTRP_DIFF - real(kind=avl_real) CA_LSTRP_DIFF - real(kind=avl_real) CD_LSTRP_DIFF - real(kind=avl_real) CL_LSTRP_DIFF real(kind=avl_real) CDV_LSTRP_DIFF - real(kind=avl_real) CLT_LSTRP_DIFF - real(kind=avl_real) CLA_LSTRP_DIFF - real(kind=avl_real) CMC4_LSTRP_DIFF - real(kind=avl_real) CMLE_LSTRP_DIFF - real(kind=avl_real) CNC_DIFF - real(kind=avl_real) DWWAKE_DIFF - real(kind=avl_real) CNC_U_DIFF - real(kind=avl_real) CNC_D_DIFF - real(kind=avl_real) CNC_G_DIFF COMMON /STRP_R_DIFF/ & RLE_DIFF(3,NSMAX), & CHORD_DIFF(NSMAX), @@ -523,11 +276,7 @@ C & RLE2_DIFF(3,NSMAX), & CHORD2_DIFF(NSMAX), & WSTRIP_DIFF(NSMAX), - & TANLE_DIFF(NSMAX), - & TANTE_DIFF(NSMAX), - & GINCSTRIP_DIFF(NSMAX), & CLCD_DIFF(6,NSMAX), - & SAXFR_DIFF, & ESS_DIFF(3,NSMAX), & ENSY_DIFF(NSMAX), & ENSZ_DIFF(NSMAX), @@ -540,7 +289,6 @@ C & CYSTRP_DIFF(NSMAX), & CLSTRP_DIFF(NSMAX), & CDST_A_DIFF(NSMAX), - & CYST_A_DIFF(NSMAX), & CLST_A_DIFF(NSMAX), & CDST_U_DIFF(NSMAX,NUMAX), & CYST_U_DIFF(NSMAX,NUMAX), @@ -548,33 +296,13 @@ C & CDST_D_DIFF(NSMAX,NDMAX), & CYST_D_DIFF(NSMAX,NDMAX), & CLST_D_DIFF(NSMAX,NDMAX), - & CDST_G_DIFF(NSMAX,NGMAX), - & CYST_G_DIFF(NSMAX,NGMAX), - & CLST_G_DIFF(NSMAX,NGMAX), & CFSTRP_DIFF(3,NSMAX), & CFST_U_DIFF(3,NSMAX,NUMAX), & CFST_D_DIFF(3,NSMAX,NDMAX), - & CFST_G_DIFF(3,NSMAX,NGMAX), & CMSTRP_DIFF(3,NSMAX), & CMST_U_DIFF(3,NSMAX,NUMAX), & CMST_D_DIFF(3,NSMAX,NDMAX), - & CMST_G_DIFF(3,NSMAX,NGMAX), - & CF_LSTRP_DIFF(3,NSMAX), - & CM_LSTRP_DIFF(3,NSMAX), - & CN_LSTRP_DIFF(NSMAX), - & CA_LSTRP_DIFF(NSMAX), - & CD_LSTRP_DIFF(NSMAX), - & CL_LSTRP_DIFF(NSMAX), - & CDV_LSTRP_DIFF(NSMAX), - & CLT_LSTRP_DIFF(NSMAX), - & CLA_LSTRP_DIFF(NSMAX), - & CMC4_LSTRP_DIFF(NSMAX), - & CMLE_LSTRP_DIFF(NSMAX), - & CNC_DIFF(NSMAX), - & DWWAKE_DIFF(NSMAX), - & CNC_U_DIFF(NSMAX,NUMAX), - & CNC_D_DIFF(NSMAX,NDMAX), - & CNC_G_DIFF(NSMAX,NGMAX) + & CDV_LSTRP_DIFF(NSMAX) C real(kind=avl_real) RV1_DIFF real(kind=avl_real) RV2_DIFF @@ -582,42 +310,23 @@ C real(kind=avl_real) RC_DIFF real(kind=avl_real) RS_DIFF real(kind=avl_real) RL_DIFF - real(kind=avl_real) RADL_DIFF real(kind=avl_real) DXV_DIFF - real(kind=avl_real) DXSTRPV_DIFF real(kind=avl_real) CHORDV_DIFF real(kind=avl_real) SLOPEV_DIFF real(kind=avl_real) SLOPEC_DIFF real(kind=avl_real) DCONTROL_DIFF real(kind=avl_real) VHINGE_DIFF - real(kind=avl_real) PHINGE_DIFF - real(kind=avl_real) VREFL_DIFF real(kind=avl_real) ENC_DIFF real(kind=avl_real) ENV_DIFF real(kind=avl_real) ENC_D_DIFF - real(kind=avl_real) ENC_G_DIFF - real(kind=avl_real) ENV_D_DIFF - real(kind=avl_real) ENV_G_DIFF - real(kind=avl_real) DCP_DIFF - real(kind=avl_real) DCP_U_DIFF - real(kind=avl_real) DCP_D_DIFF - real(kind=avl_real) DCP_G_DIFF real(kind=avl_real) GAM_DIFF - real(kind=avl_real) GAM_U_0_DIFF - real(kind=avl_real) GAM_U_D_DIFF - real(kind=avl_real) GAM_U_G_DIFF real(kind=avl_real) GAM_U_DIFF real(kind=avl_real) GAM_D_DIFF - real(kind=avl_real) GAM_G_DIFF real(kind=avl_real) SRC_DIFF - real(kind=avl_real) DBL_DIFF real(kind=avl_real) SRC_U_DIFF real(kind=avl_real) DBL_U_DIFF - real(kind=avl_real) WCSRD_DIFF real(kind=avl_real) WCSRD_U_DIFF real(kind=avl_real) WVSRD_DIFF - real(kind=avl_real) WVSRD_U_DIFF - real(kind=avl_real) CPT_DIFF real(kind=avl_real) RHS_DIFF real(kind=avl_real) RHS_U_DIFF real(kind=avl_real) RES_DIFF @@ -630,117 +339,54 @@ C & RC_DIFF(3,NVMAX), & RS_DIFF(3,NVMAX), & RL_DIFF(3,NLMAX), - & RADL_DIFF(NLMAX), & DXV_DIFF(NVMAX), - & DXSTRPV_DIFF(NVMAX), & CHORDV_DIFF(NVMAX), & SLOPEV_DIFF(NVMAX), & SLOPEC_DIFF(NVMAX), & DCONTROL_DIFF(NVMAX,NDMAX), & VHINGE_DIFF(3,NSMAX,NDMAX), - & PHINGE_DIFF(3,NSMAX,NDMAX), - & VREFL_DIFF(NSMAX,NDMAX), & ENC_DIFF(3,NVMAX), & ENV_DIFF(3,NVMAX), & ENC_D_DIFF(3,NVMAX,NDMAX), - & ENC_G_DIFF(3,NVMAX,NGMAX), - & ENV_D_DIFF(3,NVMAX,NDMAX), - & ENV_G_DIFF(3,NVMAX,NGMAX), - & DCP_DIFF(NVMAX), - & DCP_U_DIFF(NVMAX,NUMAX), - & DCP_D_DIFF(NVMAX,NDMAX), - & DCP_G_DIFF(NVMAX,NGMAX), & GAM_DIFF(NVMAX), - & GAM_U_0_DIFF(NVMAX,NUMAX), - & GAM_U_D_DIFF(NVMAX,NUMAX,NDMAX), - & GAM_U_G_DIFF(NVMAX,NUMAX,NGMAX), & GAM_U_DIFF(NVMAX,NUMAX), & GAM_D_DIFF(NVMAX,NDMAX), - & GAM_G_DIFF(NVMAX,NGMAX), & SRC_DIFF(NLMAX), - & DBL_DIFF(3,NLMAX), & SRC_U_DIFF(NLMAX,NUMAX), & DBL_U_DIFF(3,NLMAX,NUMAX), - & WCSRD_DIFF(3,NVMAX), & WCSRD_U_DIFF(3,NVMAX,NUMAX), & WVSRD_DIFF(3,NVMAX), - & WVSRD_U_DIFF(3,NVMAX,NUMAX), - & CPT_DIFF(NVMAX), & RHS_DIFF(NVMAX), & RHS_U_DIFF(NVMAX,NUMAX), & RES_DIFF(NVMAX), & RES_U_DIFF(NVMAX,NUMAX), & RES_D_DIFF(NVMAX,NDMAX) C - real(kind=avl_real) ELBDY_DIFF - real(kind=avl_real) SRFBDY_DIFF - real(kind=avl_real) VOLBDY_DIFF - real(kind=avl_real) DCPB_DIFF real(kind=avl_real) CDBDY_DIFF real(kind=avl_real) CYBDY_DIFF real(kind=avl_real) CLBDY_DIFF real(kind=avl_real) CFBDY_DIFF real(kind=avl_real) CMBDY_DIFF - real(kind=avl_real) CMBDYBAX_DIFF COMMON /BODY_R_DIFF/ - & ELBDY_DIFF(NBMAX), - & SRFBDY_DIFF(NBMAX), - & VOLBDY_DIFF(NBMAX), - & DCPB_DIFF(3,NLMAX), & CDBDY_DIFF(NBMAX), & CYBDY_DIFF(NBMAX), & CLBDY_DIFF(NBMAX), & CFBDY_DIFF(3,NBMAX), - & CMBDY_DIFF(3,NBMAX), - & CMBDYBAX_DIFF(3,NBMAX) + & CMBDY_DIFF(3,NBMAX) C real(kind=avl_real) AMACH_DIFF - real(kind=avl_real) VC_DIFF - real(kind=avl_real) VC_U_DIFF - real(kind=avl_real) VC_D_DIFF - real(kind=avl_real) VC_G_DIFF - real(kind=avl_real) WC_DIFF - real(kind=avl_real) WC_U_DIFF - real(kind=avl_real) WC_D_DIFF - real(kind=avl_real) WC_G_DIFF real(kind=avl_real) VV_DIFF real(kind=avl_real) VV_U_DIFF real(kind=avl_real) VV_D_DIFF - real(kind=avl_real) VV_G_DIFF real(kind=avl_real) WV_DIFF real(kind=avl_real) WV_U_DIFF real(kind=avl_real) WV_D_DIFF - real(kind=avl_real) WV_G_DIFF COMMON /SOLV_R_DIFF/ & AMACH_DIFF, - & VC_DIFF(3,NVMAX), - & VC_U_DIFF(3,NVMAX,NUMAX), - & VC_D_DIFF(3,NVMAX,NDMAX), - & VC_G_DIFF(3,NVMAX,NGMAX), - & WC_DIFF(3,NVMAX), - & WC_U_DIFF(3,NVMAX,NUMAX), - & WC_D_DIFF(3,NVMAX,NDMAX), - & WC_G_DIFF(3,NVMAX,NGMAX), & VV_DIFF(3,NVMAX), & VV_U_DIFF(3,NVMAX,NUMAX), & VV_D_DIFF(3,NVMAX,NDMAX), - & VV_G_DIFF(3,NVMAX,NGMAX), & WV_DIFF(3,NVMAX), & WV_U_DIFF(3,NVMAX,NUMAX), - & WV_D_DIFF(3,NVMAX,NDMAX), - & WV_G_DIFF(3,NVMAX,NGMAX) -C - real(kind=avl_real) ROB_DIFF - real(kind=avl_real) VOB_DIFF - real(kind=avl_real) WOB_DIFF - real(kind=avl_real) VOB_GAM_DIFF - real(kind=avl_real) WOBSRD_DIFF - real(kind=avl_real) WOBSRD_U_DIFF - COMMON /OFFBODY_R_DIFF/ - & ROB_DIFF(3,NOBMAX), - & VOB_DIFF(3,NOBMAX), - & WOB_DIFF(3,NOBMAX), - & VOB_GAM_DIFF(3,NOBMAX,NVMAX), - & WOBSRD_DIFF(3,NOBMAX), - & WOBSRD_U_DIFF(3,NOBMAX,NUMAX) + & WV_D_DIFF(3,NVMAX,NDMAX) C diff --git a/src/includes/AVL_oml.INC b/src/includes/AVL_oml.INC new file mode 100644 index 00000000..6cbc72aa --- /dev/null +++ b/src/includes/AVL_oml.INC @@ -0,0 +1,16 @@ +cc#ifdef USE_CPOML + REAL(kind=avl_real) XYN1 + REAL(kind=avl_real) XYN2 + REAL(kind=avl_real) ZLON1 + REAL(kind=avl_real) ZLON2 + REAL(kind=avl_real) ZUPN1 + REAL(kind=avl_real) ZUPN2 + REAL(kind=avl_real) XYZSURF + REAL(kind=avl_real) CPSURF + COMMON /VRTX_S/ + & XYN1(2,NVMAX), XYN2(2,NVMAX), ! left/right aft-node of element + & ZLON1(NVMAX), ZLON2(NVMAX), ! left/right lower z-coord of aft-node + & ZUPN1(NVMAX), ZUPN2(NVMAX), ! left/right upper z-coord of aft-node + & XYZSURF(3, NVMAX, NFMAX), + & CPSURF(NVMAX,NFMAX) +cc#endif \ No newline at end of file diff --git a/src/includes/AVL_surf.INC b/src/includes/AVL_surf.INC new file mode 100644 index 00000000..ffefd190 --- /dev/null +++ b/src/includes/AVL_surf.INC @@ -0,0 +1,137 @@ +C !start added variables for python geometry minipulation + + LOGICAL LDUPL_B + COMMON /BODY_GEOM_L/ + & LDUPL_B(NBMAX) ! T if body is a duplicated +c + COMMON /BODY_GEOM_I/ + & NSEC_B(NBMAX), ! number of sections in body + & NVB(NFMAX), ! number of body elements + & NBOD(NBMAX) ! number of body oml points + + REAL(kind=avl_real) XYZSCAL_B + REAL(kind=avl_real) XYZTRAN_B + REAL(kind=avl_real) YDUPL_B + REAL(kind=avl_real) XYZLES_B + REAL(kind=avl_real) BSPACE + REAL(kind=avl_real) XBOD + REAL(kind=avl_real) YBOD + REAL(kind=avl_real) TBOD + REAL(kind=avl_real) XBOD_R + REAL(kind=avl_real) YBOD_R + COMMON /BODY_GEOM_R/ + & XYZSCAL_B(3, NBMAX), ! scaling factors for XYZ coordinates + & XYZTRAN_B(3, NBMAX), ! translation factors for XYZ coordinates + & YDUPL_B(NBMAX), ! y duplicate + & XYZLES_B(3, NSMAX, NBMAX), ! leading edge cordinate vector + & BSPACE(NBMAX), ! body spacing + & XBOD(IBX, NBMAX), ! body oml x-coordinates + & YBOD(IBX, NBMAX), ! body oml y-coordinates + & TBOD(IBX, NBMAX), ! body oml thickness + & XBOD_R(IBX, NBMAX), ! raw input oml x-coordinates + & YBOD_R(IBX, NBMAX) ! raw input oml y-coordinates +c + LOGICAL LDUPL + LOGICAL LSURFSPACING + COMMON /SURF_GEOM_L/ + & LDUPL(NFMAX), ! T if surface is a duplicated + & LSURFSPACING(NFMAX) ! surface spacing set under the surface heeading + + COMMON /SURF_GEOM_I/ + & NSEC(NFMAX), ! number of sections in surface + & NVC(NFMAX), ! number of chordwise elements + & NVS(NFMAX), ! number of spanwise elements + & NSPANS(NSECMAX, NFMAX), ! number of spanwise elements vector + & NASEC(NSECMAX, NFMAX), ! number of points used to represent section geometry + & NAPTSSEC(NSECMAX, IBX), ! number of points of input airfoil geometry + & ICONTD(ICONX, NSECMAX, NFMAX), ! control variable index + & NSCON(NSECMAX, NFMAX), ! number of control variables + & IDESTD(ICONX, NSECMAX, NFMAX), ! design variable index + & NSDES(NSECMAX, NFMAX) ! number of design variables + + REAL(kind=avl_real) XYZSCAL + REAL(kind=avl_real) XYZTRAN + REAL(kind=avl_real) YDUPL + REAL(kind=avl_real) ADDINC + REAL(kind=avl_real) CSPACE + REAL(kind=avl_real) SSPACE + REAL(kind=avl_real) SSPACES + REAL(kind=avl_real) XYZLES + REAL(kind=avl_real) XSEC + REAL(kind=avl_real) YSEC + REAL(kind=avl_real) XFMIN_R + REAL(kind=avl_real) XFMAX_R +c REAL(kind=avl_real) XLASEC +c REAL(kind=avl_real) ZLASEC +c REAL(kind=avl_real) XUASEC +c REAL(kind=avl_real) ZUASEC + REAL(kind=avl_real) CHORDS + REAL(kind=avl_real) AINCS + REAL(kind=avl_real) XASEC + REAL(kind=avl_real) CASEC + REAL(kind=avl_real) SASEC + REAL(kind=avl_real) TASEC + REAL(kind=avl_real) CLCDSEC + REAL(kind=avl_real) CLCDSRF + REAL(kind=avl_real) CLAF + REAL(kind=avl_real) XHINGED + REAL(kind=avl_real) VHINGED + REAL(kind=avl_real) GAIND + REAL(kind=avl_real) REFLD + REAL(kind=avl_real) GAING + COMMON /SURF_GEOM_R/ + & XYZSCAL(3, NFMAX), ! scaling factors for XYZ coordinates + & XYZTRAN(3, NFMAX), ! translation factors for XYZ coordinates + & YDUPL(NFMAX), ! y position of duplicate plane + & ADDINC(NFMAX), ! additional incidence angle + & CSPACE(NFMAX), ! chordwise spacing + & SSPACE(NFMAX), ! spanwise spacing + & SSPACES(NSECMAX, NFMAX), ! spanwise spacing vector + & XYZLES(3, NSECMAX, NFMAX), ! leading edge cordinate vector + & XSEC(IBX, NSECMAX, NFMAX), ! raw section input x coordinates + & YSEC(IBX, NSECMAX, NFMAX), ! raw section input y coordinates + & XFMIN_R(NSECMAX, NFMAX), ! section input lower x bound + & XFMAX_R(NSECMAX, NFMAX), ! section input upper x bound +c & XLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface x-coords +c & ZLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface z-coords +c & XUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface x-coords +c & ZUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface z-coords + & CHORDS(NSECMAX, NFMAX), ! chord length vectorm + & AINCS(NSECMAX, NFMAX), ! incidence angle vector + & XASEC(NASMAX, NSECMAX, NFMAX), ! the x coordinate aifoil section + & CASEC(NASMAX, NSECMAX, NFMAX), ! camber line at xasec + & SASEC(NASMAX, NSECMAX, NFMAX), ! slope of camber line at xasec + & TASEC(NASMAX, NSECMAX, NFMAX), ! thickness at xasec + & CLCDSEC(6, NSECMAX, NFMAX), + & CLCDSRF(6, NFMAX), + & CLAF(NSECMAX, NFMAX), ! CL alpha (dCL/da) scaling factor + & XHINGED(ICONX, NSECMAX, NFMAX), ! hinge location + & VHINGED(3, ICONX, NSECMAX, NFMAX), ! hinge vector + & GAIND(ICONX, NSECMAX, NFMAX), ! control surface gain + & REFLD(ICONX, NSECMAX, NFMAX), ! control surface reflection + & GAING(ICONX, NSECMAX, NFMAX) ! desgin variable gain + + COMMON /SURF_MESH_I/ + & MFRST(NFMAX), ! stores the index in the MSHBLK where each surface's mesh begins + & IPTSEC(NSMAX,NFMAX) ! stores the iptloc vector for each surface + + REAL(kind=avl_real) MSHBLK + REAL(kind=avl_real) RV1MSH + REAL(kind=avl_real) RV2MSH + REAL(kind=avl_real) RVMSH + REAL(kind=avl_real) RCMSH + COMMON /SURF_MESH_R/ + & MSHBLK(3, 4*NVMAX), ! block to store all surface meshes + & RV1MSH(3,NVMAX), ! mesh h.v. vortex left points + & RV2MSH(3,NVMAX), ! mesh h.v. vortex right points + & RVMSH(3,NVMAX), ! mesh h.v. vortex center points + & RCMSH(3,NVMAX) ! mesh h.v. control points + + LOGICAL LSURFMSH + LOGICAL LMESHFLAT + COMMON /SURF_MESH_L/ + & LSURFMSH(NFMAX), ! T if surface uses a mesh cordinates for its geometry + & LMESHFLAT(NFMAX) ! T if the surface's mesh should be flattened when computing vorticies and control points + + +C !!--- end added variables for python geometry manipulation --- diff --git a/src/includes/AVL_surf_ad_seeds.inc b/src/includes/AVL_surf_ad_seeds.inc new file mode 100644 index 00000000..168d30bc --- /dev/null +++ b/src/includes/AVL_surf_ad_seeds.inc @@ -0,0 +1,38 @@ +C automatically generated by gen_ad_inc.py +C must be precended by 'AVL_surf.INC' + real(kind=avl_real) XYZSCAL_DIFF + real(kind=avl_real) XYZTRAN_DIFF + real(kind=avl_real) ADDINC_DIFF + real(kind=avl_real) XYZLES_DIFF + real(kind=avl_real) CHORDS_DIFF + real(kind=avl_real) AINCS_DIFF + real(kind=avl_real) XASEC_DIFF + real(kind=avl_real) SASEC_DIFF + real(kind=avl_real) CLCDSEC_DIFF + real(kind=avl_real) CLCDSRF_DIFF + real(kind=avl_real) CLAF_DIFF + COMMON /SURF_GEOM_R_DIFF/ + & XYZSCAL_DIFF(3, NFMAX), + & XYZTRAN_DIFF(3, NFMAX), + & ADDINC_DIFF(NFMAX), + & XYZLES_DIFF(3, NSECMAX, NFMAX), + & CHORDS_DIFF(NSECMAX, NFMAX), + & AINCS_DIFF(NSECMAX, NFMAX), + & XASEC_DIFF(NASMAX, NSECMAX, NFMAX), + & SASEC_DIFF(NASMAX, NSECMAX, NFMAX), + & CLCDSEC_DIFF(6, NSECMAX, NFMAX), + & CLCDSRF_DIFF(6, NFMAX), + & CLAF_DIFF(NSECMAX, NFMAX) +C + real(kind=avl_real) MSHBLK_DIFF + real(kind=avl_real) RV1MSH_DIFF + real(kind=avl_real) RV2MSH_DIFF + real(kind=avl_real) RVMSH_DIFF + real(kind=avl_real) RCMSH_DIFF + COMMON /SURF_MESH_R_DIFF/ + & MSHBLK_DIFF(3, 4*NVMAX), + & RV1MSH_DIFF(3,NVMAX), + & RV2MSH_DIFF(3,NVMAX), + & RVMSH_DIFF(3,NVMAX), + & RCMSH_DIFF(3,NVMAX) +C diff --git a/src/includes/gen_ad_inc.py b/src/includes/gen_ad_inc.py index a0373b55..05c1db42 100644 --- a/src/includes/gen_ad_inc.py +++ b/src/includes/gen_ad_inc.py @@ -1,7 +1,26 @@ +import os import re -def process_fortran_include_file(input_filename, output_filename, ad_ext="_DIFF"): +def get_used_variables(ad_src_dir, ad_ext="_DIFF"): + """Scan all .f files under ad_src_dir and return the set of base variable + names (uppercased) that appear with the given ad_ext suffix.""" + used = set() + pattern = re.compile( + r'\b([A-Za-z][A-Za-z0-9_]*)' + re.escape(ad_ext) + r'\b', + re.IGNORECASE, + ) + for root, _, files in os.walk(ad_src_dir): + for fname in files: + if fname.lower().endswith('.f'): + with open(os.path.join(root, fname)) as fh: + for match in pattern.finditer(fh.read()): + used.add(match.group(1).upper()) + return used + + +def process_fortran_include_file(input_filename, output_filename, + ad_ext="_DIFF", ad_src_dir=None): variables = {} include_lines = [] @@ -53,13 +72,13 @@ def get_line(idx_line): regex = r",(?![^(]*\))" # Matches commas that are not within parentheses _vars = re.split(regex, line) - + for _var in _vars: # import pdb; pdb.set_trace() _var = _var.strip() if _var != "" and _var not in ignore_vars: variables[block_name].append(_var) - + if _var in ignore_vars: print(f"ignoring {_var}") # import pdb; pdb.set_trace() @@ -70,10 +89,26 @@ def get_line(idx_line): idx_line += 1 idx_line += 1 - + + # Filter to only variables actually used in the AD source files. + if ad_src_dir is not None: + used = get_used_variables(ad_src_dir, ad_ext) + for block in list(variables.keys()): + before = variables[block] + variables[block] = [ + v for v in before + if v.split("(")[0].strip().upper() in used + ] + removed = ( + {v.split("(")[0].strip() for v in before} + - {v.split("(")[0].strip() for v in variables[block]} + ) + for r in sorted(removed): + print(f" skipping unused variable: {r}") + with open(output_filename, "w") as output_file: output_file.write("C automatically generated by gen_ad_inc.py\n") - output_file.write(f"C must be precended by 'include AVL_ad_seed'\n") + output_file.write(f"C must be precended by '{input_filename}'\n") # add the include lines # for line in include_lines: @@ -81,6 +116,8 @@ def get_line(idx_line): # output_file.write(f"\n\n") for block in variables: + if not variables[block]: + continue # declare the type of the variables in the block for idx_var, var in enumerate(variables[block]): if "(" in var: @@ -115,4 +152,6 @@ def get_line(idx_line): if __name__ == "__main__": - process_fortran_include_file("AVL.INC", "AVL_ad_seeds.inc") + AD_SRC = "../ad_src" + process_fortran_include_file("AVL.INC", "AVL_ad_seeds.inc", ad_src_dir=AD_SRC) + process_fortran_include_file("AVL_surf.INC", "AVL_surf_ad_seeds.inc", ad_src_dir=AD_SRC) diff --git a/src/includes/parse_common_blocks.py b/src/includes/parse_common_blocks.py new file mode 100644 index 00000000..a4013d2f --- /dev/null +++ b/src/includes/parse_common_blocks.py @@ -0,0 +1,825 @@ +#!/usr/bin/env python3 +""" +Parse a Fortran common-block include file and report the memory footprint +of every variable and every COMMON block. + +Parameter constants (NVMAX, NSMAX, …) are read automatically from +INCLUDE'd files (ADIMEN.INC, AINDEX.INC, etc.) and from PARAMETER +statements in the main file itself. + +Supports a multiplier per file — use FILE:COUNT syntax to indicate how +many times a file's storage is replicated. + +Usage: + python parse_common_blocks.py AVL.INC + python parse_common_blocks.py AVL.INC:30 AVL_surf.INC:4 + python parse_common_blocks.py AVL_oml.INC:4 AVL_surf.INC:4 AVL.INC:30 AVL_ad_seeds.INC:14 + python parse_common_blocks.py AVL.INC -I /path/to/includes + python parse_common_blocks.py AVL.INC -p ADIMEN.INC AINDEX.INC + python parse_common_blocks.py AVL.INC --avl-real 4 + python parse_common_blocks.py AVL.INC --set NVMAX 10000 --set NOBMAX 5000 + + Here is what I use currently to get the total common block allocation with llvm-flang + There seems to be some issue with llvm-flang that causes all the commonblocks in each file adds data independently + python parse_common_blocks.py AVL_oml.INC:4 AVL_surf.INC:4 AVL.INC:30 AVL_ad_seeds.INC:14 AVL_surf_ad_seeds.INC:2 +""" + +import argparse +import re +import sys +import math +from collections import OrderedDict +from pathlib import Path + +# ============================================================ +# Fallback defaults — used only when an include file cannot be +# found or a parameter is referenced but never defined. +# ============================================================ +FALLBACK_PARAMS: dict[str, int] = { + "NVMAX": 5000, + "NSMAX": 500, + "NSECMAX": 301, + "NFMAX": 100, + "NLMAX": 500, + "NBMAX": 20, + + "NUMAX": 6, + "NDMAX": 30, + "NGMAX": 21, + + "NRMAX": 25, + "NTMAX": 503, + + "NOBMAX": 1, + + "ICONX": 20, + "IBX": 200, + + "IVTOT": 20, + "ICTOT": 20, + "IPTOT": 20, + "JETOT": 12, + "KPTOT": 20, + "KPVTOT": 10, +} + + +# ============================================================ +# Type → bytes mapping +# ============================================================ + +def make_type_bytes(avl_real_bytes: int) -> dict[str, int]: + return { + "integer": 4, + "real": avl_real_bytes, + "real*4": 4, + "real*8": 8, + "real(kind=avl_real)": avl_real_bytes, + "double precision": 8, + "logical": 4, + "complex": 2 * avl_real_bytes, + "complex*8": 8, + "complex*16": 16, + "complex(kind=avl_real)": 2 * avl_real_bytes, + } + + +# ============================================================ +# Helpers +# ============================================================ + +def eval_dim(expr: str, params: dict[str, int]) -> int: + """Safely evaluate a dimension expression using *params*.""" + expr = expr.strip() + safe = re.sub(r"[A-Za-z_]\w*", lambda m: str(params.get(m.group(), m.group())), expr) + try: + return int(eval(safe, {"__builtins__": {}})) + except Exception as e: + print(f" *** WARNING: cannot evaluate dimension '{expr}' -> '{safe}': {e}", + file=sys.stderr) + return 0 + + +def human_bytes(n: int) -> str: + """Return a human-readable byte string.""" + n = int(n) + if n == 0: + return "0 B" + units = ["B", "KiB", "MiB", "GiB", "TiB"] + idx = min(int(math.log2(n) // 10), len(units) - 1) if n > 0 else 0 + value = n / (1 << (10 * idx)) + if idx == 0: + return f"{n} B" + return f"{value:,.2f} {units[idx]}" + + +# ============================================================ +# Fortran source reader — joins continuation lines +# ============================================================ + +def read_fortran_lines(path: str) -> list[str]: + """Read a fixed-form Fortran file, join continuation lines, + strip comments, and return a list of logical lines.""" + raw = Path(path).read_text(errors="replace").splitlines() + logical: list[str] = [] + + for line in raw: + # Completely blank or comment line (C/c/*/! in column 1) + if not line or line[0] in ("C", "c", "!", "*"): + continue + + # Inline comment after '!' + bang = line.find("!") + if bang > 0: + line = line[:bang] + + # Fixed-form continuation: non-blank, non-zero char in column 6 + if len(line) > 5 and line[0] == " " and line[5] not in (" ", "0", "\t"): + if logical: + logical[-1] += " " + line[6:].strip() + continue + + # Tab-form continuation + if line and line[0] == "\t" and len(line) > 1 and line[1].isdigit(): + if logical: + logical[-1] += " " + line[2:].strip() + continue + + logical.append(line.rstrip()) + + return logical + + +# ============================================================ +# Parse PARAMETER statements from Fortran source +# ============================================================ + +_RE_PARAMETER = re.compile( + r"^\s*PARAMETER\s*\((.+)\)\s*$", re.IGNORECASE +) + +def parse_parameters_from_lines(lines: list[str], params: dict[str, int]): + for line in lines: + m = _RE_PARAMETER.match(line) + if not m: + continue + + body = m.group(1) + depth = 0 + parts: list[str] = [] + cur: list[str] = [] + for ch in body: + if ch == "(": + depth += 1; cur.append(ch) + elif ch == ")": + depth -= 1; cur.append(ch) + elif ch == "," and depth == 0: + parts.append("".join(cur).strip()); cur = [] + else: + cur.append(ch) + parts.append("".join(cur).strip()) + + for part in parts: + if "=" not in part: + continue + name, expr = part.split("=", 1) + name = name.strip().upper() + expr = expr.strip() + val = eval_dim(expr, params) + if val != 0 or expr.strip() == "0": + params[name] = val + + +def parse_parameters_from_file(path: str, params: dict[str, int]): + lines = read_fortran_lines(path) + parse_parameters_from_lines(lines, params) + + +# ============================================================ +# Resolve INCLUDE directives +# ============================================================ + +_RE_INCLUDE = re.compile( + r"""^\s*INCLUDE\s+['"]([^'"]+)['"]\s*$""", re.IGNORECASE +) + +def find_include_file(name: str, search_dirs: list[Path]) -> Path: + for d in search_dirs: + candidate = d / name + if candidate.is_file(): + return candidate + return None + + +def resolve_includes( + main_path: str, + extra_search_dirs: list[str] = None, +) -> list[str]: + main_dir = Path(main_path).resolve().parent + search_dirs = [main_dir] + if extra_search_dirs: + search_dirs.extend(Path(d).resolve() for d in extra_search_dirs) + + raw = Path(main_path).read_text(errors="replace").splitlines() + found: list[str] = [] + not_found: list[str] = [] + + for line in raw: + stripped = line.lstrip() + if stripped and stripped[0] in ("C", "c", "!", "*"): + continue + m = _RE_INCLUDE.match(line) + if m: + inc_name = m.group(1) + inc_path = find_include_file(inc_name, search_dirs) + if inc_path: + found.append(str(inc_path)) + else: + not_found.append(inc_name) + + return found, not_found + + +# ============================================================ +# Parse type declarations +# ============================================================ + +_RE_TYPE_DECL = re.compile( + r""" + ^\s* + (?: + (INTEGER|LOGICAL) + | (REAL|COMPLEX)\s*\(\s*kind\s*=\s*avl_real\s*\) + | (REAL\*\d+|COMPLEX\*\d+|DOUBLE\s+PRECISION) + | (REAL|COMPLEX|INTEGER|LOGICAL) + | (CHARACTER)\*(\d+) + ) + \s+ + (.+) + """, + re.IGNORECASE | re.VERBOSE, +) + +def parse_type_declarations(lines: list[str], type_bytes: dict[str, int], + avl_real_bytes: int) -> dict[str, tuple[str, int]]: + var_types: dict[str, tuple[str, int]] = {} + + for line in lines: + m = _RE_TYPE_DECL.match(line) + if not m: + continue + + if m.group(1): + tkey = m.group(1).upper().strip() + elif m.group(2): + tkey = m.group(2).lower() + "(kind=avl_real)" + elif m.group(3): + tkey = re.sub(r"\s+", " ", m.group(3)).lower().strip() + elif m.group(4): + tkey = m.group(4).lower().strip() + elif m.group(5): + tkey = f"character*{m.group(6)}" + else: + continue + + bpe = type_bytes.get(tkey.lower()) + if bpe is None: + cm = re.match(r"character\*(\d+)", tkey, re.IGNORECASE) + if cm: + bpe = int(cm.group(1)) + else: + bpe = avl_real_bytes + + names_str = m.group(7) if m.group(7) else "" + if not names_str: + continue + depth = 0 + parts: list[str] = [] + cur: list[str] = [] + for ch in names_str: + if ch == "(": + depth += 1; cur.append(ch) + elif ch == ")": + depth -= 1; cur.append(ch) + elif ch == "," and depth == 0: + parts.append("".join(cur).strip()); cur = [] + else: + cur.append(ch) + parts.append("".join(cur).strip()) + + for part in parts: + name = re.split(r"[(\s]", part)[0].strip().upper() + if name: + var_types[name] = (tkey, bpe) + + return var_types + + +# ============================================================ +# Parse CHARACTER declarations (old-style) +# ============================================================ + +_RE_CHAR_DECL = re.compile( + r"^\s*CHARACTER\*(\d+)\s+(.+)", re.IGNORECASE +) + +def parse_character_decls(lines: list[str], var_types: dict[str, tuple[str, int]]): + for line in lines: + m = _RE_CHAR_DECL.match(line) + if not m: + continue + charlen = int(m.group(1)) + names_str = m.group(2) + depth = 0 + parts: list[str] = [] + cur: list[str] = [] + for ch in names_str: + if ch == "(": + depth += 1; cur.append(ch) + elif ch == ")": + depth -= 1; cur.append(ch) + elif ch == "," and depth == 0: + parts.append("".join(cur).strip()); cur = [] + else: + cur.append(ch) + parts.append("".join(cur).strip()) + for part in parts: + name = re.split(r"[(\s]", part)[0].strip().upper() + if name: + var_types[name] = (f"character*{charlen}", charlen) + + +# ============================================================ +# Parse COMMON blocks +# ============================================================ + +_RE_COMMON = re.compile(r"^\s*COMMON\s*/\s*(\w+)\s*/\s*(.+)", re.IGNORECASE) + +def parse_common_blocks(lines: list[str]) -> OrderedDict[str, list[tuple[str, str]]]: + blocks: OrderedDict[str, list[tuple[str, str]]] = OrderedDict() + current_block: str = None + current_tail = "" + + for line in lines: + m = _RE_COMMON.match(line) + if m: + if current_block is not None: + _parse_var_list(current_tail, blocks[current_block]) + current_block = m.group(1).upper() + if current_block not in blocks: + blocks[current_block] = [] + current_tail = m.group(2) + continue + + if current_block is not None: + _parse_var_list(current_tail, blocks[current_block]) + + return blocks + + +def _parse_var_list(text: str, out: list[tuple[str, str]]): + text = re.sub(r"!.*", "", text) + depth = 0 + parts: list[str] = [] + cur: list[str] = [] + for ch in text: + if ch == "(": + depth += 1; cur.append(ch) + elif ch == ")": + depth -= 1; cur.append(ch) + elif ch == "," and depth == 0: + parts.append("".join(cur).strip()); cur = [] + else: + cur.append(ch) + tail = "".join(cur).strip() + if tail: + parts.append(tail) + + for part in parts: + part = part.strip() + if not part: + continue + pm = re.match(r"(\w+)\s*\(([^)]*)\)", part) + if pm: + out.append((pm.group(1).upper(), pm.group(2))) + else: + name = part.strip().upper() + if name and re.match(r"[A-Z_]\w*$", name): + out.append((name, "")) + + +# ============================================================ +# Compute sizes +# ============================================================ + +def compute_sizes( + blocks: OrderedDict[str, list[tuple[str, str]]], + var_types: dict[str, tuple[str, int]], + params: dict[str, int], + avl_real_bytes: int, +) -> OrderedDict[str, list[dict]]: + result: OrderedDict[str, list[dict]] = OrderedDict() + + for bname, varlist in blocks.items(): + entries: list[dict] = [] + for vname, dims_str in varlist: + tinfo = var_types.get(vname) + if tinfo: + tkey, bpe = tinfo + else: + if bname.endswith("_I"): + tkey, bpe = "integer", 4 + elif bname.endswith("_L"): + tkey, bpe = "logical", 4 + else: + tkey, bpe = "real(kind=avl_real)", avl_real_bytes + + if dims_str: + dim_parts = [] + depth = 0 + cur: list[str] = [] + for ch in dims_str: + if ch == "(": + depth += 1; cur.append(ch) + elif ch == ")": + depth -= 1; cur.append(ch) + elif ch == "," and depth == 0: + dim_parts.append("".join(cur).strip()); cur = [] + else: + cur.append(ch) + dim_parts.append("".join(cur).strip()) + + nelems = 1 + dim_vals: list[int] = [] + for d in dim_parts: + v = eval_dim(d, params) + dim_vals.append(v) + nelems *= v + else: + nelems = 1 + dim_vals = [] + + total_bytes = nelems * bpe + entries.append({ + "name": vname, + "type": tkey, + "bpe": bpe, + "dims_str": dims_str, + "dim_vals": dim_vals, + "nelems": nelems, + "bytes": total_bytes, + }) + + result[bname] = entries + + return result + + +# ============================================================ +# Pretty-print — weighted report +# ============================================================ + +def print_weighted_report( + file_results: list[dict], + params: dict[str, int], +): + """Print a combined report across all files, weighted by their counts. + + file_results is a list of dicts: + { + "path": str, + "count": int, + "result": OrderedDict[block_name -> [entries]], + } + """ + + # ── Per-file summary ── + print(f"\n{'='*90}") + print(f" PER-FILE SUMMARY") + print(f"{'='*90}") + print(f" {'File':<30s} {'Count':>6s} {'Single copy':>14s} {'Total (× count)':>18s}") + print(f" {'-'*30} {'-'*6} {'-'*14} {'-'*18}") + + overall_total = 0 + for fr in file_results: + single = sum(e["bytes"] for entries in fr["result"].values() for e in entries) + weighted = single * fr["count"] + overall_total += weighted + print(f" {Path(fr['path']).name:<30s} {fr['count']:>6d} {human_bytes(single):>14s} {human_bytes(weighted):>18s}") + + print(f" {'-'*30} {'-'*6} {'-'*14} {'-'*18}") + print(f" {'GRAND TOTAL':<30s} {'':>6s} {'':>14s} {human_bytes(overall_total):>18s}") + + # ── Per-file block breakdown ── + for fr in file_results: + fname = Path(fr["path"]).name + count = fr["count"] + result = fr["result"] + file_single = sum(e["bytes"] for entries in result.values() for e in entries) + file_total = file_single * count + + print(f"\n{'='*90}") + print(f" COMMON BLOCKS IN: {fname} (×{count}) single={human_bytes(file_single)} total={human_bytes(file_total)}") + print(f"{'='*90}") + print(f" {'Block':<30s} {'Single copy':>14s} {'× {:<4d}'.format(count):>12s} {'% of grand':>10s}") + print(f" {'-'*30} {'-'*14} {'-'*12} {'-'*10}") + + block_list = [] + for bname, entries in result.items(): + bsize = sum(e["bytes"] for e in entries) + block_list.append((bname, bsize)) + + block_list.sort(key=lambda x: -x[1]) + for bname, bsize in block_list: + weighted = bsize * count + pct = 100.0 * weighted / overall_total if overall_total else 0 + print(f" {bname:<30s} {human_bytes(bsize):>14s} {human_bytes(weighted):>12s} {pct:>9.1f}%") + + print(f" {'-'*30} {'-'*14} {'-'*12} {'-'*10}") + pct_file = 100.0 * file_total / overall_total if overall_total else 0 + print(f" {'FILE TOTAL':<30s} {human_bytes(file_single):>14s} {human_bytes(file_total):>12s} {pct_file:>9.1f}%") + + # ── Weighted COMMON block ranking (across all files) ── + # A block name may appear in multiple files; aggregate weighted sizes. + block_weighted: dict[str, int] = {} + block_files: dict[str, list[str]] = {} + for fr in file_results: + fname = Path(fr["path"]).name + count = fr["count"] + for bname, entries in fr["result"].items(): + bsize = sum(e["bytes"] for e in entries) + weighted = bsize * count + block_weighted[bname] = block_weighted.get(bname, 0) + weighted + block_files.setdefault(bname, []).append(f"{fname}×{count}") + + sorted_blocks = sorted(block_weighted.items(), key=lambda x: -x[1]) + + print(f"\n{'='*90}") + print(f" WEIGHTED COMMON BLOCK RANKING (all files combined)") + print(f"{'='*90}") + print(f" {'#':<4s} {'Block':<24s} {'Weighted size':>14s} {'% of total':>10s} {'Sources':<30s}") + print(f" {'-'*4} {'-'*24} {'-'*14} {'-'*10} {'-'*30}") + + for rank, (bname, wsize) in enumerate(sorted_blocks[:40], 1): + pct = 100.0 * wsize / overall_total if overall_total else 0 + sources = ", ".join(block_files[bname]) + print(f" {rank:<4d} {bname:<24s} {human_bytes(wsize):>14s} {pct:>9.1f}% {sources:<30s}") + + # ── Top variables by weighted size ── + # Collect (file, count, block, var_entry) tuples + all_vars: list[tuple[str, int, str, dict]] = [] + for fr in file_results: + fname = Path(fr["path"]).name + count = fr["count"] + for bname, entries in fr["result"].items(): + for e in entries: + all_vars.append((fname, count, bname, e)) + + # Sort by weighted bytes (size × count) + all_vars.sort(key=lambda x: -(x[3]["bytes"] * x[1])) + top_n = all_vars[:50] + + if top_n: + print(f"\n{'='*110}") + print(f" TOP {len(top_n)} LARGEST VARIABLES (by weighted size = single size × count)") + print(f"{'='*110}") + print(f" {'#':<4s} {'Variable':<20s} {'Block':<18s} {'File':<20s} {'Cnt':>4s} {'Type':<22s} " + f"{'Single size':>12s} {'Weighted':>14s} {'% total':>8s}") + print(f" {'-'*4} {'-'*20} {'-'*18} {'-'*20} {'-'*4} {'-'*22} " + f"{'-'*12} {'-'*14} {'-'*8}") + + for rank, (fname, count, bname, e) in enumerate(top_n, 1): + weighted = e["bytes"] * count + pct = 100.0 * weighted / overall_total if overall_total else 0 + print(f" {rank:<4d} {e['name']:<20s} {bname:<18s} {fname:<20s} {count:>4d} {e['type']:<22s} " + f"{human_bytes(e['bytes']):>12s} {human_bytes(weighted):>14s} {pct:>7.1f}%") + + print(f" {'-'*4} {'-'*20} {'-'*18} {'-'*20} {'-'*4} {'-'*22} " + f"{'-'*12} {'-'*14} {'-'*8}") + top_total = sum(x[3]["bytes"] * x[1] for x in top_n) + top_pct = 100.0 * top_total / overall_total if overall_total else 0 + print(f" {'':4s} {'Top ' + str(len(top_n)) + ' total':<20s} {'':18s} {'':20s} {'':>4s} {'':22s} " + f"{'':>12s} {human_bytes(top_total):>14s} {top_pct:>7.1f}%") + + print(f" {'-'*80}") + print(f" {'GRAND TOTAL':<30s} {'':>6s} {'':>14s} {human_bytes(overall_total):>18s}") + print(f" {'-'*80}") + + + print() + + +# ============================================================ +# Main +# ============================================================ + +def parse_file_spec(spec: str) -> tuple[str, int]: + """Parse 'filename:count' or just 'filename' (count defaults to 1).""" + # Handle Windows paths with drive letters like C:\...:4 + # Try splitting from the right + parts = spec.rsplit(":", 1) + if len(parts) == 2 and parts[1].isdigit(): + return parts[0], int(parts[1]) + return spec, 1 + + +def main(): + ap = argparse.ArgumentParser( + description="Report weighted memory footprint of Fortran COMMON blocks. " + "Use FILE:COUNT syntax to specify how many times each file's " + "storage is replicated (e.g. AVL.INC:30). " + "The report shows totals weighted by count.", + epilog="Example:\n" + " python parse_common_blocks.py AVL_oml.INC:4 AVL_surf.INC:4 " + "AVL.INC:30 AVL_ad_seeds.INC:14\n", + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + ap.add_argument( + "include_files", + nargs="+", + help="Fortran include file(s) with optional :COUNT suffix " + "(e.g. AVL.INC:30). Count defaults to 1.", + ) + ap.add_argument( + "-p", "--params", + nargs="+", metavar="FILE", + help="Fortran file(s) to read PARAMETER values from.", + ) + ap.add_argument( + "-I", "--include-dir", + action="append", default=[], metavar="DIR", + help="Additional directory to search for INCLUDE'd files.", + ) + ap.add_argument( + "--avl-real", + type=int, default=8, choices=[4, 8], + help="Byte width of avl_real kind (default: 8 = double precision).", + ) + ap.add_argument( + "-d", "--details", + action="store_true", + help="Print per-variable breakdown for each COMMON block per file.", + ) + ap.add_argument( + "--set", + nargs=2, action="append", default=[], metavar=("NAME", "VALUE"), + help="Override or add a parameter value, e.g. --set NVMAX 10000.", + ) + args = ap.parse_args() + + avl_real_bytes: int = args.avl_real + type_bytes = make_type_bytes(avl_real_bytes) + + # Parse file:count specs + file_specs: list[tuple[str, int]] = [] + for spec in args.include_files: + path, count = parse_file_spec(spec) + file_specs.append((path, count)) + + input_paths = [p for p, _ in file_specs] + + for p in input_paths: + if not Path(p).is_file(): + print(f"Error: file not found: {p}", file=sys.stderr) + sys.exit(1) + + # ---------------------------------------------------------- + # 1. Build a shared parameter table + # ---------------------------------------------------------- + params: dict[str, int] = dict(FALLBACK_PARAMS) + + param_files_found: list[str] = [] + param_files_missing: list[str] = [] + + if args.params: + search_dirs: list[Path] = [] + for ip in input_paths: + d = Path(ip).resolve().parent + if d not in search_dirs: + search_dirs.append(d) + search_dirs.extend(Path(d).resolve() for d in args.include_dir) + + for pf in args.params: + p = Path(pf) + if p.is_file(): + param_files_found.append(str(p)) + else: + found_path = find_include_file(p.name, search_dirs) + if found_path: + param_files_found.append(str(found_path)) + else: + param_files_missing.append(pf) + else: + seen: set[str] = set() + for ip in input_paths: + found, not_found = resolve_includes(ip, args.include_dir) + for f in found: + if f not in seen: + param_files_found.append(f) + seen.add(f) + for nf in not_found: + if nf not in seen: + param_files_missing.append(nf) + seen.add(nf) + + for pf in param_files_found: + parse_parameters_from_file(pf, params) + + all_lines: dict[str, list[str]] = {} + for ip in input_paths: + lines = read_fortran_lines(ip) + all_lines[ip] = lines + parse_parameters_from_lines(lines, params) + + for name, value in args.set: + params[name.upper()] = int(value) + + # ---------------------------------------------------------- + # 2. Process each file independently, keeping its count + # ---------------------------------------------------------- + file_results: list[dict] = [] + + for ip, count in file_specs: + lines = all_lines[ip] + + # Types + var_types: dict[str, tuple[str, int]] = {} + parse_character_decls(lines, var_types) + var_types.update(parse_type_declarations(lines, type_bytes, avl_real_bytes)) + + # COMMON blocks + blocks = parse_common_blocks(lines) + result = compute_sizes(blocks, var_types, params, avl_real_bytes) + + file_results.append({ + "path": ip, + "count": count, + "result": result, + }) + + # ---------------------------------------------------------- + # 3. Print header + # ---------------------------------------------------------- + print(f"\nFortran Common Block Weighted Memory Report") + print(f"{'='*60}") + print(f" {'File':<30s} {'Count':>6s}") + print(f" {'-'*30} {'-'*6}") + for ip, count in file_specs: + print(f" {Path(ip).name:<30s} {count:>6d}") + print(f"avl_real = {avl_real_bytes} bytes") + + if param_files_found: + print(f"\nParameter files read:") + for pf in param_files_found: + print(f" {pf}") + + if param_files_missing: + print(f"\nWARNING — include files not found (using fallback values):") + for pf in param_files_missing: + print(f" {pf}") + + print(f"\nParameter values used:") + for k, v in sorted(params.items()): + print(f" {k:>12s} = {v:>8,d}") + + # ---------------------------------------------------------- + # 4. Per-file detailed breakdown (if requested) + # ---------------------------------------------------------- + if args.details: + for fr in file_results: + fname = Path(fr["path"]).name + count = fr["count"] + result = fr["result"] + + for bname, entries in result.items(): + block_total = sum(e["bytes"] for e in entries) + print(f"\n{'='*100}") + print(f" COMMON /{bname}/ from {fname} ×{count} " + f"(single={human_bytes(block_total)}, weighted={human_bytes(block_total * count)})") + print(f"{'='*100}") + print(f" {'Variable':<20s} {'Type':<26s} {'Dimensions':<30s} " + f"{'Elements':>12s} {'Single':>14s} {'Weighted':>14s}") + print(f" {'-'*20} {'-'*26} {'-'*30} {'-'*12} {'-'*14} {'-'*14}") + + for e in entries: + dims_display = "" + if e["dim_vals"]: + dims_display = "(" + ", ".join(str(v) for v in e["dim_vals"]) + ")" + if e["dims_str"] and len(dims_display) < 20: + dims_display += f" [{e['dims_str']}]" + + weighted = e["bytes"] * count + print(f" {e['name']:<20s} {e['type']:<26s} {dims_display:<30s} " + f"{e['nelems']:>12,d} {human_bytes(e['bytes']):>14s} {human_bytes(weighted):>14s}") + + print(f" {'':->116s}") + print(f" {'Block total':>90s}: {human_bytes(block_total):>14s} → {human_bytes(block_total * count)}") + + # ---------------------------------------------------------- + # 5. Weighted summary report + # ---------------------------------------------------------- + print_weighted_report(file_results, params) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/second.f b/src/second.f index a2daf9f1..2ddaf210 100755 --- a/src/second.f +++ b/src/second.f @@ -1,12 +1,10 @@ - SUBROUTINE SECONDS(TSEC) REAL*8 TSEC C -C...SECNDS is a real*4 function that returns seconds. -C The value is modified by subtracting the supplied argument. -C It acts as in the VMS FORTRAN Manual. +C...Returns elapsed wall-clock time in seconds +C Replacement for non-standard SECNDS intrinsic C - REAL*4 SECNDS, TIME - TIME = 0.0 - TSEC = SECNDS(TIME) - END + INTEGER COUNT, COUNT_RATE + CALL SYSTEM_CLOCK(COUNT, COUNT_RATE) + TSEC = DBLE(COUNT) / DBLE(COUNT_RATE) + END \ No newline at end of file diff --git a/tests/ovl_avl_comparison_utils.py b/tests/ovl_avl_comparison_utils.py index 89394c64..bbc29ddb 100644 --- a/tests/ovl_avl_comparison_utils.py +++ b/tests/ovl_avl_comparison_utils.py @@ -144,7 +144,7 @@ def check_vals( ) -def run_comparison(ovl, ref_data_cases, **kwargs): +def run_comparison(ovl, ref_data_cases, print_val=False, **kwargs): for ref_data in [ref_data_cases[0]]: set_inputs(ovl, ref_data) @@ -214,7 +214,7 @@ def run_comparison(ovl, ref_data_cases, **kwargs): for idx_surf, surf in enumerate(strip_data): - print(f"-------- {surf} --------") + if print_val: print(f"-------- {surf} --------") for key in strip_data[surf]: @@ -237,7 +237,7 @@ def run_comparison(ovl, ref_data_cases, **kwargs): avl_val = ref_data["outputs"]["strip_forces"][surf][avl_key] - print(key, avl_key, strip_data[surf][key][:3], avl_val[:3]) + if print_val: print(key, avl_key, strip_data[surf][key][:3], avl_val[:3]) check_vals(strip_data[surf][key], avl_val, key, **kwargs) diff --git a/tests/test_import.py b/tests/test_import.py index 71c643df..69066f12 100644 --- a/tests/test_import.py +++ b/tests/test_import.py @@ -1,15 +1,7 @@ -# ============================================================================= -# Extension modules -# ============================================================================= - -# ============================================================================= # Standard Python Modules -# ============================================================================= import os -# ============================================================================= # External Python modules -# ============================================================================= import unittest base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder @@ -25,7 +17,6 @@ def test_instances(self): from optvl import OVLSolver ovl_solver1 = OVLSolver(geo_file=geom_file) - ovl_solver2 = OVLSolver(geo_file=geom_file) assert ovl_solver1.avl is not ovl_solver2.avl diff --git a/tests/test_input_dict.py b/tests/test_input_dict.py index 477db9a5..1b8024f5 100644 --- a/tests/test_input_dict.py +++ b/tests/test_input_dict.py @@ -8,6 +8,7 @@ # ============================================================================= import os from copy import deepcopy +import copy # ============================================================================= # External Python modules @@ -539,11 +540,13 @@ def test_comp_input_file(self): coefs[i]["CL"], coefs_f[i]["CL"], rtol=1e-8, + err_msg="CL" ) np.testing.assert_allclose( coefs[i]["CD"], coefs_f[i]["CD"], rtol=1e-8, + err_msg="CD" ) else: diff --git a/tests/test_io.py b/tests/test_io.py index 877de069..98c38345 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -9,6 +9,7 @@ import os import psutil import re +import platform # ============================================================================= # External Python modules @@ -186,7 +187,7 @@ def test_get_scalar(self): def test_get_array(self): chords = self.ovl.get_avl_fort_arr("SURF_GEOM_R", "CHORDS") - self.assertEqual(chords.shape, (100, 301)) + self.assertEqual(chords.shape, (self.ovl.NFMAX, self.ovl.NSECMAX)) np.testing.assert_array_equal(chords[0, :5], np.array([0.45, 0.45, 0.4, 0.3, 0.2])) def parse_constants_file(filepath: str) -> dict[str, int]: @@ -208,16 +209,5 @@ def parse_constants_file(filepath: str) -> dict[str, int]: constants[name] = value return constants - -class TestConstants(unittest.TestCase): - def setUp(self): - self.ovl = OVLSolver(geo_file=geom_file, mass_file=mass_file) - - def test_constants(self): - # read the constants from src - constants = parse_constants_file(os.path.join(base_dir, "..", "src", "includes", "ADIMEN.INC")) - for var in constants: - assert getattr(self.ovl, var) == constants[var] - if __name__ == "__main__": unittest.main() diff --git a/tools/wheels/cibw_test_command.sh b/tools/wheels/cibw_test_command.sh index 99ebf954..e55055dc 100644 --- a/tools/wheels/cibw_test_command.sh +++ b/tools/wheels/cibw_test_command.sh @@ -3,9 +3,12 @@ set -xe PROJECT_DIR="$1" cd $PROJECT_DIR/tests +python -c "import platform; print(platform.machine())" +ls ../optvl # install tesing dependencies -pip install "scipy<=1.16" psutil "openmdao!=3.38" +pip install --only-binary :all: scipy +pip install psutil "openmdao!=3.38" #HACK: if the tests are not split up the CI runs out of memory... @@ -15,7 +18,6 @@ pip install "scipy<=1.16" psutil "openmdao!=3.38" # test package built and installed correctly -python test_import.py python -m unittest -v test_import.py python -m unittest -v test_io.py python -m unittest -v test_input_dict.py diff --git a/tools/wheels/repair_windows.sh b/tools/wheels/repair_windows.sh index 583e589f..58a638d6 100644 --- a/tools/wheels/repair_windows.sh +++ b/tools/wheels/repair_windows.sh @@ -1,32 +1,45 @@ +# from scipy + set -xe WHEEL="$1" DEST_DIR="$2" -# create a temporary directory in the destination folder and unpack the wheel -# into there -pushd $DEST_DIR -mkdir -p tmp -pushd tmp -wheel unpack $WHEEL -pushd optvl* -# To avoid DLL hell, the file name of libopenblas that's being vendored with -# the wheel has to be name-mangled. delvewheel is unable to name-mangle PYD -# containing extra data at the end of the binary, which frequently occurs when -# building with mingw. -# We therefore find each PYD in the directory structure and strip them. +# Skip the strip command based on TARGET_ARCH +# TARGET_ARCH should be set by the CI environment (e.g., ARM64, AMD64) +TARGET_ARCH="${TARGET_ARCH:-}" # Default to empty string if not set + +if [ "$TARGET_ARCH" = "ARM64" ]; then + echo "Skipping stripping for ARM64 target." + +else + echo "Performing stripping for AMD64 target." + + # create a temporary directory in the destination folder and unpack the wheel + # into there + pushd $DEST_DIR + mkdir -p tmp + pushd tmp + wheel unpack $WHEEL + pushd optvl* + + # To avoid DLL hell, the file name of libopenblas that's being vendored with + # the wheel has to be name-mangled. delvewheel is unable to name-mangle PYD + # containing extra data at the end of the binary, which frequently occurs when + # building with mingw. + # We therefore find each PYD in the directory structure and strip them. + + for f in $(find ./optvl* -name '*.pyd'); do strip $f; done -for f in $(find ./optvl* -name '*.pyd'); do strip $f; done + # now repack the wheel and overwrite the original + wheel pack . + mv -fv *.whl $WHEEL -# now repack the wheel and overwrite the original -wheel pack . -mv -fv *.whl $WHEEL + cd $DEST_DIR + rm -rf tmp +fi -cd $DEST_DIR -rm -rf tmp +delvewheel repair -w $DEST_DIR $WHEEL -# the libopenblas.dll is placed into this directory in the cibw_before_build -# script. -delvewheel repair --add-path /c/opt/openblas/openblas_dll -w $DEST_DIR $WHEEL