From 1272ca868b2782713d8e665a4179e294aa273c75 Mon Sep 17 00:00:00 2001 From: James Elliott Date: Thu, 14 May 2020 13:58:04 -0700 Subject: [PATCH 1/5] Scripts used for collecting build statistics --- commonTools/build_stats/wrapper/NMParser.py | 80 ++++++++ .../wrapper/WrapperCommandLineParser.py | 88 +++++++++ .../build_stats/wrapper/WrapperOpTimer.py | 185 ++++++++++++++++++ .../build_stats/wrapper/magic_wrapper.py | 62 ++++++ 4 files changed, 415 insertions(+) create mode 100644 commonTools/build_stats/wrapper/NMParser.py create mode 100644 commonTools/build_stats/wrapper/WrapperCommandLineParser.py create mode 100644 commonTools/build_stats/wrapper/WrapperOpTimer.py create mode 100755 commonTools/build_stats/wrapper/magic_wrapper.py diff --git a/commonTools/build_stats/wrapper/NMParser.py b/commonTools/build_stats/wrapper/NMParser.py new file mode 100644 index 000000000000..d85eb2acee73 --- /dev/null +++ b/commonTools/build_stats/wrapper/NMParser.py @@ -0,0 +1,80 @@ +""" + +Note: +Try to by python2 and python3 compliant +""" +import subprocess # spawning nm +import re # re matching +import os # line seperator + +class NMParser: + """Simple NM parser that""" + + # the values are + nm_option_csv_map = { + 'N' : 'symbol_debug', + 'p' : 'symbol_stack_unwind', + 'R' : 'symbol_ro_data_global', + 'r' : 'symbol_ro_data_local', + 'T' : 'symbol_text_global', + 't' : 'symbol_text_local', + 'u' : 'symbol_unique_global', + } + + nm_option_desc_map = { + 'N' : 'debugging symbol', + 'p' : 'stack unwind section', + 'R' : 'read only global data', + 'r' : 'read only local data', + 'T' : 'global text section', + 't' : 'local text section', + 'u' : 'unique global symbol', + } + + nm_re_type_expr = ''.join(nm_option_desc_map) + nm_re_str = r'^[a-zA-Z0-9]+\s+(?P[a-zA-Z0-9]{2,})\s+(?P[' + nm_re_type_expr + '])\s+' + nm_re = re.compile(nm_re_str) + + @staticmethod + def parse_object(filename): + """ + Simple NM parsing of an object file + Given an object file, we call nm -aS file + + Next, we parse stdout and match symbol lines corresponding to types + from nm_option_desc_map. + + Data are aggregated into a dict using the keys from nm_option_desc_map + + The keys are obtained from nm_option_desc_map and enforced inside the regex used + See nm_re_type_expr, nm_re_str, and nm_re in the static fields of this class + """ + p = subprocess.Popen(['nm', '-aS', filename], + stdout=subprocess.PIPE) + output = p.communicate()[0] + + nm_counts = dict() + + for line in output.split(os.linesep): + m = NMParser.nm_re.match(line) + if m: + nm_counts[m.group('type')] = nm_counts.get(m.group('type'), 0) + 1 + # return what we found + return nm_counts + + @staticmethod + def print_counts(nm_counts, + cvs_line=False, + csv_header=False): + for k,v in nm_counts.items(): + print("\"{key}\",{value}".format(key=NMParser.nm_option_desc_map[k], + value=v)) + @staticmethod + def get_csv_map (nm_counts): + # create a map of the form: csv_header_str : value + # loop over the csv_map, which will guarantee we always return the same columns. + # otherwise, looping over nm_counts will only return csv columns found in this specific file + # , while the wrapper needs consistent output from all files parsed + csv_map = { v : nm_counts.get(k,0) for k,v in NMParser.nm_option_csv_map.items() } + return csv_map + diff --git a/commonTools/build_stats/wrapper/WrapperCommandLineParser.py b/commonTools/build_stats/wrapper/WrapperCommandLineParser.py new file mode 100644 index 000000000000..0fa53da8205e --- /dev/null +++ b/commonTools/build_stats/wrapper/WrapperCommandLineParser.py @@ -0,0 +1,88 @@ +import os +import sys + +class WrapperCommandLineParser: + """ + Commandline parsing find any wrapper args, determine any output names + """ + def __init__(self, cmdline_args): + # if we write anything out it goes here + self.output_stats_file = '' + # if op generates an output file (-o ...) + self.op_output_file = '' + # if we perform an operation this is it + self.op = '' + # whether to gather and print a csv_banner + self.print_csv_banner = False + # whatever the op's args should be + self.op_args = [] + self.parse_cmdline_args(cmdline_args) + + def __repr__(self): + return self.lcl_print() + + def __str__(self): + return self.lcl_print() + + def lcl_print(self): + fmt_string = [ + 'output_stats_file : {output_stats_file}', + 'op : {op}', + 'op_output_file : {op_output_file}', + 'print_csv_banner : {print_csv_banner}', + + ] + return '\n'.join(fmt_string).format( + output_stats_file=self.output_stats_file, + op_output_file=self.op_output_file, + op=self.op, + print_csv_banner=self.print_csv_banner) + + + def parse_cmdline_args(self, cmdline_args): + wrapper_header_arg = '----get_header' + wrapper_op_arg_prefix = '----op=' + print_csv_header=False + have_op=False + # require that any wrapper arg be the first + try: + wrapper_arg = cmdline_args[1] + if wrapper_arg == wrapper_header_arg: + self.print_csv_banner=True + elif wrapper_arg.startswith(wrapper_op_arg_prefix): + self.op = wrapper_arg.split('=', 1)[1] + # find the output arg (will raise an exception if not found) + # we use -o blah.o or -o /path/to/blah.o or none at all + # we name the output as: blah.o.op.timing + # this will result in blah.ar.timing, blah.mpicc.timing blah.ld.timing... + short_op = os.path.basename(self.op) + self.output_stats_file = short_op + '.timing' + try: + output_idx = cmdline_args.index('-o') + self.op_output_file = cmdline_args[output_idx+1] + self.output_stats_file = self.op_output_file + '.' + self.output_stats_file + except: + pass + + else: + raise Exception('unparseable arguments') + + # remove the first 2 args (script name + wrapper arg) + self.op_args = cmdline_args[2:] + + except: + # any error and we give up + help_msg = ["Compiler wrapper:", + " Usage: wrapper ----op= [args] | ----get_header", + "", + " ----op=/path/to/compiler", + " path to the compiler we are wrapping", + " ----get_header", + " may not be combined with ----op, prints the csv_header generated", + "", + " Tool depends on finding a -o option in args", + " statistics will be written to .timing", + ] + print('\n'.join(help_msg)) + sys.exit(0) + diff --git a/commonTools/build_stats/wrapper/WrapperOpTimer.py b/commonTools/build_stats/wrapper/WrapperOpTimer.py new file mode 100644 index 000000000000..e5087ac6e130 --- /dev/null +++ b/commonTools/build_stats/wrapper/WrapperOpTimer.py @@ -0,0 +1,185 @@ +import subprocess +import csv +import os + +class WrapperOpTimer: + # the values are + usr_bin_time_csv_map = { + "E": + "elapsed_real_time_fmt", + "e": + "elapsed_real_time_sec", + "S": + "cpu_sec_kernel_mode", + "U": + "cpu_sec_user_mode", + "P": + "perc_cpu_used", + "M": + "max_resident_size_Kb", + "t": + "avg_resident_size_Kb", + "K": + "avg_total_memory_used_Kb", + "D": + "avg_size_unshared_data_area_Kb", + "p": + "avg_size_unshared_stack_area_Kb", + "X": + "avg_size_unshared_text_area_Kb", + "Z": + 'page_size_bytes', + "F": + "num_major_page_faults", + "R": + "num_minor_page_faults", + "W": + "num_swapped", + "c": + "num_involuntary_context_switch", + "w": + "num_waits", + "I": + "num_filesystem_inputs", + "O": + "num_filesystem_outputs", + "r": + "num_socket_msg_recv", + "s": + "num_socket_msg_sent", + "k": + "num_signals", + "x": + "exit_status", + } + + usr_bin_time_desc_map = { + "E": + "Elapsed real time ([h:]m:s)", + "e": + "Elapsed real time (s)", + "S": + "Total number of CPU-seconds that the process spent in kernel mode", + "U": + "Total number of CPU-seconds that the process spent in user mode", + "P": + "Percentage of the CPU that this job got", + "M": + "Maximum resident set size of the process during its lifetime (Kb)", + "t": + "(Not in tcsh.) Average resident set size of the process (Kb)", + "K": + "Average total (data+stack+text) memory use of the process (Kb)", + "D": + "Average size of unshared data area (Kb)", + "p": + "Average size of unshared stack space (Kb)", + "X": + "Average size of shared text space (Kb)", + "Z": + "System page size (bytes)", + "F": + "Number of major page faults", + "R": + "Number of minor or recoverable page faults", + "W": + "Number of times the process was swapped out of main memory", + "c": + "Number of times the process was context-switched involuntarily", + "w": + "Number of waits", + "I": + "Number of file system inputs by the process", + "O": + "Number of file system outputs by the process", + "r": + "Number of socket messages received by the process", + "s": + "Number of socket messages sent by the process", + "k": + "Number of signals delivered to the process", + "x": + "(Not in tcsh.) Exit status of the command", + } + + default_fields = [ + "e", + "M", + "K", + "D", + "X", + "F", + "R", + "W", + "w", + "c", + "S", + "U", + "P", + "I", + "O", + "r", + "s", + "k", + "x", + ] + + field_header_full = ','.join([ usr_bin_time_csv_map[f] for f in default_fields ]) + field_header_short = ','.join(default_fields) + field_arg = '--format=' + field_header_full + '\n' + ','.join([ '%{}'.format(f) for f in default_fields] ) + + @staticmethod + def time_op(op, + op_output_file, + output_stats_file, + op_args): + """ + evaluate 'op' with 'op_args', and gather stats into output_stats_file + """ + cmd = [ + '/usr/bin/time', + # '--append', + '--output=' + output_stats_file, + WrapperOpTimer.field_arg, + op ] + op_args + + # print(' '.join(cmd)) + p = subprocess.Popen(cmd) + p.communicate() + + # initializing the titles and rows list + fields = [] + csv_row = {} + + # reading csv file + with open(output_stats_file, 'r') as csvfile: + # creating a csv reader object + csvreader = csv.reader(csvfile) + + # extracting field names through first row + fields = next(csvreader) + + # extracting each data row one by one + # we effectively retain on the last row. + # it isn't clear if we should expect multiple rows per file + for row in csvreader: + csv_row = dict(zip(fields, row)) + + # markup the output + csv_row['FileSize'] = WrapperOpTimer.get_file_size(op_output_file) + csv_row['FileName'] = op_output_file + + return csv_row + + + # returns the file size in bytes + @staticmethod + def get_file_size(filename): + sz = -1 + try: + sz = os.stat(filename).st_size + except: + pass + return sz + + diff --git a/commonTools/build_stats/wrapper/magic_wrapper.py b/commonTools/build_stats/wrapper/magic_wrapper.py new file mode 100755 index 000000000000..b4d2380dd579 --- /dev/null +++ b/commonTools/build_stats/wrapper/magic_wrapper.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python +''' + +Note: +Try to by python2 and python3 compliant + +Ideally, keep this file minimal to reduce parsing overhead +Most 'work' is done in the Classes imported +''' + +# we need to run subcommands +import os +import csv +import sys +from NMParser import NMParser +from WrapperCommandLineParser import WrapperCommandLineParser +from WrapperOpTimer import WrapperOpTimer + +# given a dict of key/val pairs, write them as a CSV line +def write_csv_map(filename,csv_map): + try: + with open(filename, 'w') as csvfile: + writer = csv.DictWriter(csvfile, + fieldnames=[ k for k in csv_map ]) + writer.writeheader() + writer.writerow(csv_map) + except IOError: + print("I/O error") + +def main(cmdline_args): + # parse the command line args, find wrapper args and organize the + # the info into fields in this class + wcp = WrapperCommandLineParser(cmdline_args) + # you can 'print()' a WrapperCommandLineParser + # we could add a verbose mode (----verbose) + #print(wcp) + + # keep a dict of field : value + # first do the operation + # this must be first, as it generates the output file + csv_map = WrapperOpTimer.time_op(op=wcp.op, + op_output_file=wcp.op_output_file, + output_stats_file=wcp.output_stats_file, + op_args=wcp.op_args) + # test nm + # we probably need some to handle the case the .o isn't created + # as-us, the following will return empty dicts (we parse/validate) the output + # from NM, so we won't return garbage + nmp = NMParser.parse_object(wcp.op_output_file) + # NMParser.print_counts(nmp) + # add NM's output to our csv data + # we could move this into the NM parser + csv_map.update(NMParser.get_csv_map(nmp)) + + # ultimately, print the csv data to a file + # make sure to quote csv columns + write_csv_map(wcp.output_stats_file, csv_map) + + +if __name__ == '__main__': + main(sys.argv) + From f3bd07fba489f321411e9227064a9b235d0c698a Mon Sep 17 00:00:00 2001 From: Kyungjoo Kim Date: Thu, 14 May 2020 19:47:02 -0700 Subject: [PATCH 2/5] Tacho - Clark's request. --- .../example/Tacho_ExampleExternalInterface.cpp | 6 ++++++ .../example/Tacho_ExampleExternalInterface.hpp | 10 ++++++++++ .../shylu/shylu_node/tacho/src/Tacho_Solver.hpp | 2 ++ .../tacho/src/impl/Tacho_Solver_Impl.hpp | 14 ++++++++++++++ 4 files changed, 32 insertions(+) diff --git a/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.cpp b/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.cpp index 03fc7cd9ddd6..38f253547d57 100644 --- a/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.cpp +++ b/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.cpp @@ -80,6 +80,12 @@ void testTachoSolver(int numRows, Kokkos::fence(); + /// + /// Export supernodes + /// + std::vector supernodes; + solver.exportSupernodes(supernodes); + /// /// Export matrix and permutation vector /// diff --git a/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.hpp b/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.hpp index 3e45800cece6..98329e60b0b4 100644 --- a/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.hpp +++ b/packages/shylu/shylu_node/tacho/example/Tacho_ExampleExternalInterface.hpp @@ -117,6 +117,16 @@ namespace tacho { return 0; } + void exportSupernodes(std::vector &supernodes) { + const auto supernodes_device = m_Solver.getSupernodes(); + auto supernodes_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), supernodes_device); + { + const int n = supernodes_host.extent(0); + supernodes.resize(n); + std::copy(supernodes_host.data(), supernodes_host.data()+n, supernodes.data()); + } + } + void exportUpperTriangularFactorsToCrsMatrix(std::vector &rowBeginU, std::vector &columnsU, std::vector &valuesU, diff --git a/packages/shylu/shylu_node/tacho/src/Tacho_Solver.hpp b/packages/shylu/shylu_node/tacho/src/Tacho_Solver.hpp index e356c9397731..7c04e51b7ed2 100644 --- a/packages/shylu/shylu_node/tacho/src/Tacho_Solver.hpp +++ b/packages/shylu/shylu_node/tacho/src/Tacho_Solver.hpp @@ -171,6 +171,8 @@ namespace Tacho { /// /// get interface /// + ordinal_type getNumSupernodes() const; + ordinal_type_array getSupernodes() const; ordinal_type_array getPermutationVector() const; ordinal_type_array getInversePermutationVector() const; diff --git a/packages/shylu/shylu_node/tacho/src/impl/Tacho_Solver_Impl.hpp b/packages/shylu/shylu_node/tacho/src/impl/Tacho_Solver_Impl.hpp index 2e3aecee3ec0..d65f41f02d37 100644 --- a/packages/shylu/shylu_node/tacho/src/impl/Tacho_Solver_Impl.hpp +++ b/packages/shylu/shylu_node/tacho/src/impl/Tacho_Solver_Impl.hpp @@ -175,6 +175,20 @@ namespace Tacho { /// /// get interface /// + template + ordinal_type + Solver + ::getNumSupernodes() const { + return _nsupernodes; + } + + template + typename Solver::ordinal_type_array + Solver + ::getSupernodes() const { + return _supernodes; + } + template typename Solver::ordinal_type_array Solver From 4c16cfad8664f634729064c22df6fc85ded8da4d Mon Sep 17 00:00:00 2001 From: "Heidi K. Thornquist" Date: Thu, 14 May 2020 21:34:03 -0600 Subject: [PATCH 3/5] Fixes errors in OrthoManagerTest and enables Epetra testing The OrthoManagerTest had several typos in it, which were not found because all tests that used it were essentially disabled. The OrthoManagerTest has been fixed and all TSQR-specific tests have been placed behind an ifdef. The Epetra interface test for the OrthoManager has been enabled and the dependency on TSQR has been removed due to the change in the tester. --- .../belos/epetra/test/OrthoManager/CMakeLists.txt | 8 ++++++-- packages/belos/src/BelosOrthoManagerTest.hpp | 15 ++++++++------- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/packages/belos/epetra/test/OrthoManager/CMakeLists.txt b/packages/belos/epetra/test/OrthoManager/CMakeLists.txt index fb3707d36d66..e4e3a24536b4 100644 --- a/packages/belos/epetra/test/OrthoManager/CMakeLists.txt +++ b/packages/belos/epetra/test/OrthoManager/CMakeLists.txt @@ -2,11 +2,15 @@ ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) -IF (${PACKAGE_NAME}_ENABLE_Triutils AND ${PACKAGE_NAME}_ENABLE_TSQR) +IF (${PACKAGE_NAME}_ENABLE_Triutils) + TRIBITS_ADD_EXECUTABLE_AND_TEST( Epetra_OrthoManager_test SOURCES belos_orthomanager_epetra.cpp - ARGS "" + ARGS "--ortho=ICGS --verbose" + "--ortho=DGKS --verbose" + "--ortho=IMGS --verbose" COMM serial mpi ) + ENDIF() diff --git a/packages/belos/src/BelosOrthoManagerTest.hpp b/packages/belos/src/BelosOrthoManagerTest.hpp index 53da89c1edee..97589c149d13 100644 --- a/packages/belos/src/BelosOrthoManagerTest.hpp +++ b/packages/belos/src/BelosOrthoManagerTest.hpp @@ -463,7 +463,7 @@ namespace Belos { debugOut << "done: || || = " << err << endl; } - +#ifdef HAVE_BELOS_TSQR // // If OM is an OutOfPlaceNormalizerMixin, exercise the // out-of-place normalization routines. @@ -560,6 +560,7 @@ namespace Belos { << "=== Done with OutOfPlaceNormalizerMixin tests ===" << endl << endl; } +#endif // HAVE_BELOS_TSQR { // @@ -640,7 +641,7 @@ namespace Belos { std::vector ind(1); ind[0] = i; RCP Si = MVT::CloneViewNonConst(*S,ind); - MVT::MvAddMv(scaleS[i],*one,ZERO,*one,*Si); + MVT::MvAddMv(scaleS(i,0),*one,ZERO,*one,*Si); } debugOut << "Testing normalize() on a rank-1 multivector " << endl; const int thisNumFailed = testNormalize(OM,S,MyOM); @@ -725,7 +726,7 @@ namespace Belos { std::vector ind(1); ind[0] = i; RCP Si = MVT::CloneViewNonConst(*S,ind); - MVT::MvAddMv(scaleS[i],*one,ZERO,*one,*Si); + MVT::MvAddMv(scaleS(i,0),*one,ZERO,*one,*Si); } debugOut << "Testing projectAndNormalize() on a rank-1 multivector " << endl; bool constantStride = true; @@ -940,9 +941,9 @@ namespace Belos { // copies of S,MS Scopy = MVT::CloneCopy(*S); // randomize this data, it should be overwritten - Teuchos::randomSyncedMatrix(B); + Teuchos::randomSyncedMatrix(*B); for (size_type i=0; inormalize (*S_copy, B); sout << "normalize() returned rank " << reportedRank << endl; From 51b7fd47859f36ea15a423e13331f1671eac810d Mon Sep 17 00:00:00 2001 From: James Elliott Date: Fri, 15 May 2020 03:04:45 -0700 Subject: [PATCH 4/5] Add an ENV variable inside CMake so that tools can determine if they are being run inside cmake (e.g., detect configure vs build time) --- CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 520d6ac81d38..11219b520b39 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -72,6 +72,9 @@ PROJECT(${PROJECT_NAME} NONE) # Set default C++ standard to C++11 SET(Trilinos_ENABLE_CXX11_DEFAULT ON) +## set an env so we know we are in configure +set(ENV{CMAKE_IS_IN_CONFIGURE_MODE} 1) + # # B) Pull in the TriBITS system and execute # From 669b52d6a76bf520d239674afc6e3efbc64468cb Mon Sep 17 00:00:00 2001 From: "Irina K. Tezaur" Date: Fri, 15 May 2020 07:46:05 -0700 Subject: [PATCH 5/5] Implementing forward sensitivities for the Piro::TempusSolver class for use by Albany (#7359) * Piro: adding a couple member functions to Piro::TransientSolver class. * Piro: cleanups to Piro::TempusSolver. Removing unneeded logic/variables. Renaming member variables in a consistent way. * Piro: more cleanups of Piro::TempusSolver. Adding some unit tests involving transient sensitivities, which will fail for now since I have not written the code for the sensitivities in the Piro::TempusSolver class. * Piro: more cleanups. * Piro: more refactoring/cleanup towards adding transient sensitivities. Moving response stuff to Piro::TransientSolver. Extending code to work with >1 parameters and responses (caveat: this case has not been tested yet). * Piro: adding flags for forward/adjoint sensitivities. Further refactors. * Tempus: replacing DERIV_MV_BY_COL with DERIV_MV_JACOBIAN_FORM and DERIV_TRANS_MV_BY_ROW with DERIV_MV_GRADIENT_FORM. DERIV_MV_BY_COL and DERIV_TRANS_MV_BY_ROW are deprecated in Thyra. * Piro: adding ifdef to allow option of disabling sensivity tests. * Piro: work on transient sensitivities. * Creating Piro::TempusIntegrator class, which is effectively a wrapper class of various integrators from Tempus (IntegratorBasic, IntegratorFowardSensitivity, IntegratorAdjointSensitivity). * Piro: forgot to check this in with previous commit. * Piro: filling in TempusIntegrator routines. * Piro: switching TempusSolver to use new Piro::TempusIntegrator class. This is an interface to the various integrators in Piro, including the sensitivity ones. Working tests pass. * Piro: hooking up setting of sensitivity method based on input parameters for transient sensitivities. * Piro: continuing hooking up of relevant Tempus routines for transient sensitivities. * Piro: filling in response sensitivities code. * Piro: continuing work on transient forward sensitivities. Filling in/cleaning up/documenting response sensitivity calculation. * Piro: adding getDxDp(), getDxdotDp() and getDxdotdotDp() routines to Piro::TempusIntegrator class, as well as Piro::TempusIntegrator member variable to Piro::TransientSolver class. Populating dxdp_mv in evalConvergedModel() using getDxDp(). * Piro: adding getDgDp() routine to Piro::TempusIntegrator - looking forward to adjoint sensitivities. * Piro: correcting comment about clearObservers() routine not existing for sensitivity integrators in Tempus. * Piro: cleaning up verbose output. Removing exception throwing when requesting sensitivities. Refactoring Tempus unit tests into sensitivity ones and those w/o sensitivities. Creating an EpetraExt and Thyra version of the sensitivity tests, for debugging. Sensitivity tests currently fail - need to debug. * Piro: removing LinearOp sensitivity tests for Tempus - forward sensitivities do not work with LinearOp. * Piro: removing unneeded support. * Piro: bug fix affecting tempus solver forward sensitivity unit tests. Now the sensitivity tests all pass on 4 procs but 2 of them fail on 1 proc, which is super weird... Removing forward sensitivity unit tests using the Tpetra version of the MockModel_A, as these were redundant. * Piro: creating SinCos example, analogous to the one in Tempus, for testing forward sensitivities. Test works in serial and parallel - yay! * Piro: some cleanup and renaming of tests for forward sensitivities. * Piro: refactoring transient sensitivity sincos unit tests - creating 4 tests out of one for 'Staggered' vs. 'Combined' sensitivities, and sensitivities w/ and w/o Tangent operator. * Piro: rearranging sensitivity calculation to avoid calling getDxDp for each parameter and response. * Piro: fixing failing sensitivity tests. All tests pass once again, including new sensitivity ones :). * Piro: adding comments and optional print statements checking response sensitivities. * Piro: removing obsolete types DERIV_TRANS_MV_BY_ROW and DERIV_MV_BY_COL for Thyra model evaluator. * Piro: bug fix - sensitivity type was set too late, causing Albany transient tests to fail. * Piro: cleaning up comments. * Piro: adding back missing bracket that got deleted accidentally. * Piro: several improvements as follows. - Removing logic preventing user to start time-integration at time other than 0, which is possible to do now in Tempus. Writing corresponding tests. - Extending Piro::ObserverToTempusIntegrationObserverAdapter class so that it handles correctly sensitivities. Adding a test for sensitivities + observer. * Piro: removing __PRETTY_FUNCTION__ debug output. Fixing a couple of the errors caught by the auto-tester. * Piro: removing debug output from tests. * Piro: removing initialConditionModel argument/member variable from Tempus classes. It is not needed/used and a remnant from the Rythmos solver. * Piro: corrections to createOutArgs, which was not implemented earlier for >1 parameter and response. * Piro: changing final difference check to relative difference. * Piro: changing int sens_method to enum. * Piro: removing unused variable. * Piro: fixing Thyra type mismatch that was causing some of the auto-tester tests to fail. * Piro: renaming Piro_TempusHelpers.hpp file and removing unneeded enums * Piro: correcting include. Co-authored-by: Irina Tezaur --- packages/piro/src/CMakeLists.txt | 3 + packages/piro/src/Piro_Helpers.hpp | 52 ++ packages/piro/src/Piro_NOXSolver_Def.hpp | 14 +- ...rverToTempusIntegrationObserverAdapter.hpp | 8 +- ...ToTempusIntegrationObserverAdapter_Def.hpp | 39 +- packages/piro/src/Piro_TempusIntegrator.hpp | 120 ++++ .../piro/src/Piro_TempusIntegrator_Def.hpp | 321 +++++++++ packages/piro/src/Piro_TempusSolver.hpp | 53 +- packages/piro/src/Piro_TempusSolver_Def.hpp | 473 +++++-------- packages/piro/src/Piro_TransientSolver.hpp | 37 +- .../piro/src/Piro_TransientSolver_Def.hpp | 401 ++++++++--- .../src/TriKota_MPDirectApplicInterface.cpp | 14 +- packages/piro/test/CMakeLists.txt | 28 + packages/piro/test/MockModelEval_A_Tpetra.cpp | 6 +- ...empusSolver_SensitivitySinCosUnitTests.cpp | 311 +++++++++ ...Piro_TempusSolver_SensitivityUnitTests.cpp | 520 +++++++++++++++ .../piro/test/Piro_TempusSolver_UnitTests.cpp | 61 +- packages/piro/test/SinCosModel.cpp | 630 ++++++++++++++++++ packages/piro/test/SinCosModel.hpp | 236 +++++++ .../_input_Tempus_BackwardEuler_SinCos.xml | 110 +++ .../test/TestModels/SinCosModel_impl.hpp | 6 +- .../test/TestModels/VanDerPolModel_impl.hpp | 2 +- .../VanDerPol_IMEXPart_ImplicitModel_impl.hpp | 4 +- .../VanDerPol_IMEX_ExplicitModel_impl.hpp | 2 +- .../VanDerPol_IMEX_ImplicitModel_impl.hpp | 2 +- 25 files changed, 2968 insertions(+), 485 deletions(-) create mode 100644 packages/piro/src/Piro_Helpers.hpp create mode 100644 packages/piro/src/Piro_TempusIntegrator.hpp create mode 100644 packages/piro/src/Piro_TempusIntegrator_Def.hpp create mode 100644 packages/piro/test/Piro_TempusSolver_SensitivitySinCosUnitTests.cpp create mode 100644 packages/piro/test/Piro_TempusSolver_SensitivityUnitTests.cpp create mode 100644 packages/piro/test/SinCosModel.cpp create mode 100644 packages/piro/test/SinCosModel.hpp create mode 100644 packages/piro/test/_input_Tempus_BackwardEuler_SinCos.xml diff --git a/packages/piro/src/CMakeLists.txt b/packages/piro/src/CMakeLists.txt index 6966f1af47b4..014cbf843e78 100644 --- a/packages/piro/src/CMakeLists.txt +++ b/packages/piro/src/CMakeLists.txt @@ -130,8 +130,11 @@ ENDIF() ENDIF() IF (Piro_ENABLE_Tempus) APPEND_SET(HEADERS + Piro_Helpers.hpp Piro_TempusSolver.hpp Piro_TempusSolver_Def.hpp + Piro_TempusIntegrator.hpp + Piro_TempusIntegrator_Def.hpp Piro_TransientSolver.hpp Piro_TransientSolver_Def.hpp Piro_TempusStepperFactory.hpp diff --git a/packages/piro/src/Piro_Helpers.hpp b/packages/piro/src/Piro_Helpers.hpp new file mode 100644 index 000000000000..6b4a3ca5caf4 --- /dev/null +++ b/packages/piro/src/Piro_Helpers.hpp @@ -0,0 +1,52 @@ +// @HEADER +// ************************************************************************ +// +// Piro: Strategy package for embedded analysis capabilitites +// Copyright (2010) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Andy Salinger (agsalin@sandia.gov), Sandia +// National Laboratories. +// +// ************************************************************************ +// @HEADER + +#ifndef PIRO_HELPERS_H +#define PIRO_HELPERS_H + +namespace Piro { + + enum SENS_METHOD {NONE, FORWARD, ADJOINT}; + +} // namespace Piro + +#endif /*PIRO_HELPERS_H*/ diff --git a/packages/piro/src/Piro_NOXSolver_Def.hpp b/packages/piro/src/Piro_NOXSolver_Def.hpp index d3eef0c4632b..0df3052708a1 100644 --- a/packages/piro/src/Piro_NOXSolver_Def.hpp +++ b/packages/piro/src/Piro_NOXSolver_Def.hpp @@ -305,7 +305,7 @@ void Piro::NOXSolver::evalModelImpl( // choice to minimze the number of solves. However this choice depends // also on what layout of dg/dx is supported (e.g., if only the operator // form is supported for forward sensitivities, then df/dp must be - // DERIV_MV_BY_COL). For simplicity, we order the conditional tests + // DERIV_MV_JACOBIAN_FORM). For simplicity, we order the conditional tests // to get the right layout in most situations. DerivativeLayout dfdp_layout; { // if (sensitivity_method == "Adjoint") @@ -321,8 +321,8 @@ void Piro::NOXSolver::evalModelImpl( std::endl << "Piro::NOXSolver::evalModel(): " << "For df/dp(" << i <<") with adjoint sensitivities, " << "underlying ModelEvaluator must support DERIV_LINEAR_OP, " << - "DERIV_MV_BY_COL with p not distributed, or " - "DERIV_TRANS_MV_BY_ROW with f not distributed." << + "DERIV_MV_JACOBIAN_FORM with p not distributed, or " + "DERIV_MV_GRADIENT_FORM with f not distributed." << std::endl); } if (dfdp_layout == COL) { @@ -390,8 +390,8 @@ void Piro::NOXSolver::evalModelImpl( std::endl << "Piro::NOXSolver::evalModel(): " << "For dg/dx(" << j <<") with adjoint sensitivities, " << "underlying ModelEvaluator must support DERIV_LINEAR_OP, " << - "DERIV_MV_BY_COL with x not distributed, or " - "DERIV_TRANS_MV_BY_ROW with g not distributed." << + "DERIV_MV_JACOBIAN_FORM with x not distributed, or " + "DERIV_MV_GRADIENT_FORM with g not distributed." << std::endl); } @@ -455,8 +455,8 @@ void Piro::NOXSolver::evalModelImpl( "For dg/dp(" << j << "," << i << ") with operator sensitivities, "<< "underlying ModelEvaluator must support DERIV_LINEAR_OP, " << - "DERIV_MV_BY_COL with p not distributed, or " - "DERIV_TRANS_MV_BY_ROW with g not distributed." << + "DERIV_MV_JACOBIAN_FORM with p not distributed, or " + "DERIV_MV_GRADIENT_FORM with g not distributed." << std::endl); } else diff --git a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp index d19886318533..d610e3dc1d5a 100644 --- a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp +++ b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp @@ -45,6 +45,7 @@ #include "Tempus_IntegratorObserverBasic.hpp" +#include "Piro_Helpers.hpp" #include "Piro_ObserverBase.hpp" #include "Teuchos_RCP.hpp" @@ -64,7 +65,8 @@ class ObserverToTempusIntegrationObserverAdapter : public Tempus::IntegratorObse const Teuchos::RCP >& timeStepControl, const Teuchos::RCP > &wrappedObserver, const bool supports_x_dotdot = false, - const bool abort_on_fail_at_min_dt = false); + const bool abort_on_fail_at_min_dt = false, + const SENS_METHOD sens_method = NONE); // Overridden from Tempus::IntegratorObserver @@ -108,7 +110,9 @@ class ObserverToTempusIntegrationObserverAdapter : public Tempus::IntegratorObse Teuchos::RCP > wrappedObserver_; bool supports_x_dotdot_; Scalar previous_dt_; - bool abort_on_fail_at_min_dt_; + bool abort_on_fail_at_min_dt_; + + SENS_METHOD sens_method_; }; } // namespace Piro diff --git a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp index 084063938934..730832da5b17 100644 --- a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp +++ b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp @@ -51,16 +51,16 @@ Piro::ObserverToTempusIntegrationObserverAdapter::ObserverToTempusIntegr const Teuchos::RCP >& solutionHistory, const Teuchos::RCP >& timeStepControl, const Teuchos::RCP > &wrappedObserver, - const bool supports_x_dotdot, const bool abort_on_fail_at_min_dt) + const bool supports_x_dotdot, const bool abort_on_fail_at_min_dt, + const SENS_METHOD sens_method) : solutionHistory_(solutionHistory), timeStepControl_(timeStepControl), out_(Teuchos::VerboseObjectBase::getDefaultOStream()), wrappedObserver_(wrappedObserver), supports_x_dotdot_(supports_x_dotdot), - abort_on_fail_at_min_dt_(abort_on_fail_at_min_dt) + abort_on_fail_at_min_dt_(abort_on_fail_at_min_dt), + sens_method_(sens_method) { - //Currently, sensitivities are not supported in Tempus. - hasSensitivities_ = false; previous_dt_ = 0.0; } @@ -224,16 +224,37 @@ Piro::ObserverToTempusIntegrationObserverAdapter::observeTimeStep() } //Get solution - Teuchos::RCP > solution = solutionHistory_->getCurrentState()->getX(); - solution.assert_not_null(); - //Get solution_dot - Teuchos::RCP > solution_dot = solutionHistory_->getCurrentState()->getXDot(); + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + Teuchos::RCP> x = solutionHistory_->getCurrentState()->getX(); + Teuchos::RCP X = Teuchos::rcp_dynamic_cast(x); + //IKT, 5/12/2020: getX() returns a vector containing the sensitivities for the case of sensitivity calculations for sensitivity integrator. + //Hence, we extract the first column, which is the solution. + Teuchos::RCP> solution = (sens_method_ == NONE) ? x : X->getMultiVector()->col(0); + //IKT FIXME? Extract the remaining columns which would represent dxdp for the forward sensitivities? + //This is if we want to observe them/write them to Exodus file in Albany. + //if (sens_method == FORWARD) { + //const int num_param = X->getMultiVector()->domain()->dim()-1; + //const Teuchos::Range1D rng(1,num_param); + //Teuchos::RCP> solution_dxdp = X->getMultiVector()->subView(rng); + //} + + //Get solution_dot + Teuchos::RCP> xdot = solutionHistory_->getCurrentState()->getXDot(); + Teuchos::RCP XDot = Teuchos::rcp_dynamic_cast(xdot); + //IKT, 5/11/2020: getXDot() returns a vector containing the sensitivities for the case of sensitivity calculations for sentivity integrator + //Hence, we extract the first column, which is the solution_dot. + Teuchos::RCP> solution_dot = (sens_method_ == NONE) ? xdot : XDot->getMultiVector()->col(0); const Scalar scalar_time = solutionHistory_->getCurrentState()->getTime(); typedef typename Teuchos::ScalarTraits::magnitudeType StampScalar; const StampScalar time = Teuchos::ScalarTraits::real(scalar_time); - Teuchos::RCP > solution_dotdot = solutionHistory_->getCurrentState()->getXDotDot(); + //Get solution_dotdot + Teuchos::RCP> xdotdot = solutionHistory_->getCurrentState()->getXDotDot(); + Teuchos::RCP XDotDot = Teuchos::rcp_dynamic_cast(xdotdot); + //IKT, 5/11/2020: getXDotDot() returns a vector containing the sensitivities for the case of sensitivity calculations. + //Hence, we extract the first column, which is the solution_dotdot. + Teuchos::RCP> solution_dotdot = (sens_method_ == NONE) ? xdot : XDotDot->getMultiVector()->col(0); if (Teuchos::nonnull(solution_dot)) { if (supports_x_dotdot_) { diff --git a/packages/piro/src/Piro_TempusIntegrator.hpp b/packages/piro/src/Piro_TempusIntegrator.hpp new file mode 100644 index 000000000000..69467ccdc109 --- /dev/null +++ b/packages/piro/src/Piro_TempusIntegrator.hpp @@ -0,0 +1,120 @@ +// @HEADER +// ************************************************************************ +// +// Piro: Strategy package for embedded analysis capabilitites +// Copyright (2010) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Andy Salinger (agsalin@sandia.gov), Sandia +// National Laboratories. +// +// ************************************************************************ +// @HEADER + +#ifndef PIRO_TRANSIENTINTEGRATOR_H +#define PIRO_TRANSIENTINTEGRATOR_H + +#include "Piro_ConfigDefs.hpp" + +#include "Tempus_IntegratorBasic.hpp" +#include "Tempus_IntegratorForwardSensitivity.hpp" +#include "Tempus_IntegratorAdjointSensitivity.hpp" +#include "Tempus_StepperObserverBasic.hpp" +#include "Piro_Helpers.hpp" + +#include +#include + +namespace Piro { + +/** \brief Thyra-based Model Evaluator for Tempus solves using Tempus + * */ +template +class TempusIntegrator +{ +public: + /** \name Constructors/initializers */ + //@{ + + /** \brief . */ + TempusIntegrator(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, + const SENS_METHOD sens_method = NONE); + + Teuchos::RCP> getStepper() const; + + bool advanceTime(const Scalar time_final); + + Scalar getTime() const; + + Teuchos::RCP> getX() const; + + Teuchos::RCP> getSolutionHistory() const; + + Teuchos::RCP> getTimeStepControl() const; + + void clearObservers(); + + void setObserver(Teuchos::RCP> obs = Teuchos::null); + + void initialize(); + + void initializeSolutionHistory(Scalar t0, + Teuchos::RCP< const Thyra::VectorBase< Scalar > > x0, + Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdot0 = Teuchos::null, + Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdotdot0 = Teuchos::null, + Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxDp0 = Teuchos::null, + Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotDp0 = Teuchos::null, + Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotdotDp0 = Teuchos::null); + + Tempus::Status getStatus() const; + + // The following 3 routines are only for forward sensitivities + Teuchos::RCP> getDxDp() const; + Teuchos::RCP> getDxdotDp() const; + Teuchos::RCP> getDxdotdotDp() const; + + //The following routine is only for adjoint sensitivities + Teuchos::RCP> getDgDp() const; + +private: + + Teuchos::RCP > basicIntegrator_; + Teuchos::RCP > fwdSensIntegrator_; + Teuchos::RCP > adjSensIntegrator_; + Teuchos::RCP out_; + +}; + +} + +#include "Piro_TempusIntegrator_Def.hpp" +#endif diff --git a/packages/piro/src/Piro_TempusIntegrator_Def.hpp b/packages/piro/src/Piro_TempusIntegrator_Def.hpp new file mode 100644 index 000000000000..8ddd885e1a7d --- /dev/null +++ b/packages/piro/src/Piro_TempusIntegrator_Def.hpp @@ -0,0 +1,321 @@ +// @HEADER +// ************************************************************************ +// +// Piro: Strategy package for embedded analysis capabilitites +// Copyright (2010) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Andy Salinger (agsalin@sandia.gov), Sandia +// National Laboratories. +// +// ************************************************************************ +// @HEADER + +#include "Piro_TempusIntegrator.hpp" + +#include "Teuchos_ScalarTraits.hpp" +#include "Teuchos_Array.hpp" +#include "Teuchos_Tuple.hpp" +#include "Teuchos_XMLParameterListHelpers.hpp" +#include "Teuchos_Assert.hpp" + +#include +#include +#include + +template +Piro::TempusIntegrator::TempusIntegrator(Teuchos::RCP< Teuchos::ParameterList > pList, + const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, + const SENS_METHOD sens_method) : + out_(Teuchos::VerboseObjectBase::getDefaultOStream()) +{ + if (sens_method == NONE) { + //no sensitivities + basicIntegrator_ = Tempus::integratorBasic(pList, model); + fwdSensIntegrator_ = Teuchos::null; + adjSensIntegrator_ = Teuchos::null; + } + else if (sens_method == FORWARD) { + //forward sensitivities + basicIntegrator_ = Teuchos::null; + fwdSensIntegrator_ = Tempus::integratorForwardSensitivity(pList, model); + adjSensIntegrator_ = Teuchos::null; + } + else if (sens_method == ADJOINT) { + //adjoint sensitivities + basicIntegrator_ = Teuchos::null; + fwdSensIntegrator_ = Teuchos::null; + adjSensIntegrator_ = Tempus::integratorAdjointSensitivity(pList, model); + } +} + +template +Teuchos::RCP> +Piro::TempusIntegrator::getStepper() const +{ + Teuchos::RCP> stepper; + if (basicIntegrator_ != Teuchos::null) { + stepper = basicIntegrator_->getStepper(); + } + if (fwdSensIntegrator_ != Teuchos::null) { + stepper = fwdSensIntegrator_->getStepper(); + } + if (adjSensIntegrator_ != Teuchos::null) { + stepper = adjSensIntegrator_->getStepper(); + } + return stepper; +} + +template +bool +Piro::TempusIntegrator::advanceTime(const Scalar time_final) +{ + bool out; + if (basicIntegrator_ != Teuchos::null) { + out = basicIntegrator_->advanceTime(time_final); + } + if (fwdSensIntegrator_ != Teuchos::null) { + out = fwdSensIntegrator_->advanceTime(time_final); + } + if (adjSensIntegrator_ != Teuchos::null) { + out = adjSensIntegrator_->advanceTime(time_final); + } + return out; +} + +template +Scalar +Piro::TempusIntegrator::getTime() const +{ + Scalar time; + if (basicIntegrator_ != Teuchos::null) { + time = basicIntegrator_->getTime(); + } + if (fwdSensIntegrator_ != Teuchos::null) { + time = fwdSensIntegrator_->getTime(); + } + if (adjSensIntegrator_ != Teuchos::null) { + time = adjSensIntegrator_->getTime(); + } + return time; +} + +template +Teuchos::RCP> +Piro::TempusIntegrator::getX() const +{ + Teuchos::RCP> x; + if (basicIntegrator_ != Teuchos::null) { + x = basicIntegrator_->getX(); + } + if (fwdSensIntegrator_ != Teuchos::null) { + x = fwdSensIntegrator_->getX(); + } + if (adjSensIntegrator_ != Teuchos::null) { + x = adjSensIntegrator_->getX(); + } + return x; +} + +template +Teuchos::RCP> +Piro::TempusIntegrator::getSolutionHistory() const +{ + Teuchos::RCP> soln_history; + if (basicIntegrator_ != Teuchos::null) { + soln_history = basicIntegrator_->getSolutionHistory(); + } + if (fwdSensIntegrator_ != Teuchos::null) { + soln_history = fwdSensIntegrator_->getSolutionHistory(); + } + if (adjSensIntegrator_ != Teuchos::null) { + soln_history = adjSensIntegrator_->getSolutionHistory(); + } + return soln_history; +} + +template +Teuchos::RCP> +Piro::TempusIntegrator::getTimeStepControl() const +{ + Teuchos::RCP> ts_control; + if (basicIntegrator_ != Teuchos::null) { + ts_control = basicIntegrator_->getTimeStepControl(); + } + if (fwdSensIntegrator_ != Teuchos::null) { + ts_control = fwdSensIntegrator_->getTimeStepControl(); + } + if (adjSensIntegrator_ != Teuchos::null) { + ts_control = adjSensIntegrator_->getTimeStepControl(); + } + return ts_control; +} + +template +void +Piro::TempusIntegrator::clearObservers() +{ + if (basicIntegrator_ != Teuchos::null) { + basicIntegrator_->getObserver()->clearObservers(); + } + //IKT: fwdSensIntegrator and adjSensIntegrator do not support composite + //observer so clearObservers() routine is not relevant. If a + //composite observer is ever added to these integrators, we'd need + //to uncomment the code below + /*if (fwdSensIntegrator_ != Teuchos::null) { + fwdSensIntegrator_->getObserver()->clearObservers(); + } + if (adjSensIntegrator_ != Teuchos::null) { + adjSensIntegrator_->getObserver()->clearObservers(); + }*/ +} + +template +void +Piro::TempusIntegrator::setObserver(Teuchos::RCP> obs) +{ + if (basicIntegrator_ != Teuchos::null) { + basicIntegrator_->setObserver(obs); + } + if (fwdSensIntegrator_ != Teuchos::null) { + fwdSensIntegrator_->setObserver(obs); + } + //IKT, FIXME: adjSensIntegrator_ has no routine setObsever(obs) + //Look into. + /*if (adjSensIntegrator_ != Teuchos::null) { + adjSensIntegrator_->setObserver(obs); + }*/ +} + +template +void +Piro::TempusIntegrator::initialize() +{ + if (basicIntegrator_ != Teuchos::null) { + basicIntegrator_->initialize(); + } + if (fwdSensIntegrator_ != Teuchos::null) { + fwdSensIntegrator_->initialize(); + } + //IKT FIXME: adjSensIntegrator_ has no initialize() routine. + //Look into. + /*if (adjSensIntegrator_ != Teuchos::null) { + adjSensIntegrator_->initialize(); + }*/ +} + +template +void +Piro::TempusIntegrator::initializeSolutionHistory(Scalar t0, + Teuchos::RCP< const Thyra::VectorBase< Scalar > > x0, + Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdot0, + Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdotdot0, + Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxDp0, + Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotDp0, + Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotdotDp0) +{ + if (basicIntegrator_ != Teuchos::null) { + basicIntegrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0); + } + if (fwdSensIntegrator_ != Teuchos::null) { + fwdSensIntegrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0, DxDp0, DxdotDp0, DxdotdotDp0); + } + if (adjSensIntegrator_ != Teuchos::null) { + adjSensIntegrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0, DxDp0, DxdotDp0, DxdotdotDp0); + } +} + +template +Tempus::Status +Piro::TempusIntegrator::getStatus() const +{ + Tempus::Status status; + if (basicIntegrator_ != Teuchos::null) { + status = basicIntegrator_->getStatus(); + } + if (fwdSensIntegrator_ != Teuchos::null) { + status = fwdSensIntegrator_->getStatus(); + } + if (adjSensIntegrator_ != Teuchos::null) { + status = adjSensIntegrator_->getStatus(); + } + return status; +} + +template +Teuchos::RCP> +Piro::TempusIntegrator::getDxDp() const +{ + if (fwdSensIntegrator_ != Teuchos::null) { + return fwdSensIntegrator_->getDxDp(); + } + else { + return Teuchos::null; + } +} + +template +Teuchos::RCP> +Piro::TempusIntegrator::getDxdotDp() const +{ + if (fwdSensIntegrator_ != Teuchos::null) { + return fwdSensIntegrator_->getDxdotDp(); + } + else { + return Teuchos::null; + } +} + +template +Teuchos::RCP> +Piro::TempusIntegrator::getDxdotdotDp() const +{ + if (fwdSensIntegrator_ != Teuchos::null) { + return fwdSensIntegrator_->getDxdotdotDp(); + } + else { + return Teuchos::null; + } +} + + +template +Teuchos::RCP> +Piro::TempusIntegrator::getDgDp() const +{ + if (adjSensIntegrator_ != Teuchos::null) { + return adjSensIntegrator_->getDgDp(); + } + else { + return Teuchos::null; + } +} diff --git a/packages/piro/src/Piro_TempusSolver.hpp b/packages/piro/src/Piro_TempusSolver.hpp index cc28c33fea6c..318c287d391a 100644 --- a/packages/piro/src/Piro_TempusSolver.hpp +++ b/packages/piro/src/Piro_TempusSolver.hpp @@ -46,14 +46,12 @@ #include "Piro_ConfigDefs.hpp" #include "Thyra_ResponseOnlyModelEvaluatorBase.hpp" -#include "Tempus_IntegratorBasic.hpp" -#include "Tempus_IntegratorObserverBasic.hpp" - #include "Piro_ObserverBase.hpp" #include "Piro_TempusStepperFactory.hpp" #include "Piro_TempusStepControlFactory.hpp" #include "Piro_TransientSolver.hpp" +#include "Piro_Helpers.hpp" #include #include @@ -78,29 +76,29 @@ class TempusSolver TempusSolver( const Teuchos::RCP &appParams, const Teuchos::RCP > &model, - bool computeSensitivities, const Teuchos::RCP > &piroObserver = Teuchos::null); /** \brief Initialize using prebuilt objects. */ TempusSolver( - const Teuchos::RCP > &stateIntegrator, + const Teuchos::RCP > &stateIntegrator, const Teuchos::RCP > &stateStepper, const Teuchos::RCP > &timeStepSolver, const Teuchos::RCP > &model, Scalar finalTime, - const Teuchos::RCP > &initialConditionModel = Teuchos::null, + const std::string sens_method_string = "None", Teuchos::EVerbosityLevel verbosityLevel = Teuchos::VERB_DEFAULT); - //@} + + //@} /** \brief Initialize using prebuilt objects - supplying initial time value. */ TempusSolver( - const Teuchos::RCP > &stateIntegrator, + const Teuchos::RCP > &stateIntegrator, const Teuchos::RCP > &stateStepper, const Teuchos::RCP > &timeStepSolver, const Teuchos::RCP > &model, Scalar initialTime, Scalar finalTime, - const Teuchos::RCP > &initialConditionModel = Teuchos::null, + const std::string sens_method_string = "None", Teuchos::EVerbosityLevel verbosityLevel = Teuchos::VERB_DEFAULT); //@} @@ -157,6 +155,9 @@ class TempusSolver Tempus::Status getTempusIntegratorStatus() const; + + Teuchos::RCP> + getPiroTempusIntegrator() const {return piroTempusIntegrator_;} private: /** \name Overridden from Thyra::ModelEvaluatorDefaultBase. */ @@ -171,36 +172,33 @@ class TempusSolver /** \brief . */ Teuchos::RCP getValidTempusParameters() const; - Teuchos::RCP > fwdStateIntegrator; - Teuchos::RCP > fwdStateStepper; - Teuchos::RCP > fwdTimeStepSolver; + Teuchos::RCP> piroTempusIntegrator_; + Teuchos::RCP > fwdStateStepper_; + Teuchos::RCP > fwdTimeStepSolver_; - Teuchos::RCP > model; - Teuchos::RCP > thyraModel; - Teuchos::RCP > initialConditionModel; + Teuchos::RCP > model_; + Teuchos::RCP > thyraModel_; - Scalar t_initial; - Scalar t_final; + Scalar t_initial_; + Scalar t_final_; - int num_p; - int num_g; + int num_p_; + int num_g_; - bool computeSensitivities_; - - Teuchos::RCP out; - Teuchos::EVerbosityLevel solnVerbLevel; + Teuchos::RCP out_; + Teuchos::EVerbosityLevel solnVerbLevel_; // used for adding user defined steppers externally, this gives us "the open-close principal" - std::map > > stepperFactories; + std::map > > stepperFactories_; - std::map > > stepControlFactories; + std::map > > stepControlFactories_; - bool isInitialized; + bool isInitialized_; Teuchos::RCP > piroObserver_; bool supports_x_dotdot_; - + //! Set observer void setObserver(); @@ -212,6 +210,7 @@ class TempusSolver //from Albany. bool abort_on_fail_at_min_dt_; + SENS_METHOD sens_method_; }; /** \brief Non-member constructor function */ diff --git a/packages/piro/src/Piro_TempusSolver_Def.hpp b/packages/piro/src/Piro_TempusSolver_Def.hpp index a4585c913aac..1d5a5ae2197e 100644 --- a/packages/piro/src/Piro_TempusSolver_Def.hpp +++ b/packages/piro/src/Piro_TempusSolver_Def.hpp @@ -76,7 +76,6 @@ # include "Thyra_NonlinearSolver_NOX.hpp" #endif -//#define DEBUG_OUTPUT #include #include @@ -86,18 +85,16 @@ template Piro::TempusSolver::TempusSolver( const Teuchos::RCP &appParams, const Teuchos::RCP > &in_model, - bool computeSensitivities, const Teuchos::RCP > &piroObserver): TransientSolver(in_model), - computeSensitivities_(computeSensitivities), - out(Teuchos::VerboseObjectBase::getDefaultOStream()), - isInitialized(false), + out_(Teuchos::VerboseObjectBase::getDefaultOStream()), + isInitialized_(false), piroObserver_(piroObserver), supports_x_dotdot_(false) { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif + std::string sens_method_string = appParams->get("Sensitivity Method","None"); + this->setSensitivityMethod(sens_method_string); + sens_method_ = this->getSensitivityMethod(); std::string jacobianSource = appParams->get("Jacobian Operator", "Have Jacobian"); if (jacobianSource == "Matrix-Free") { Teuchos::RCP > model; @@ -108,30 +105,27 @@ Piro::TempusSolver::TempusSolver( else model = Teuchos::rcp(new Piro::MatrixFreeDecorator(in_model)); initialize(appParams, model); } - else + else { initialize(appParams, in_model); + } } template void Piro::TempusSolver::initialize( const Teuchos::RCP &appParams, - const Teuchos::RCP< Thyra::ModelEvaluator > &in_model) + const Teuchos::RCP< Thyra::ModelEvaluator > &in_model) { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif using Teuchos::ParameterList; using Teuchos::parameterList; using Teuchos::RCP; using Teuchos::rcp; - // set some internals - model = in_model; - num_p = in_model->Np(); - num_g = in_model->Ng(); + model_ = in_model; + num_p_ = in_model->Np(); + num_g_ = in_model->Ng(); // - *out << "\nA) Get the base parameter list ...\n"; + *out_ << "\nA) Get the base parameter list ...\n"; // if (appParams->isSublist("Tempus")) { @@ -139,18 +133,16 @@ void Piro::TempusSolver::initialize( RCP tempusPL = sublist(appParams, "Tempus", true); abort_on_failure_ = tempusPL->get("Abort on Failure", true); - //*out << "tempusPL = " << *tempusPL << "\n"; RCP integratorPL = sublist(tempusPL, "Tempus Integrator", true); - //*out << "integratorPL = " << *integratorPL << "\n"; //IKT, 10/31/16, FIXME: currently there is no Verbosity Sublist in Tempus, but //Curt will add this at some point. When this option is added, set Verbosity //based on that sublist, rather than hard-coding it here. - solnVerbLevel = Teuchos::VERB_DEFAULT; + solnVerbLevel_ = Teuchos::VERB_DEFAULT; RCP timeStepControlPL = Teuchos::null; RCP albTimeStepControlPL = Teuchos::null; if (tempusPL->isSublist("Albany Time Step Control Options")) { - *out << "\n Using 'Albany Time Step Control Options'.\n"; + *out_ << "\n Using 'Albany Time Step Control Options'.\n"; abort_on_fail_at_min_dt_ = true; albTimeStepControlPL = sublist(tempusPL, "Albany Time Step Control Options"); if (integratorPL->isSublist("Time Step Control")) { @@ -164,8 +156,8 @@ void Piro::TempusSolver::initialize( << "Please re-run with 'Abort on Failure = false' or use Tempus Time Step Control.\n"); } abort_on_failure_ = false; - t_initial = albTimeStepControlPL->get("Initial Time", 0.0); - t_final = albTimeStepControlPL->get("Final Time"); + t_initial_ = albTimeStepControlPL->get("Initial Time", 0.0); + t_final_ = albTimeStepControlPL->get("Final Time"); Scalar dt_initial; dt_initial = albTimeStepControlPL->get("Initial Time Step"); Scalar dt_min = albTimeStepControlPL->get("Minimum Time Step", dt_initial); @@ -173,8 +165,8 @@ void Piro::TempusSolver::initialize( Scalar reduc_factor = albTimeStepControlPL->get("Reduction Factor", 1.0); Scalar ampl_factor = albTimeStepControlPL->get("Amplification Factor", 1.0); timeStepControlPL = sublist(integratorPL, "Time Step Control", false); - timeStepControlPL->set("Initial Time", t_initial); - timeStepControlPL->set("Final Time", t_final); + timeStepControlPL->set("Initial Time", t_initial_); + timeStepControlPL->set("Final Time", t_final_); timeStepControlPL->set("Initial Time Step", dt_initial); timeStepControlPL->set("Minimum Time Step", dt_min); timeStepControlPL->set("Maximum Time Step", dt_max); @@ -194,19 +186,19 @@ void Piro::TempusSolver::initialize( basic_vs_PL->set("Maximum Value Monitoring Function", 1.0e20); } else { - *out << "\n Using Tempus 'Time Step Control'.\n"; + *out_ << "\n Using Tempus 'Time Step Control'.\n"; RCP timeStepControlPL = sublist(integratorPL, "Time Step Control", true); - t_initial = timeStepControlPL->get("Initial Time", 0.0); - t_final = timeStepControlPL->get("Final Time", t_initial); + t_initial_ = timeStepControlPL->get("Initial Time", 0.0); + t_final_ = timeStepControlPL->get("Final Time", t_initial_); } - //*out << "tempusPL = " << *tempusPL << "\n"; + //*out_ << "tempusPL = " << *tempusPL << "\n"; RCP stepperPL = sublist(tempusPL, "Tempus Stepper", true); - //*out << "stepperPL = " << *stepperPL << "\n"; + //*out_ << "stepperPL = " << *stepperPL << "\n"; const std::string stepperType = stepperPL->get("Stepper Type", "Backward Euler"); - //*out << "Stepper Type = " << stepperType << "\n"; + //*out_ << "Stepper Type = " << stepperType << "\n"; // - // *out << "\nB) Create the Stratimikos linear solver factory ...\n"; + // *out_ << "\nB) Create the Stratimikos linear solver factory ...\n"; // // This is the linear solve strategy that will be used to solve for the // linear system with the W. @@ -227,7 +219,7 @@ void Piro::TempusSolver::initialize( RCP > lowsFactory = createLinearSolveStrategy(linearSolverBuilder); // - *out << "\nC) Create and initalize the forward model ...\n"; + *out_ << "\nC) Create and initalize the forward model ...\n"; // // C.1) Create the underlying Thyra::ModelEvaluator @@ -250,13 +242,13 @@ void Piro::TempusSolver::initialize( bool invertMassMatrix = tempusPL->get("Invert Mass Matrix", false); if (!invertMassMatrix) { - *out << "\n WARNING in Piro::TempusSolver! You are attempting to run \n" + *out_ << "\n WARNING in Piro::TempusSolver! You are attempting to run \n" << "Explicit Stepper (" << stepperType << ") with 'Invert Mass Matrix' set to 'false'. \n" << "This option should be set to 'true' unless your mass matrix is the identiy.\n"; } else { - Teuchos::RCP > origModel = model; - model = Teuchos::rcp(new Piro::InvertMassMatrixDecorator( + Teuchos::RCP > origModel = model_; + model_ = Teuchos::rcp(new Piro::InvertMassMatrixDecorator( sublist(tempusPL,"Stratimikos", true), origModel, true, tempusPL->get("Lump Mass Matrix", false),false)); } } @@ -266,33 +258,34 @@ void Piro::TempusSolver::initialize( else if (stepperType == "Newmark Explicit a-Form") { bool invertMassMatrix = tempusPL->get("Invert Mass Matrix", false); if (!invertMassMatrix) { - *out << "\n WARNING in Piro::TempusSolver! You are attempting to run \n" + *out_ << "\n WARNING in Piro::TempusSolver! You are attempting to run \n" << "'Newmark Explicit a-Form' Stepper with 'Invert Mass Matrix' set to 'false'. \n" << "This option should be set to 'true' unless your mass matrix is the identiy.\n"; } else { - Teuchos::RCP > origModel = model; - model = Teuchos::rcp(new Piro::InvertMassMatrixDecorator( + Teuchos::RCP > origModel = model_; + model_ = Teuchos::rcp(new Piro::InvertMassMatrixDecorator( sublist(tempusPL,"Stratimikos", true), origModel, true, tempusPL->get("Lump Mass Matrix", false),true)); } } // C.2) Create the Thyra-wrapped ModelEvaluator - thyraModel = rcp(new Thyra::DefaultModelEvaluatorWithSolveFactory(model, lowsFactory)); + thyraModel_ = rcp(new Thyra::DefaultModelEvaluatorWithSolveFactory(model_, lowsFactory)); - const RCP > x_space = thyraModel->get_x_space(); + const RCP > x_space = thyraModel_->get_x_space(); // - *out << "\nD) Create the stepper and integrator for the forward problem ...\n"; + *out_ << "\nD) Create the stepper and integrator for the forward problem ...\n"; - //Create Tempus integrator with observer using tempusPL and model. - fwdStateIntegrator = Tempus::integratorBasic(tempusPL, model); + //Create Tempus integrator with observer using tempusPL, model_ and sensitivity method + piroTempusIntegrator_ = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, model_, sens_method_)); + this->setPiroTempusIntegrator(piroTempusIntegrator_); //Get stepper from integrator - fwdStateStepper = fwdStateIntegrator->getStepper(); + fwdStateStepper_ = piroTempusIntegrator_->getStepper(); //Set observer - supports_x_dotdot_ = model->createInArgs().supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot_dot); + supports_x_dotdot_ = model_->createInArgs().supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot_dot); setObserver(); } @@ -304,80 +297,68 @@ void Piro::TempusSolver::initialize( } - isInitialized = true; + isInitialized_ = true; } template Piro::TempusSolver::TempusSolver( - const Teuchos::RCP > &stateIntegrator, + const Teuchos::RCP > &stateIntegrator, const Teuchos::RCP > &stateStepper, const Teuchos::RCP > &timeStepSolver, const Teuchos::RCP > &underlyingModel, Scalar finalTime, - const Teuchos::RCP > &icModel, + const std::string sens_method_string, Teuchos::EVerbosityLevel verbosityLevel) : - TransientSolver(underlyingModel, icModel), - fwdStateIntegrator(stateIntegrator), - fwdStateStepper(stateStepper), - fwdTimeStepSolver(timeStepSolver), - model(underlyingModel), - initialConditionModel(icModel), - t_initial(0.0), - t_final(finalTime), - num_p(model->Np()), - num_g(model->Ng()), - computeSensitivities_(false), - out(Teuchos::VerboseObjectBase::getDefaultOStream()), - solnVerbLevel(verbosityLevel), - isInitialized(true) + TransientSolver(underlyingModel), + piroTempusIntegrator_(stateIntegrator), + fwdStateStepper_(stateStepper), + fwdTimeStepSolver_(timeStepSolver), + model_(underlyingModel), + t_initial_(0.0), + t_final_(finalTime), + num_p_(model_->Np()), + num_g_(model_->Ng()), + out_(Teuchos::VerboseObjectBase::getDefaultOStream()), + solnVerbLevel_(verbosityLevel), + isInitialized_(true) { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif - if (fwdStateStepper->getModel() != underlyingModel) { - fwdStateStepper->setModel(underlyingModel); + if (fwdStateStepper_->getModel() != underlyingModel) { + fwdStateStepper_->setModel(underlyingModel); } + this->setSensitivityMethod(sens_method_string); + sens_method_ = this->getSensitivityMethod(); + this->setPiroTempusIntegrator(piroTempusIntegrator_); } template Piro::TempusSolver::TempusSolver( - const Teuchos::RCP > &stateIntegrator, + const Teuchos::RCP > &stateIntegrator, const Teuchos::RCP > &stateStepper, const Teuchos::RCP > &timeStepSolver, const Teuchos::RCP > &underlyingModel, Scalar initialTime, Scalar finalTime, - const Teuchos::RCP > &icModel, + const std::string sens_method_string, Teuchos::EVerbosityLevel verbosityLevel) : - TransientSolver(underlyingModel, icModel), - fwdStateIntegrator(stateIntegrator), - fwdStateStepper(stateStepper), - fwdTimeStepSolver(timeStepSolver), - model(underlyingModel), - initialConditionModel(icModel), - t_initial(initialTime), - t_final(finalTime), - num_p(model->Np()), - num_g(model->Ng()), - computeSensitivities_(false), - out(Teuchos::VerboseObjectBase::getDefaultOStream()), - solnVerbLevel(verbosityLevel), - isInitialized(true) + TransientSolver(underlyingModel), + piroTempusIntegrator_(stateIntegrator), + fwdStateStepper_(stateStepper), + fwdTimeStepSolver_(timeStepSolver), + model_(underlyingModel), + t_initial_(initialTime), + t_final_(finalTime), + num_p_(model_->Np()), + num_g_(model_->Ng()), + out_(Teuchos::VerboseObjectBase::getDefaultOStream()), + solnVerbLevel_(verbosityLevel), + isInitialized_(true) { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif - //IKT, 12/5/16: the following exception check is needed until setInitialCondition method - //is added to the Tempus::IntegratorBasic class. - if (initialTime > 0.0) { - TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, - "\n Error in Piro::TempusSolver: the constructor employed does not support initialTime > 0.0. " << - "You have set initialTime = " << initialTime << "\n"); - } - - if (fwdStateStepper->getModel() != underlyingModel) { - fwdStateStepper->setModel(underlyingModel); + if (fwdStateStepper_->getModel() != underlyingModel) { + fwdStateStepper_->setModel(underlyingModel); } + this->setSensitivityMethod(sens_method_string); + sens_method_ = this->getSensitivityMethod(); + this->setPiroTempusIntegrator(piroTempusIntegrator_); } template @@ -385,199 +366,94 @@ void Piro::TempusSolver::evalModelImpl( const Thyra::ModelEvaluatorBase::InArgs& inArgs, const Thyra::ModelEvaluatorBase::OutArgs& outArgs) const { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif using Teuchos::RCP; using Teuchos::rcp; - // TODO: Support more than 1 parameter and 1 response - const int j = 0; - const int l = 0; - - // Parse InArgs - RCP > p_in; - if (num_p > 0) { - p_in = inArgs.get_p(l); + // Set initial time and initial condition + Thyra::ModelEvaluatorBase::InArgs state_ic = model_->getNominalValues(); + if(t_initial_ > 0.0 && state_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_t)) { + state_ic.set_t(t_initial_); } - RCP > p_in2; //JF add for multipoint - if (num_p > 1) { - p_in2 = inArgs.get_p(l+1); - } - - // Parse OutArgs - RCP > g_out; - if (num_g > 0) { - g_out = outArgs.get_g(j); - } - const RCP > gx_out = outArgs.get_g(num_g); - - Thyra::ModelEvaluatorBase::InArgs state_ic = model->getNominalValues(); - - // Set initial time in ME if needed - - if(t_initial > 0.0 && state_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_t)) { - state_ic.set_t(t_initial); - } - - if (Teuchos::nonnull(initialConditionModel)) { - // The initial condition depends on the parameter - // It is found by querying the auxiliary model evaluator as the last response - const RCP > initialState = - Thyra::createMember(model->get_x_space()); - - { - Thyra::ModelEvaluatorBase::InArgs initCondInArgs = initialConditionModel->createInArgs(); - if (num_p > 0) { - initCondInArgs.set_p(l, inArgs.get_p(l)); - } - - Thyra::ModelEvaluatorBase::OutArgs initCondOutArgs = initialConditionModel->createOutArgs(); - initCondOutArgs.set_g(initCondOutArgs.Ng() - 1, initialState); - - initialConditionModel->evalModel(initCondInArgs, initCondOutArgs); - } - - state_ic.set_x(initialState); - } - - // Set paramters p_in as part of initial conditions - if (num_p > 0) { + + // Set parameters as part of initial conditions + for (int l = 0; l < num_p_; ++l) { + auto p_in = inArgs.get_p(l); if (Teuchos::nonnull(p_in)) { state_ic.set_p(l, p_in); } } - if (num_p > 1) { //JF added for multipoint - if (Teuchos::nonnull(p_in2)) { - state_ic.set_p(l+1, p_in2); - } - } - //*out << "\nstate_ic:\n" << Teuchos::describe(state_ic, solnVerbLevel); - - //JF may need a version of the following for multipoint, i.e. num_p>1, l+1, if we want sensitivities - RCP > dgxdp_out; - Thyra::ModelEvaluatorBase::Derivative dgdp_deriv_out; - if (num_p > 0) { - const Thyra::ModelEvaluatorBase::DerivativeSupport dgxdp_support = - outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, num_g, l); - if (dgxdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) { - const Thyra::ModelEvaluatorBase::Derivative dgxdp_deriv = - outArgs.get_DgDp(num_g, l); - dgxdp_out = dgxdp_deriv.getMultiVector(); - } - - if (num_g > 0) { - const Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support = - outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l); - if (!dgdp_support.none()) { - dgdp_deriv_out = outArgs.get_DgDp(j, l); - } - } - } - - bool requestedSensitivities = true; - if (computeSensitivities_ == true) - requestedSensitivities = Teuchos::nonnull(dgxdp_out) || !dgdp_deriv_out.isEmpty(); - else - requestedSensitivities = false; + //*out_ << "\nstate_ic:\n" << Teuchos::describe(state_ic, solnVerbLevel_); + // + *out_ << "\nE) Solve the forward problem ...\n"; + // RCP > finalSolution; RCP > solutionState; RCP > solutionHistory; - if (!requestedSensitivities) - { - // - *out << "\nE) Solve the forward problem ...\n"; - // - // - *out << "T final requested: " << t_final << " \n"; - - fwdStateIntegrator->advanceTime(t_final); - - double time = fwdStateIntegrator->getTime(); - - *out << "T final actual: " << time << "\n"; - - if (abs(time-t_final) > 1.0e-10) { - if (abort_on_failure_ == true) { - TEUCHOS_TEST_FOR_EXCEPTION( - true, - Teuchos::Exceptions::InvalidParameter, - "\n Error! Piro::TempusSolver: time-integrator did not make it to final time " << - "specified in Input File. Final time in input file is " << t_final << - ", whereas actual final time is " << time << ". If you'd like to " << - "suppress this exception, run with 'Abort on Failure' set to 'false' in " << - "Tempus sublist.\n" ); - } - else { - *out << "\n WARNING: Piro::TempusSolver did not make it to final time, but " - << "solver will not abort since you have specified 'Abort on Failure' = 'false'.\n"; - } + + *out_ << "T final requested: " << t_final_ << " \n"; + piroTempusIntegrator_->advanceTime(t_final_); + double time = piroTempusIntegrator_->getTime(); + *out_ << "T final actual: " << time << "\n"; + + Scalar diff = 0.0; + if (abs(t_final_) == 0) diff = abs(time-t_final_); + else diff = abs(time-t_final_)/abs(t_final_); + if (diff > 1.0e-10) { + if (abort_on_failure_ == true) { + TEUCHOS_TEST_FOR_EXCEPTION( + true, + Teuchos::Exceptions::InvalidParameter, + "\n Error! Piro::TempusSolver: time-integrator did not make it to final time " << + "specified in Input File. Final time in input file is " << t_final_ << + ", whereas actual final time is " << time << ". If you'd like to " << + "suppress this exception, run with 'Abort on Failure' set to 'false' in " << + "Tempus sublist.\n" ); } - - finalSolution = fwdStateIntegrator->getX(); - - solutionHistory = fwdStateIntegrator->getSolutionHistory(); - auto numStates = solutionHistory->getNumStates(); - solutionState = (*solutionHistory)[numStates-1]; - //Get final solution from solutionHistory. - finalSolution = solutionState->getX(); - - if (Teuchos::VERB_MEDIUM <= solnVerbLevel) { - *out << "Final Solution\n" << *finalSolution << "\n"; + else { + *out_ << "\n WARNING: Piro::TempusSolver did not make it to final time, but " + << "solver will not abort since you have specified 'Abort on Failure' = 'false'.\n"; } - } - else { - // - *out << "\nE) Solve the forward problem with Sensitivities...\n"; - // - TEUCHOS_TEST_FOR_EXCEPTION( - true, - Teuchos::Exceptions::InvalidParameter, - "\n Error! Piro::TempusSolver: sensitivities with Tempus are not yet supported!"); + + solutionHistory = piroTempusIntegrator_->getSolutionHistory(); + auto numStates = solutionHistory->getNumStates(); + solutionState = (*solutionHistory)[numStates-1]; + //Get final solution from solutionHistory. + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + Teuchos::RCP> x = solutionState->getX(); + Teuchos::RCP X = Teuchos::rcp_dynamic_cast(x); + finalSolution = (sens_method_ == NONE) ? x : X->getMultiVector()->col(0); + + if (Teuchos::VERB_MEDIUM <= solnVerbLevel_) { + *out_ << "Final Solution\n" << *finalSolution << "\n"; } - *out << "\nF) Check the solution to the forward problem ...\n"; // As post-processing step, calculate responses at final solution - { - Thyra::ModelEvaluatorBase::InArgs modelInArgs = model->createInArgs(); - { - modelInArgs.set_x(finalSolution); - if (num_p > 0) { - modelInArgs.set_p(l, p_in); - } - if (num_p > 1) { //JF added for multipoint - modelInArgs.set_p(l+1, p_in2); - } - //Set time to be final time at which the solve occurs (< t_final in the case we don't make it to t_final). - //IKT: get final time from solutionHistory workingSpace, which is different than how it is done in Piro::RythmosSolver class. - //IKT, 11/1/16, FIXME? workingState pointer is null right now, so the following - //code is commented out for now. Use t_final and soln_dt in set_t instead for now. - /*RCP > workingState = solutionHistory->getWorkingState(); - const Scalar time = workingState->getTime(); - const Scalar dt = workingState->getTimeStep(); - const Scalar t = time + dt; - modelInArgs.set_t(t);*/ - const Scalar soln_dt = solutionState->getTimeStep(); - modelInArgs.set_t(t_final - soln_dt); - } - - Thyra::ModelEvaluatorBase::OutArgs modelOutArgs = model->createOutArgs(); - if (Teuchos::nonnull(g_out)) { - Thyra::put_scalar(Teuchos::ScalarTraits::zero(), g_out.ptr()); - modelOutArgs.set_g(j, g_out); - } - - model->evalModel(modelInArgs, modelOutArgs); + Thyra::ModelEvaluatorBase::InArgs modelInArgs = model_->createInArgs(); + + modelInArgs.set_x(finalSolution); + for (int l=0; l < num_p_; ++l) { + auto p_in = inArgs.get_p(l); + modelInArgs.set_p(l, p_in); } + //Set time to be final time at which the solve occurs (< t_final_ in the case we don't make it to t_final_). + //IKT: get final time from solutionHistory workingSpace, which is different than how it is done in Piro::RythmosSolver class. + //IKT, 11/1/16, FIXME? workingState pointer is null right now, so the following + //code is commented out for now. Use t_final_ and soln_dt in set_t instead for now. + /*RCP > workingState = solutionHistory->getWorkingState(); + const Scalar time = workingState->getTime(); + const Scalar dt = workingState->getTimeStep(); + const Scalar t = time + dt; + modelInArgs.set_t(t);*/ + const Scalar soln_dt = solutionState->getTimeStep(); + modelInArgs.set_t(t_final_ - soln_dt); + + //Calculate responses and sensitivities + this->evalConvergedModel(modelInArgs, outArgs); - // Return the final solution as an additional g-vector, if requested - if (Teuchos::nonnull(gx_out)) { - Thyra::copy(*finalSolution, gx_out.ptr()); - } } @@ -613,7 +489,7 @@ template void Piro::TempusSolver:: addStepperFactory(const std::string & stepperName,const Teuchos::RCP > & factory) { - stepperFactories[stepperName] = factory; + stepperFactories_[stepperName] = factory; } template @@ -621,14 +497,14 @@ void Piro::TempusSolver:: addStepControlFactory(const std::string & stepControlName, const Teuchos::RCP> & step_control_strategy) { - stepControlFactories[stepControlName] = step_control_strategy; + stepControlFactories_[stepControlName] = step_control_strategy; } template void Piro::TempusSolver:: setStartTime(const Scalar start_time) { - Teuchos::RCP > tsc_const = fwdStateIntegrator->getTimeStepControl(); + Teuchos::RCP > tsc_const = piroTempusIntegrator_->getTimeStepControl(); Teuchos::RCP > tsc = Teuchos::rcp_const_cast >(tsc_const); tsc->setInitTime(start_time); } @@ -637,7 +513,7 @@ template Scalar Piro::TempusSolver:: getStartTime() const { - Teuchos::RCP > tsc = fwdStateIntegrator->getTimeStepControl(); + Teuchos::RCP > tsc = piroTempusIntegrator_->getTimeStepControl(); Scalar start_time = tsc->getInitTime(); return start_time; } @@ -646,9 +522,9 @@ template void Piro::TempusSolver:: setFinalTime(const Scalar final_time) { - Teuchos::RCP > tsc_const = fwdStateIntegrator->getTimeStepControl(); + Teuchos::RCP > tsc_const = piroTempusIntegrator_->getTimeStepControl(); Teuchos::RCP > tsc = Teuchos::rcp_const_cast >(tsc_const); - t_final = final_time; + t_final_ = final_time; tsc->setFinalTime(final_time); } @@ -656,7 +532,7 @@ template Scalar Piro::TempusSolver:: getFinalTime() const { - Teuchos::RCP > tsc = fwdStateIntegrator->getTimeStepControl(); + Teuchos::RCP > tsc = piroTempusIntegrator_->getTimeStepControl(); Scalar final_time = tsc->getFinalTime(); return final_time; } @@ -665,7 +541,7 @@ template void Piro::TempusSolver:: setInitTimeStep(const Scalar init_time_step) { - Teuchos::RCP > tsc_const = fwdStateIntegrator->getTimeStepControl(); + Teuchos::RCP > tsc_const = piroTempusIntegrator_->getTimeStepControl(); Teuchos::RCP > tsc = Teuchos::rcp_const_cast >(tsc_const); tsc->setInitTimeStep(init_time_step); } @@ -675,7 +551,7 @@ template Scalar Piro::TempusSolver:: getInitTimeStep() const { - Teuchos::RCP > tsc = fwdStateIntegrator->getTimeStepControl(); + Teuchos::RCP > tsc = piroTempusIntegrator_->getTimeStepControl(); auto init_time_step = tsc->getInitTimeStep(); return init_time_step; } @@ -686,19 +562,19 @@ setObserver() Teuchos::RCP > observer = Teuchos::null; if (Teuchos::nonnull(piroObserver_)) { //Get solutionHistory from integrator - const Teuchos::RCP > solutionHistory = fwdStateIntegrator->getSolutionHistory(); - const Teuchos::RCP > timeStepControl = fwdStateIntegrator->getTimeStepControl(); + const Teuchos::RCP > solutionHistory = piroTempusIntegrator_->getSolutionHistory(); + const Teuchos::RCP > timeStepControl = piroTempusIntegrator_->getTimeStepControl(); //Create Tempus::IntegratorObserverBasic object observer = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter(solutionHistory, - timeStepControl, piroObserver_, supports_x_dotdot_, abort_on_fail_at_min_dt_)); + timeStepControl, piroObserver_, supports_x_dotdot_, abort_on_fail_at_min_dt_, sens_method_)); } if (Teuchos::nonnull(observer)) { //Set observer in integrator - fwdStateIntegrator->getObserver()->clearObservers(); - fwdStateIntegrator->setObserver(observer); - fwdStateStepper->initialize(); + piroTempusIntegrator_->clearObservers(); + piroTempusIntegrator_->setObserver(observer); + fwdStateStepper_->initialize(); //Reinitialize everything in integrator class, since we have changed the observer. - fwdStateIntegrator->initialize(); + piroTempusIntegrator_->initialize(); } } @@ -709,7 +585,7 @@ setInitialState(Scalar t0, Teuchos::RCP > xdot0, Teuchos::RCP > xdotdot0) { - fwdStateIntegrator->initializeSolutionHistory(t0, x0, xdot0, xdotdot0); + piroTempusIntegrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0); //Reset observer. This is necessary for correct observation of solution //since initializeSolutionHistory modifies the solutionHistory object. setObserver(); @@ -720,15 +596,15 @@ template void Piro::TempusSolver:: setInitialGuess(Teuchos::RCP< const Thyra::VectorBase > initial_guess) { - fwdStateStepper->setInitialGuess(initial_guess); - fwdStateStepper->initialize(); + fwdStateStepper_->setInitialGuess(initial_guess); + fwdStateStepper_->initialize(); } template Teuchos::RCP > Piro::TempusSolver:: getSolutionHistory() const { - Teuchos::RCP > soln_history_const = fwdStateIntegrator->getSolutionHistory(); + Teuchos::RCP > soln_history_const = piroTempusIntegrator_->getSolutionHistory(); Teuchos::RCP > soln_history = Teuchos::rcp_const_cast >(soln_history_const); return soln_history; } @@ -738,7 +614,7 @@ template Teuchos::RCP > Piro::TempusSolver:: getSolver() const { - return fwdStateStepper->getSolver(); + return fwdStateStepper_->getSolver(); } @@ -746,7 +622,7 @@ template Tempus::Status Piro::TempusSolver:: getTempusIntegratorStatus() const { - return fwdStateIntegrator->getStatus(); + return piroTempusIntegrator_->getStatus(); } @@ -757,17 +633,8 @@ Piro::tempusSolver( const Teuchos::RCP > &in_model, const Teuchos::RCP > &piroObserver) { - Teuchos::RCP out(Teuchos::VerboseObjectBase::getDefaultOStream()); -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif - bool computeSensitivities = true; - if (appParams->isSublist("Analysis")) { - Teuchos::ParameterList& analysisPL = appParams->sublist("Analysis"); - if (analysisPL.isParameter("Compute Sensitivities")) - computeSensitivities = analysisPL.get("Compute Sensitivities"); - } + Teuchos::RCP out_(Teuchos::VerboseObjectBase::getDefaultOStream()); + return Teuchos::rcp(new TempusSolver(appParams, in_model, piroObserver)); +} - return Teuchos::rcp(new TempusSolver(appParams, in_model, computeSensitivities, piroObserver)); -} diff --git a/packages/piro/src/Piro_TransientSolver.hpp b/packages/piro/src/Piro_TransientSolver.hpp index 26eae8892112..d61921829c11 100644 --- a/packages/piro/src/Piro_TransientSolver.hpp +++ b/packages/piro/src/Piro_TransientSolver.hpp @@ -45,6 +45,8 @@ #include "Piro_ConfigDefs.hpp" #include "Thyra_ResponseOnlyModelEvaluatorBase.hpp" +#include "Piro_TempusIntegrator.hpp" +#include "Piro_Helpers.hpp" #include #include @@ -61,14 +63,12 @@ class TransientSolver /** \name Constructors/initializers */ //@{ /** \brief . */ - explicit TransientSolver(const Teuchos::RCP >&model, - const Teuchos::RCP > &initialConditionModel = Teuchos::null); + explicit TransientSolver(const Teuchos::RCP >&model); /** \brief . */ TransientSolver( const Teuchos::RCP > &model, - int numParameters, - const Teuchos::RCP > &initialConditionModel = Teuchos::null); + int numParameters); //@} /** \name Overridden from Thyra::ModelEvaluatorBase. */ @@ -83,6 +83,30 @@ class TransientSolver Teuchos::RCP > get_g_space(int j) const; //@} + /** \name Getters for subbclasses */ + //@{ + /** \brief . */ + const Thyra::ModelEvaluator &getModel() const; + + /** \brief . */ + int num_p() const; + + /** \brief . */ + int num_g() const; + + /** \brief . */ + SENS_METHOD getSensitivityMethod(); + //@} + + /** \name Setters for subbclasses */ + /** \brief . */ + void setSensitivityMethod(const std::string sensitivity_method_string); + //@} + + /** \brief . */ + void setPiroTempusIntegrator(Teuchos::RCP> piroTempusIntegrator); + //@} + protected: /** \name Service methods for subclasses. */ //@{ @@ -102,13 +126,14 @@ class TransientSolver Teuchos::RCP > create_DgDp_op_impl(int j, int l) const; - Teuchos::RCP out; + Teuchos::RCP out_; Teuchos::RCP > model_; - Teuchos::RCP > initialConditionModel_; + Teuchos::RCP> piroTempusIntegrator_; int num_p_; int num_g_; + SENS_METHOD sensitivityMethod_; }; diff --git a/packages/piro/src/Piro_TransientSolver_Def.hpp b/packages/piro/src/Piro_TransientSolver_Def.hpp index dffa4162e334..563795505179 100644 --- a/packages/piro/src/Piro_TransientSolver_Def.hpp +++ b/packages/piro/src/Piro_TransientSolver_Def.hpp @@ -54,53 +54,41 @@ #include "Thyra_VectorStdOps.hpp" #include "Thyra_DefaultModelEvaluatorWithSolveFactory.hpp" -#define DEBUG_OUTPUT - #include #include #include template Piro::TransientSolver::TransientSolver( - const Teuchos::RCP > &model, - const Teuchos::RCP > &icModel): - out(Teuchos::VerboseObjectBase::getDefaultOStream()), + const Teuchos::RCP > &model) : + out_(Teuchos::VerboseObjectBase::getDefaultOStream()), model_(model), - initialConditionModel_(icModel), num_p_(model->Np()), num_g_(model->Ng()) { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif + //Nothing to do } template Piro::TransientSolver::TransientSolver( - const Teuchos::RCP > &model, int numParameters, - const Teuchos::RCP > &icModel) : - out(Teuchos::VerboseObjectBase::getDefaultOStream()), + const Teuchos::RCP > &model, int numParameters) : + out_(Teuchos::VerboseObjectBase::getDefaultOStream()), model_(model), - initialConditionModel_(icModel), num_p_(numParameters), - num_g_(model->Ng()) + num_g_(model->Ng()), + sensitivityMethod_(NONE) { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif + //Nothing to do } template Teuchos::RCP > Piro::TransientSolver::get_p_space(int l) const { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif TEUCHOS_TEST_FOR_EXCEPTION( l >= num_p_ || l < 0, Teuchos::Exceptions::InvalidParameter, - "\n Error in Piro::TransientSolver::get_p_map(): " << + "\n Error in Piro::TransientSolver::get_p_space(): " << "Invalid parameter index l = " << l << "\n"); @@ -111,13 +99,10 @@ template Teuchos::RCP > Piro::TransientSolver::get_g_space(int j) const { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif TEUCHOS_TEST_FOR_EXCEPTION( j > num_g_ || j < 0, Teuchos::Exceptions::InvalidParameter, - "\n Error in Piro::TransientSolver::get_g_map(): " << + "\n Error in Piro::TransientSolver::get_g_space(): " << "Invalid response index j = " << j << "\n"); @@ -133,9 +118,6 @@ template Thyra::ModelEvaluatorBase::InArgs Piro::TransientSolver::getNominalValues() const { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif Thyra::ModelEvaluatorBase::InArgs result = this->createInArgs(); const Thyra::ModelEvaluatorBase::InArgs modelNominalValues = model_->getNominalValues(); for (int l = 0; l < num_p_; ++l) { @@ -148,9 +130,6 @@ template Thyra::ModelEvaluatorBase::InArgs Piro::TransientSolver::createInArgs() const { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif Thyra::ModelEvaluatorBase::InArgsSetup inArgs; inArgs.setModelEvalDescription(this->description()); inArgs.set_Np(num_p_); @@ -161,9 +140,6 @@ template Thyra::ModelEvaluatorBase::OutArgs Piro::TransientSolver::createOutArgsImpl() const { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif Thyra::ModelEvaluatorBase::OutArgsSetup outArgs; outArgs.setModelEvalDescription(this->description()); @@ -172,59 +148,52 @@ Piro::TransientSolver::createOutArgsImpl() const const Thyra::ModelEvaluatorBase::OutArgs modelOutArgs = model_->createOutArgs(); - if (num_p_ > 0) { - // Only one parameter supported - const int l = 0; - - if (Teuchos::nonnull(initialConditionModel_)) { - const Thyra::ModelEvaluatorBase::OutArgs initCondOutArgs = - initialConditionModel_->createOutArgs(); - const Thyra::ModelEvaluatorBase::DerivativeSupport init_dxdp_support = - initCondOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, initCondOutArgs.Ng() - 1, l); - if (!init_dxdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) { - // Ok to return early since only one parameter supported + outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f, modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_f)); + + // Sensitivity support (Forward approach only, for now) + // Jacobian solver required for all sensitivities + if (modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_W)) { + for (int l = 0; l < num_p_; ++l) { + // Solution sensitivities: DxDp(l) + // DfDp(l) required + const Thyra::ModelEvaluatorBase::DerivativeSupport dfdp_support = + modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, l); + if (!dfdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) { return outArgs; } - } - - // Computing the DxDp sensitivity for a transient problem currently requires the evaluation of - // the mutilivector-based, Jacobian-oriented DfDp derivatives of the underlying transient model. - const Thyra::ModelEvaluatorBase::DerivativeSupport model_dfdp_support = - modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, l); - if (!model_dfdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) { - // Ok to return early since only one parameter supported - return outArgs; - } + const bool dxdp_linOpSupport = + dfdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP); + const bool dxdp_mvJacSupport = + dfdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM); + Thyra::ModelEvaluatorBase::DerivativeSupport dxdp_support; + if (dxdp_linOpSupport) { + dxdp_support.plus(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP); + } + if (dxdp_mvJacSupport) { + dxdp_support.plus(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM); + } + outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, num_g_, l, dxdp_support); - // Solution sensitivity - outArgs.setSupports( - Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, - num_g_, - l, - Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM); - - if (num_g_ > 0) { - // Only one response supported - const int j = 0; - - const Thyra::ModelEvaluatorBase::DerivativeSupport model_dgdx_support = - modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDx, j); - if (!model_dgdx_support.none()) { - const Thyra::ModelEvaluatorBase::DerivativeSupport model_dgdp_support = - modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l); - // Response sensitivity - Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support; - if (model_dgdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) { - dgdp_support.plus(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM); - } - if (model_dgdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP)) { - dgdp_support.plus(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP); + // Response sensitivities: DgDp(j, l) + // DxDp(l) required + if (dxdp_linOpSupport || dxdp_mvJacSupport) { + for (int j = 0; j < num_g_; ++j) { + const Thyra::ModelEvaluatorBase::DerivativeSupport model_dgdx_support = + modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDx, j); + if (!model_dgdx_support.none()) { + const Thyra::ModelEvaluatorBase::DerivativeSupport model_dgdp_support = + modelOutArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l); + // Response sensitivity + Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support; + if (model_dgdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) { + dgdp_support.plus(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM); + } + if (model_dgdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP)) { + dgdp_support.plus(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP); + } + outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l, dgdp_support); + } } - outArgs.setSupports( - Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, - j, - l, - dgdp_support); } } } @@ -237,20 +206,274 @@ template Teuchos::RCP > Piro::TransientSolver::create_DgDp_op_impl(int j, int l) const { - TEUCHOS_ASSERT(j != num_g_); + TEUCHOS_TEST_FOR_EXCEPTION( + j > num_g_ || j < 0, + Teuchos::Exceptions::InvalidParameter, + "\n Error in Piro::TransientSolver::create_DgDp_op_impl(): " << + "Invalid response index j = " << + j << "\n"); + TEUCHOS_TEST_FOR_EXCEPTION( + l >= num_p_ || l < 0, + Teuchos::Exceptions::InvalidParameter, + "\n Error in Piro::TransientSolver::create_DgDp_op_impl(): " << + "Invalid parameter index l = " << + l << "\n"); const Teuchos::Array > > dummy = Teuchos::tuple(Thyra::zero(this->get_g_space(j), this->get_p_space(l))); - return Teuchos::rcp(new Thyra::DefaultAddedLinearOp(dummy)); + if (j == num_g_) { + return Thyra::defaultMultipliedLinearOp(dummy); + } else { + return Teuchos::rcp(new Thyra::DefaultAddedLinearOp(dummy)); + } +} + + +template +const Thyra::ModelEvaluator & +Piro::TransientSolver::getModel() const +{ + return model_; +} + +template +int +Piro::TransientSolver::num_p() const +{ + return num_p_; +} + +template +int +Piro::TransientSolver::num_g() const +{ + return num_g_; } +template +void +Piro::TransientSolver::setSensitivityMethod(const std::string sensitivity_method_string) +{ + if (sensitivity_method_string == "None") sensitivityMethod_ = NONE; + else if (sensitivity_method_string == "Forward") sensitivityMethod_ = FORWARD; + else if (sensitivity_method_string == "Adjoint") sensitivityMethod_ = ADJOINT; + else { + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + "\n Error! Piro::TransientSolver: invalid Sensitivity Method = " << sensitivity_method_string << "! \n" + << " Valid options for Sensitivity Method are 'None', 'Forward' and 'Adjoint'.\n"); + } + //IKT, 5/8/2020: remove the following once we have support for adjoint sensitivities + if (sensitivityMethod_ == ADJOINT) { + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + "\n Error! Piro::TransientSolver: adjoint sentivities (Sensitivity Method = " + << "Adjoint) are not yet supported! Please set 'Sensitivity Method' to 'None' or 'Forward'.\n"); + } +} + +template +Piro::SENS_METHOD +Piro::TransientSolver::getSensitivityMethod() +{ + return sensitivityMethod_; +} + +template +void +Piro::TransientSolver::setPiroTempusIntegrator(Teuchos::RCP> piroTempusIntegrator) +{ + piroTempusIntegrator_ = piroTempusIntegrator; +} + + template void Piro::TransientSolver::evalConvergedModel( const Thyra::ModelEvaluatorBase::InArgs& modelInArgs, const Thyra::ModelEvaluatorBase::OutArgs& outArgs) const { -#ifdef DEBUG_OUTPUT - *out << "DEBUG: " << __PRETTY_FUNCTION__ << "\n"; -#endif - //IKT FIXME: FILL IN! + using Teuchos::RCP; + using Teuchos::rcp; + + *out_ << "\nE) Calculate responses ...\n"; + + // Solution at convergence is the response at index num_g_ + RCP > gx_out = outArgs.get_g(num_g_); + if (Teuchos::nonnull(gx_out)) { + Thyra::copy(*modelInArgs.get_x(), gx_out.ptr()); + } + + // Setup output for final evalution of underlying model + Thyra::ModelEvaluatorBase::OutArgs modelOutArgs = model_->createOutArgs(); + + //Responses + for (int j=0; j::zero(), g_out.ptr()); + modelOutArgs.set_g(j, g_out); + } + } + + // DgDx derivatives + for (int j = 0; j < num_g_; ++j) { + Thyra::ModelEvaluatorBase::DerivativeSupport dgdx_request; + for (int l = 0; l < num_p_; ++l) { + const Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support = + outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l); + if (!dgdp_support.none()) { + const Thyra::ModelEvaluatorBase::Derivative dgdp_deriv = + outArgs.get_DgDp(j, l); + if (!dgdp_deriv.isEmpty()) { + const bool dgdp_mvGrad_required = + Teuchos::nonnull(dgdp_deriv.getMultiVector()) && + dgdp_deriv.getMultiVectorOrientation() == Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM; + if (dgdp_mvGrad_required) { + dgdx_request.plus(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM); + } + else { + dgdx_request.plus(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP); + } + } + } + } + + if (!dgdx_request.none()) { + Thyra::ModelEvaluatorBase::Derivative dgdx_deriv; + if (dgdx_request.supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM)) { + dgdx_deriv = Thyra::create_DgDx_mv(*model_, j, Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM); + } + else if (dgdx_request.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP)) { + dgdx_deriv = model_->create_DgDx_op(j); + } + modelOutArgs.set_DgDx(j, dgdx_deriv); + } + } + + // DgDp derivatives + for (int l = 0; l < num_p_; ++l) { + for (int j = 0; j < num_g_; ++j) { + const Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support = + outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l); + if (!dgdp_support.none()) { + const Thyra::ModelEvaluatorBase::Derivative dgdp_deriv = + outArgs.get_DgDp(j, l); + Thyra::ModelEvaluatorBase::Derivative model_dgdp_deriv; + const RCP > dgdp_op = dgdp_deriv.getLinearOp(); + if (Teuchos::nonnull(dgdp_op)) { + model_dgdp_deriv = model_->create_DgDp_op(j, l); + } + else { + model_dgdp_deriv = dgdp_deriv; + } + if (!model_dgdp_deriv.isEmpty()) { + modelOutArgs.set_DgDp(j, l, model_dgdp_deriv); + } + } + } + } + + model_->evalModel(modelInArgs, modelOutArgs); + + // Check if sensitivities are requested + bool requestedSensitivities = false; + for (int i=0; i > dxdp_mv = piroTempusIntegrator_->getDxDp(); + //IMPORTANT REMARK: we are currently NOT using DxdotDp and DxdotdotDp in transient sensitivities! + //The capability to use them can be added at a later point in time, if desired. + //IKT, 5/10/20: throw error if dxdp_mv returned by Tempus is null. Not sure if this can happen in practice or not... + if (dxdp_mv == Teuchos::null) { + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + "\n Error! Piro::TransientSolver: DxDp returned by Tempus::IntegratorForwardSensitivity::getDxDp() routine is null!\n"); + } + //IKT FIXME: probably a lot of this can be reused for adjoint sensitivities. Will clean up later, + //when adjoint sensitivities are implemented. + for (int l = 0; l < num_p_; ++l) { + for (int j = 0; j < num_g_; ++j) { + //Get DgDp and DgDx + const Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support = + outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l); + if (!dgdp_support.none()) { + const Thyra::ModelEvaluatorBase::Derivative dgdp_deriv = outArgs.get_DgDp(j, l); + if (!dgdp_deriv.isEmpty()) { + const Thyra::ModelEvaluatorBase::Derivative dgdx_deriv = modelOutArgs.get_DgDx(j); + const RCP > dgdx_mv = dgdx_deriv.getMultiVector(); + RCP > dgdx_op = dgdx_deriv.getLinearOp(); + if (Teuchos::is_null(dgdx_op)) { + //NOTE: dgdx_mv is the transpose, so by calling Thyra::adjoint on dgdx_mv, + //we get the untransposed operator back as dgdx_op + dgdx_op = Thyra::adjoint(dgdx_mv); + } + const RCP > dgdp_op = dgdp_deriv.getLinearOp(); + if (Teuchos::nonnull(dgdp_op)) { + //Case 1: DgDp, DgDx and DxDp are linear ops. This corresponds to a non-scalar + //response and distributed parameters. Tempus::ForwardSensitivityIntegrator + //cannot return a LinearOp for DxDp. Therefore this case is not relevant here. + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + "\n Error! Piro::TransientSolver: DgDp = DERIV_LINEAR_OP is not supported with forward sensitivities!"); + } + //Cases 2 and 3 below correspond to scalar responses and scalar parameters. These can happen + //with dgdp = DERIV_MV_GRADIENT_FORM and dgpg = DERIV_MV_JACOBIAN_FORM. For + //DERIV_MV_GRADIENT_FORM, the map is the responses, and the columns are the parameters; for + //DERIV_MV_JACOBIAN_FORM, the map is the parameters, and the columns are the responses. + const RCP > dgdp_mv = dgdp_deriv.getMultiVector(); + if (Teuchos::nonnull(dgdp_mv)) { + if (dgdp_deriv.getMultiVectorOrientation() == Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM) { + //Case 2: DgDp = DERIV_MV_GRADIENT_FORM, DgDx is MV, DxDp is MV. + //This corresponds to a scalar response and distributed parameters + //[dgdp_mv]^T = [dx/dp_mv]^T*dg/dx_mv + [dg/dp_mv]^T + //Note: Gradient form stores transpose of derivative in dgdx_mv (DERIV_MV_GRADIENT_FORM is transposed!) + Thyra::apply(*dxdp_mv, Thyra::TRANS, *dgdx_mv, dgdp_mv.ptr(), Teuchos::ScalarTraits::one(), + Teuchos::ScalarTraits::one()); + } + else { + //Case 3: DgDp = DERIV_MV_JACOBIAN_FORM (the alternate to DERIV_MV_GRADIENT_FORM for getMultiVectorOrientation), + //DgDx = DERIV_LINEAR_OP and DxDp is MV. Note that DgDx implementes a DERIV_LINEAR_OP for MVs, + //so there is no contradiction here in the type. + //dgdp_mv = dg/dx_op*dx/dp_mv + dg/dp_mv + Thyra::apply(*dgdx_op, Thyra::NOTRANS, *dxdp_mv, dgdp_mv.ptr(), Teuchos::ScalarTraits::one(), + Teuchos::ScalarTraits::one()); + } + } + } + } + } + } + break; + } + case ADJOINT: //adjoint sensitivities - not yet implemented + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + "\n Error! Piro::TransientSolver: adjoint sentivities (Sensitivity Method = " + << "Adjoint) are not yet supported! Please set 'Sensitivity Method' to 'None' or 'Forward'.\n"); + break; + } + } } diff --git a/packages/piro/src/TriKota_MPDirectApplicInterface.cpp b/packages/piro/src/TriKota_MPDirectApplicInterface.cpp index b6b4e9316f85..e29a99022cec 100644 --- a/packages/piro/src/TriKota_MPDirectApplicInterface.cpp +++ b/packages/piro/src/TriKota_MPDirectApplicInterface.cpp @@ -56,7 +56,7 @@ MPDirectApplicInterface( model(model_), p_index(p_index_), g_index(g_index_), - orientation(EEME::DERIV_MV_BY_COL) + orientation(EEME::DERIV_MV_JACOBIAN_FORM) { out = Teuchos::VerboseObjectBase::getDefaultOStream(); @@ -86,16 +86,16 @@ MPDirectApplicInterface( EEME::DerivativeSupport supportDgDp = model->createOutArgs().supports(EEME::OUT_ARG_DgDp_mp, g_index, p_index); supportsSensitivities = - supportDgDp.supports(EEME::DERIV_TRANS_MV_BY_ROW) || - supportDgDp.supports(EEME::DERIV_MV_BY_COL); + supportDgDp.supports(EEME::DERIV_MV_GRADIENT_FORM) || + supportDgDp.supports(EEME::DERIV_MV_JACOBIAN_FORM); if (supportsSensitivities) { *out << "TriKota:: ModeEval supports gradients calculation." << std::endl; - if (supportDgDp.supports(EEME::DERIV_TRANS_MV_BY_ROW)) { - orientation = EEME::DERIV_TRANS_MV_BY_ROW; + if (supportDgDp.supports(EEME::DERIV_MV_GRADIENT_FORM)) { + orientation = EEME::DERIV_MV_GRADIENT_FORM; model_dgdp = model->create_p_mv_mp(p_index, numResponses); } else { - orientation = EEME::DERIV_MV_BY_COL; + orientation = EEME::DERIV_MV_JACOBIAN_FORM; model_dgdp = model->create_g_mv_mp(g_index, numParameters); } } @@ -245,7 +245,7 @@ wait_local_evaluations(Dakota::PRPQueue& prp_queue) for (unsigned int j=0; j +#include + +using namespace Teuchos; +using namespace Piro; +using namespace Piro::Test; + + +namespace Thyra { + typedef ModelEvaluatorBase MEB; +} // namespace Thyra + +// Setup support + +// Floating point tolerance +const double tol = 1.0e-8; + +void test_sincos_fsa(const bool use_combined_method, + const bool use_dfdp_as_tangent, + Teuchos::FancyOStream &out, bool &success) +{ + const std::string sens_method_string = "Forward"; + std::string outfile_name; + std::string errfile_name; + + if (use_combined_method == true) { + if (use_dfdp_as_tangent == true) { + outfile_name = "Tempus_BackwardEuler_SinCos_Sens_Combined_FSA_Tangent.dat"; + } + else { + outfile_name = "Tempus_BackwardEuler_SinCos_Sens_Combined_FSA.dat"; + } + } + else { + if (use_dfdp_as_tangent == true) { + outfile_name = "Tempus_BackwardEuler_SinCos_Sens_Staggered_FSA_Tangent.dat"; + } + else { + outfile_name = "Tempus_BackwardEuler_SinCos_Sens_Staggered_FSA.dat"; + } + } + + std::vector StepSize; + std::vector ErrorNorm; + const int nTimeStepSizes = 7; + double dt = 0.2; + double order = 0.0; + Teuchos::RCP > comm = + Teuchos::DefaultComm::getComm(); + Teuchos::RCP my_out = + Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); + my_out->setProcRankAndSize(comm->getRank(), comm->getSize()); + my_out->setOutputToRootOnly(0); + + SENS_METHOD sens_method; + if (sens_method_string == "None") sens_method = Piro::NONE; + else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; + else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; + + for (int n=0; n pList = + getParametersFromXmlFile("input_Tempus_BackwardEuler_SinCos.xml"); + + // Setup the SinCosModel + RCP scm_pl = sublist(pList, "SinCosModel", true); + scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent); + RCP model = Teuchos::rcp(new SinCosModel(scm_pl)); + + dt /= 2; + + // Setup sensitivities + RCP tempus_pl = sublist(pList, "Tempus", true); + ParameterList& sens_pl = tempus_pl->sublist("Sensitivities"); + if (use_combined_method) + sens_pl.set("Sensitivity Method", "Combined"); + else { + sens_pl.set("Sensitivity Method", "Staggered"); + sens_pl.set("Reuse State Linear Solver", true); + } + sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent); + + // Setup the Integrator and reset initial time step + tempus_pl->sublist("Default Integrator") + .sublist("Time Step Control").set("Initial Time Step", dt); + Teuchos::RCP > integrator + = Teuchos::rcp(new Piro::TempusIntegrator(tempus_pl, model, sens_method)); + order = integrator->getStepper()->getOrder(); + + // Initial Conditions + double t0 = tempus_pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Initial Time"); + double tfinal = tempus_pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Final Time"); + RCP > x0 = + model->getExactSolution(t0).get_x(); + const int num_param = model->get_p_space(0)->dim(); + RCP > DxDp0 = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, t0).get_x())); + } + integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null, + DxDp0, Teuchos::null, Teuchos::null); + + const RCP > stepSolver = Teuchos::null; + + RCP stepper_pl = Teuchos::rcp(&(tempus_pl->sublist("Default Stepper")), false); + + RCP > sf = Teuchos::rcp(new Tempus::StepperFactory()); + const RCP > stepper = sf->createStepper(stepper_pl, model); + const RCP > tempus_solver = + rcp(new TempusSolver(integrator, stepper, stepSolver, model, tfinal, sens_method_string)); + + const Thyra::MEB::InArgs inArgs = tempus_solver->getNominalValues(); + Thyra::MEB::OutArgs outArgs = tempus_solver->createOutArgs(); + const int solutionResponseIndex = tempus_solver->Ng() - 1; + const int parameterIndex = 0; + const Thyra::MEB::Derivative dxdp_deriv = + Thyra::create_DgDp_mv(*tempus_solver, solutionResponseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dxdp = dxdp_deriv.getMultiVector(); + outArgs.set_DgDp(solutionResponseIndex, parameterIndex, dxdp_deriv); + + //Integrate in time + tempus_solver->evalModel(inArgs, outArgs); + + // Test if at 'Final Time' + double time = integrator->getTime(); + double timeFinal = tempus_pl->sublist("Default Integrator") + .sublist("Time Step Control").get("Final Time"); + TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14); + + // Time-integrated solution and the exact solution + RCP > x = integrator->getX(); + RCP > DxDp = integrator->getDxDp(); + RCP > x_exact = + model->getExactSolution(time).get_x(); + RCP > DxDp_exact = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; icol(i).ptr(), + *(model->getExactSensSolution(i, time).get_x())); + } + // Plot sample solution and exact solution + if (comm->getRank() == 0 && n == nTimeStepSizes-1) { + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + + std::ofstream ftmp(outfile_name); + RCP > solutionHistory = + integrator->getSolutionHistory(); + RCP< Thyra::MultiVectorBase > DxDp_exact_plot = + Thyra::createMembers(model->get_x_space(), num_param); + for (int i=0; igetNumStates(); i++) { + RCP > solutionState = (*solutionHistory)[i]; + double time_i = solutionState->getTime(); + RCP x_prod_plot = + Teuchos::rcp_dynamic_cast(solutionState->getX()); + RCP > x_plot = + x_prod_plot->getMultiVector()->col(0); + RCP > DxDp_plot = + x_prod_plot->getMultiVector()->subView(Teuchos::Range1D(1,num_param)); + RCP > x_exact_plot = + model->getExactSolution(time_i).get_x(); + for (int j=0; jcol(j).ptr(), + *(model->getExactSensSolution(j, time_i).get_x())); + } + ftmp << std::fixed << std::setprecision(7) + << time_i + << std::setw(11) << get_ele(*(x_plot), 0) + << std::setw(11) << get_ele(*(x_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1); + } + ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) + << std::setw(11) << get_ele(*(x_exact_plot), 1); + for (int j=0; jcol(j)), 0) + << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1); + } + ftmp << std::endl; + } + ftmp.close(); + } + // Calculate the error + RCP > xdiff = x->clone_v(); + RCP > DxDpdiff = DxDp->clone_mv(); + Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x)); + Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp); + StepSize.push_back(dt); + double L2norm = Thyra::norm_2(*xdiff); + L2norm *= L2norm; + Teuchos::Array L2norm_DxDp(num_param); + Thyra::norms_2(*DxDpdiff, L2norm_DxDp()); + for (int i=0; igetRank() == 0) { + std::ofstream ftmp(errfile_name); + double error0 = 0.8*ErrorNorm[0]; + for (int n=0; n +#include + +using namespace Teuchos; +using namespace Piro; +using namespace Piro::Test; + + +namespace Thyra { + typedef ModelEvaluatorBase MEB; +} // namespace Thyra + +// Setup support + +const RCP epetraModelNew() +{ +#ifdef HAVE_MPI + const MPI_Comm comm = MPI_COMM_WORLD; +#else /*HAVE_MPI*/ + const int comm = 0; +#endif /*HAVE_MPI*/ + return rcp(new MockModelEval_A(comm)); +} + +const RCP > thyraModelNew(const RCP &epetraModel) +{ + const RCP > lowsFactory(new Thyra::AmesosLinearOpWithSolveFactory); + return epetraModelEvaluator(epetraModel, lowsFactory); +} + +RCP > defaultModelNew() +{ + return thyraModelNew(epetraModelNew()); +} + +const RCP > solverNew( + const RCP > &thyraModel, + double finalTime, + const std::string sens_method_string) +{ + const RCP tempusPL(new ParameterList("Tempus")); + tempusPL->set("Integrator Name", "Demo Integrator"); + tempusPL->sublist("Demo Integrator").set("Integrator Type", "Integrator Basic"); + tempusPL->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper"); + tempusPL->sublist("Demo Integrator").sublist("Solution History").set("Storage Type", "Unlimited"); + tempusPL->sublist("Demo Integrator").sublist("Solution History").set("Storage Limit", 20); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Initial Time", 0.0); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Final Time", finalTime); + tempusPL->sublist("Demo Stepper").set("Stepper Type", "Backward Euler"); + tempusPL->sublist("Demo Stepper").set("Zero Initial Guess", false); + tempusPL->sublist("Demo Stepper").set("Solver Name", "Demo Solver"); + tempusPL->sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); + SENS_METHOD sens_method; + if (sens_method_string == "None") sens_method = Piro::NONE; + else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; + else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; + Teuchos::RCP > integrator + = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, thyraModel, sens_method)); + auto x0 = thyraModel->getNominalValues().get_x(); + const int num_param = thyraModel->get_p_space(0)->dim(); + RCP > DxDp0 = + Thyra::createMembers(thyraModel->get_x_space(), num_param); + DxDp0->assign(0.0); + integrator->initializeSolutionHistory(0.0, x0, Teuchos::null, Teuchos::null, + DxDp0, Teuchos::null, Teuchos::null); + const RCP > stepSolver = Teuchos::null; + + RCP stepperPL = Teuchos::rcp(&(tempusPL->sublist("Demo Stepper")), false); + + RCP > sf = Teuchos::rcp(new Tempus::StepperFactory()); + const RCP > stepper = sf->createStepper(stepperPL, thyraModel); + return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, finalTime, sens_method_string)); +} + +const RCP > solverNew( + const RCP > &thyraModel, + double initialTime, + double finalTime, + const RCP > &observer, + const std::string sens_method_string) +{ + const RCP tempusPL(new ParameterList("Tempus")); + tempusPL->set("Integrator Name", "Demo Integrator"); + tempusPL->sublist("Demo Integrator").set("Integrator Type", "Integrator Basic"); + tempusPL->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper"); + tempusPL->sublist("Demo Integrator").sublist("Solution History").set("Storage Type", "Unlimited"); + tempusPL->sublist("Demo Integrator").sublist("Solution History").set("Storage Limit", 20); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Initial Time", initialTime); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Final Time", finalTime); + tempusPL->sublist("Demo Stepper").set("Stepper Type", "Backward Euler"); + tempusPL->sublist("Demo Stepper").set("Zero Initial Guess", false); + tempusPL->sublist("Demo Stepper").set("Solver Name", "Demo Solver"); + tempusPL->sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); + SENS_METHOD sens_method; + if (sens_method_string == "None") sens_method = Piro::NONE; + else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; + else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; + Teuchos::RCP > integrator + = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, thyraModel, sens_method)); + auto x0 = thyraModel->getNominalValues().get_x(); + const int num_param = thyraModel->get_p_space(0)->dim(); + RCP > DxDp0 = + Thyra::createMembers(thyraModel->get_x_space(), num_param); + DxDp0->assign(0.0); + integrator->initializeSolutionHistory(initialTime, x0, Teuchos::null, Teuchos::null, + DxDp0, Teuchos::null, Teuchos::null); + const RCP > solutionHistory = integrator->getSolutionHistory(); + const RCP > timeStepControl = integrator->getTimeStepControl(); + + const Teuchos::RCP > tempusObserver = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter(solutionHistory, timeStepControl, observer, false, false, sens_method)); + integrator->setObserver(tempusObserver); + const RCP > stepSolver = Teuchos::null; + RCP stepperPL = Teuchos::rcp(&(tempusPL->sublist("Demo Stepper")), false); + RCP > sf = Teuchos::rcp(new Tempus::StepperFactory()); + const RCP > stepper = sf->createStepper(stepperPL, thyraModel); + + return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, initialTime, finalTime, sens_method_string)); +} + +const RCP > solverNew( + const RCP > &thyraModel, + double initialTime, + double finalTime, + double fixedTimeStep, + const RCP > &observer, + const std::string sens_method_string) +{ + const RCP tempusPL(new ParameterList("Tempus")); + tempusPL->set("Integrator Name", "Demo Integrator"); + tempusPL->sublist("Demo Integrator").set("Integrator Type", "Integrator Basic"); + tempusPL->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper"); + tempusPL->sublist("Demo Integrator").sublist("Solution History").set("Storage Type", "Unlimited"); + tempusPL->sublist("Demo Integrator").sublist("Solution History").set("Storage Limit", 20); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Initial Time", initialTime); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Final Time", finalTime); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Minimum Time Step", fixedTimeStep); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Initial Time Step", fixedTimeStep); + tempusPL->sublist("Demo Integrator").sublist("Time Step Control").set("Maximum Time Step", fixedTimeStep); + tempusPL->sublist("Demo Stepper").set("Stepper Type", "Backward Euler"); + tempusPL->sublist("Demo Stepper").set("Zero Initial Guess", false); + tempusPL->sublist("Demo Stepper").set("Solver Name", "Demo Solver"); + tempusPL->sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); + SENS_METHOD sens_method; + if (sens_method_string == "None") sens_method = Piro::NONE; + else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; + else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; + Teuchos::RCP > integrator + = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, thyraModel, sens_method)); + const RCP > solutionHistory = integrator->getSolutionHistory(); + const RCP > timeStepControl = integrator->getTimeStepControl(); + + const Teuchos::RCP > tempusObserver = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter(solutionHistory, timeStepControl, observer)); + integrator->setObserver(tempusObserver); + const RCP > stepSolver = Teuchos::null; + RCP stepperPL = Teuchos::rcp(&(tempusPL->sublist("Demo Stepper")), false); + RCP > sf = Teuchos::rcp(new Tempus::StepperFactory()); + const RCP > stepper = sf->createStepper(stepperPL, thyraModel); + + return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, initialTime, finalTime, sens_method_string)); +} + + +Thyra::ModelEvaluatorBase::InArgs getStaticNominalValues(const Thyra::ModelEvaluator &model) +{ + Thyra::ModelEvaluatorBase::InArgs result = model.getNominalValues(); + if (result.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) { + result.set_x_dot(Teuchos::null); + } + return result; +} + + +// Floating point tolerance +const double tol = 1.0e-8; + +TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_DefaultSolutionSensitivity) +{ + Teuchos::RCP fos = + Teuchos::VerboseObjectBase::getDefaultOStream(); + + const RCP > model = defaultModelNew(); + const double finalTime = 0.0; + + const RCP > solver = solverNew(model, finalTime, "Forward"); + + const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); + + Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); + const int solutionResponseIndex = solver->Ng() - 1; + const int parameterIndex = 0; + const Thyra::MEB::Derivative dxdp_deriv = + Thyra::create_DgDp_mv(*solver, solutionResponseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dxdp = dxdp_deriv.getMultiVector(); + outArgs.set_DgDp(solutionResponseIndex, parameterIndex, dxdp_deriv); + + solver->evalModel(inArgs, outArgs); + + // Test if at 'Final Time' + double time = solver->getPiroTempusIntegrator()->getTime(); + TEST_FLOATING_EQUALITY(time, 0.0, 1.0e-14); + + const Array > expected = tuple( + Array(tuple(0.0, 0.0, 0.0, 0.0)), + Array(tuple(0.0, 0.0, 0.0, 0.0))); + RCP > DxDp = solver->getPiroTempusIntegrator()->getDxDp(); + TEST_EQUALITY(DxDp->domain()->dim(), expected.size()); + + for (int i = 0; i < expected.size(); ++i) { + TEST_EQUALITY(DxDp->range()->dim(), expected[i].size()); + const Array actual = arrayFromVector(*DxDp->col(i)); + TEST_COMPARE_FLOATING_ARRAYS(actual, expected[i], tol); + } +} + +TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_DefaultSolutionSensitivityOp) +{ + Teuchos::RCP fos = + Teuchos::VerboseObjectBase::getDefaultOStream(); + + const RCP > model = defaultModelNew(); + const double finalTime = 0.0; + + const RCP > solver = solverNew(model, finalTime, "Forward"); + + const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); + + Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); + const int solutionResponseIndex = solver->Ng() - 1; + const int parameterIndex = 0; + const Thyra::MEB::Derivative dxdp_deriv = + solver->create_DgDp_op(solutionResponseIndex, parameterIndex); + const RCP > dxdp = dxdp_deriv.getLinearOp(); + outArgs.set_DgDp(solutionResponseIndex, parameterIndex, dxdp_deriv); + + solver->evalModel(inArgs, outArgs); + + // Test if at 'Final Time' + double time = solver->getPiroTempusIntegrator()->getTime(); + TEST_FLOATING_EQUALITY(time, 0.0, 1.0e-14); + + const Array > expected = tuple( + Array(tuple(0.0, 0.0, 0.0, 0.0)), + Array(tuple(0.0, 0.0, 0.0, 0.0))); + RCP > DxDp = solver->getPiroTempusIntegrator()->getDxDp(); + TEST_EQUALITY(DxDp->domain()->dim(), expected.size()); + for (int i = 0; i < expected.size(); ++i) { + TEST_EQUALITY(dxdp->range()->dim(), expected[i].size()); + const Array actual = arrayFromLinOp(*DxDp, i); + TEST_COMPARE_FLOATING_ARRAYS(actual, expected[i], tol); + } +} + +TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_DefaultResponseSensitivity) +{ + const RCP > model = defaultModelNew(); + + const int responseIndex = 0; + const int parameterIndex = 0; + + const Thyra::MEB::Derivative dgdp_deriv_expected = + Thyra::create_DgDp_mv(*model, responseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dgdp_expected = dgdp_deriv_expected.getMultiVector(); + { + const Thyra::MEB::InArgs modelInArgs = getStaticNominalValues(*model); + Thyra::MEB::OutArgs modelOutArgs = model->createOutArgs(); + modelOutArgs.set_DgDp(responseIndex, parameterIndex, dgdp_deriv_expected); + model->evalModel(modelInArgs, modelOutArgs); + } + + const double finalTime = 0.0; + const RCP > solver = solverNew(model, finalTime, "Forward"); + + const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); + + Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); + const Thyra::MEB::Derivative dgdp_deriv = + Thyra::create_DgDp_mv(*solver, responseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dgdp = dgdp_deriv.getMultiVector(); + outArgs.set_DgDp(responseIndex, parameterIndex, dgdp_deriv); + + solver->evalModel(inArgs, outArgs); + + TEST_EQUALITY(dgdp->domain()->dim(), dgdp_expected->domain()->dim()); + TEST_EQUALITY(dgdp->range()->dim(), dgdp_expected->range()->dim()); + for (int i = 0; i < dgdp_expected->domain()->dim(); ++i) { + const Array actual = arrayFromVector(*dgdp->col(i)); + const Array expected = arrayFromVector(*dgdp_expected->col(i)); + TEST_COMPARE_FLOATING_ARRAYS(actual, expected, tol); + } +} + +TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_DefaultResponseSensitivity_NoDgDxMv) +{ + Teuchos::RCP fos = + Teuchos::VerboseObjectBase::getDefaultOStream(); + + const RCP > model( + new WeakenedModelEvaluator_NoDgDxMv(defaultModelNew())); + + const int responseIndex = 0; + const int parameterIndex = 0; + + const Thyra::MEB::Derivative dgdp_deriv_expected = + Thyra::create_DgDp_mv(*model, responseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dgdp_expected = dgdp_deriv_expected.getMultiVector(); + { + const Thyra::MEB::InArgs modelInArgs = getStaticNominalValues(*model); + Thyra::MEB::OutArgs modelOutArgs = model->createOutArgs(); + modelOutArgs.set_DgDp(responseIndex, parameterIndex, dgdp_deriv_expected); + model->evalModel(modelInArgs, modelOutArgs); + } + + const double finalTime = 0.0; + const RCP > solver = solverNew(model, finalTime, "Forward"); + + const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); + + Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); + const Thyra::MEB::Derivative dgdp_deriv = + Thyra::create_DgDp_mv(*solver, responseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dgdp = dgdp_deriv.getMultiVector(); + outArgs.set_DgDp(responseIndex, parameterIndex, dgdp_deriv); + + solver->evalModel(inArgs, outArgs); + + TEST_EQUALITY(dgdp->domain()->dim(), dgdp_expected->domain()->dim()); + TEST_EQUALITY(dgdp->range()->dim(), dgdp_expected->range()->dim()); + for (int i = 0; i < dgdp_expected->domain()->dim(); ++i) { + const Array actual = arrayFromVector(*dgdp->col(i)); + const Array expected = arrayFromVector(*dgdp_expected->col(i)); + TEST_COMPARE_FLOATING_ARRAYS(actual, expected, tol); + } +} + + +TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_ResponseAndDefaultSensitivities) +{ + Teuchos::RCP fos = + Teuchos::VerboseObjectBase::getDefaultOStream(); + + const RCP > model = defaultModelNew(); + + const int responseIndex = 0; + const int parameterIndex = 0; + + const RCP > expectedResponse = + Thyra::createMember(model->get_g_space(responseIndex)); + { + const Thyra::MEB::InArgs modelInArgs = getStaticNominalValues(*model); + Thyra::MEB::OutArgs modelOutArgs = model->createOutArgs(); + modelOutArgs.set_g(responseIndex, expectedResponse); + model->evalModel(modelInArgs, modelOutArgs); + } + + const Thyra::MEB::Derivative dgdp_deriv_expected = + Thyra::create_DgDp_mv(*model, responseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dgdp_expected = dgdp_deriv_expected.getMultiVector(); + { + const Thyra::MEB::InArgs modelInArgs = getStaticNominalValues(*model); + Thyra::MEB::OutArgs modelOutArgs = model->createOutArgs(); + modelOutArgs.set_DgDp(responseIndex, parameterIndex, dgdp_deriv_expected); + model->evalModel(modelInArgs, modelOutArgs); + } + + const double finalTime = 0.0; + const RCP > solver = solverNew(model, finalTime, "Forward"); + + const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); + + Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); + + // Requesting response + const RCP > response = + Thyra::createMember(solver->get_g_space(responseIndex)); + outArgs.set_g(responseIndex, response); + + // Requesting response sensitivity + const Thyra::MEB::Derivative dgdp_deriv = + Thyra::create_DgDp_mv(*solver, responseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dgdp = dgdp_deriv.getMultiVector(); + outArgs.set_DgDp(responseIndex, parameterIndex, dgdp_deriv); + + // Run solver + solver->evalModel(inArgs, outArgs); + + // Checking response + { + const Array expected = arrayFromVector(*expectedResponse); + const Array actual = arrayFromVector(*response); + TEST_COMPARE_FLOATING_ARRAYS(actual, expected, tol); + } + + // Checking sensitivity + { + TEST_EQUALITY(dgdp->domain()->dim(), dgdp_expected->domain()->dim()); + TEST_EQUALITY(dgdp->range()->dim(), dgdp_expected->range()->dim()); + for (int i = 0; i < dgdp_expected->domain()->dim(); ++i) { + const Array actual = arrayFromVector(*dgdp->col(i)); + const Array expected = arrayFromVector(*dgdp_expected->col(i)); + TEST_COMPARE_FLOATING_ARRAYS(actual, expected, tol); + } + } +} + +TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveInitialConditionWhenSensitivitiesRequested) +{ + Teuchos::RCP fos = + Teuchos::VerboseObjectBase::getDefaultOStream(); + + const RCP > model = defaultModelNew(); + const RCP > observer(new MockObserver); + const double timeStamp = 2.0; + + const RCP > solver = solverNew(model, timeStamp, timeStamp, observer, "Forward"); + + const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); + Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); + { + const int solutionResponseIndex = solver->Ng() - 1; + const int parameterIndex = 0; + const Thyra::MEB::Derivative dxdp_deriv = + Thyra::create_DgDp_mv(*solver, solutionResponseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); + const RCP > dxdp = dxdp_deriv.getMultiVector(); + outArgs.set_DgDp(solutionResponseIndex, parameterIndex, dxdp_deriv); + } + solver->evalModel(inArgs, outArgs); + + { + const RCP > solution = + observer->lastSolution(); + + const RCP > initialCondition = + model->getNominalValues().get_x(); + + TEST_COMPARE_FLOATING_ARRAYS( + arrayFromVector(*solution), + arrayFromVector(*initialCondition), + tol); + } + + TEST_FLOATING_EQUALITY(observer->lastStamp(), timeStamp, tol); +} + +#endif /* HAVE_PIRO_TEMPUS */ diff --git a/packages/piro/test/Piro_TempusSolver_UnitTests.cpp b/packages/piro/test/Piro_TempusSolver_UnitTests.cpp index 15c15d147b89..47a3e59b1aec 100644 --- a/packages/piro/test/Piro_TempusSolver_UnitTests.cpp +++ b/packages/piro/test/Piro_TempusSolver_UnitTests.cpp @@ -47,7 +47,6 @@ #ifdef HAVE_PIRO_TEMPUS #include "Piro_TempusSolver.hpp" #include "Tempus_StepperFactory.hpp" -//#include "Tempus_StepperBackwardEuler.hpp" #include "Piro_ObserverToTempusIntegrationObserverAdapter.hpp" #ifdef HAVE_PIRO_NOX @@ -57,6 +56,7 @@ #include "Piro_Test_ThyraSupport.hpp" #include "Piro_Test_WeakenedModelEvaluator.hpp" #include "Piro_Test_MockObserver.hpp" +#include "Piro_TempusIntegrator.hpp" #include "MockModelEval_A.hpp" @@ -75,6 +75,7 @@ #include + using namespace Teuchos; using namespace Piro; using namespace Piro::Test; @@ -108,7 +109,8 @@ RCP > defaultModelNew() const RCP > solverNew( const RCP > &thyraModel, - double finalTime) + double finalTime, + const std::string sens_method_string) { const RCP tempusPL(new ParameterList("Tempus")); tempusPL->set("Integrator Name", "Demo Integrator"); @@ -122,22 +124,27 @@ const RCP > solverNew( tempusPL->sublist("Demo Stepper").set("Zero Initial Guess", false); tempusPL->sublist("Demo Stepper").set("Solver Name", "Demo Solver"); tempusPL->sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); - Teuchos::RCP > integrator = Tempus::integratorBasic(tempusPL, thyraModel); + SENS_METHOD sens_method; + if (sens_method_string == "None") sens_method = Piro::NONE; + else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; + else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; + Teuchos::RCP > integrator + = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, thyraModel, sens_method)); const RCP > stepSolver = Teuchos::null; RCP stepperPL = Teuchos::rcp(&(tempusPL->sublist("Demo Stepper")), false); RCP > sf = Teuchos::rcp(new Tempus::StepperFactory()); const RCP > stepper = sf->createStepper(stepperPL, thyraModel); - //const RCP > stepper = rcp(new Tempus::StepperBackwardEuler(thyraModel, stepperPL)); - return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, finalTime)); + return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, finalTime, sens_method_string)); } const RCP > solverNew( const RCP > &thyraModel, double initialTime, double finalTime, - const RCP > &observer) + const RCP > &observer, + const std::string sens_method_string) { const RCP tempusPL(new ParameterList("Tempus")); tempusPL->set("Integrator Name", "Demo Integrator"); @@ -151,7 +158,12 @@ const RCP > solverNew( tempusPL->sublist("Demo Stepper").set("Zero Initial Guess", false); tempusPL->sublist("Demo Stepper").set("Solver Name", "Demo Solver"); tempusPL->sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); - Teuchos::RCP > integrator = Tempus::integratorBasic(tempusPL, thyraModel); + SENS_METHOD sens_method; + if (sens_method_string == "None") sens_method = Piro::NONE; + else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; + else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; + Teuchos::RCP > integrator + = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, thyraModel, sens_method)); const RCP > solutionHistory = integrator->getSolutionHistory(); const RCP > timeStepControl = integrator->getTimeStepControl(); @@ -161,9 +173,8 @@ const RCP > solverNew( RCP stepperPL = Teuchos::rcp(&(tempusPL->sublist("Demo Stepper")), false); RCP > sf = Teuchos::rcp(new Tempus::StepperFactory()); const RCP > stepper = sf->createStepper(stepperPL, thyraModel); - //const RCP > stepper = rcp(new Tempus::StepperBackwardEuler(thyraModel, stepperPL)); - return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, initialTime, finalTime)); + return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, initialTime, finalTime, sens_method_string)); } const RCP > solverNew( @@ -171,7 +182,8 @@ const RCP > solverNew( double initialTime, double finalTime, double fixedTimeStep, - const RCP > &observer) + const RCP > &observer, + const std::string sens_method_string) { const RCP tempusPL(new ParameterList("Tempus")); tempusPL->set("Integrator Name", "Demo Integrator"); @@ -188,19 +200,22 @@ const RCP > solverNew( tempusPL->sublist("Demo Stepper").set("Zero Initial Guess", false); tempusPL->sublist("Demo Stepper").set("Solver Name", "Demo Solver"); tempusPL->sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); - Teuchos::RCP > integrator = Tempus::integratorBasic(tempusPL, thyraModel); + Teuchos::RCP > integrator = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, thyraModel)); const RCP > solutionHistory = integrator->getSolutionHistory(); const RCP > timeStepControl = integrator->getTimeStepControl(); - - const Teuchos::RCP > tempusObserver = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter(solutionHistory, timeStepControl, observer)); + SENS_METHOD sens_method; + if (sens_method_string == "None") sens_method = Piro::NONE; + else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; + else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; + const Teuchos::RCP > tempusObserver + = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter(solutionHistory, timeStepControl, observer, sens_method)); integrator->setObserver(tempusObserver); const RCP > stepSolver = Teuchos::null; RCP stepperPL = Teuchos::rcp(&(tempusPL->sublist("Demo Stepper")), false); RCP > sf = Teuchos::rcp(new Tempus::StepperFactory()); const RCP > stepper = sf->createStepper(stepperPL, thyraModel); - //const RCP > stepper = rcp(new Tempus::StepperBackwardEuler(thyraModel, stepperPL)); - return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, initialTime, finalTime)); + return rcp(new TempusSolver(integrator, stepper, stepSolver, thyraModel, initialTime, finalTime, sens_method_string)); } @@ -222,7 +237,7 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_Solution) const RCP > model = defaultModelNew(); const double finalTime = 0.0; - const RCP > solver = solverNew(model, finalTime); + const RCP > solver = solverNew(model, finalTime, "None"); const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); @@ -260,7 +275,7 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_Response) } const double finalTime = 0.0; - const RCP > solver = solverNew(model, finalTime); + const RCP > solver = solverNew(model, finalTime, "None"); const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); @@ -282,7 +297,7 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_NoDfDpMv_NoSensitivity) new WeakenedModelEvaluator_NoDfDpMv(defaultModelNew())); const double finalTime = 0.0; - const RCP > solver = solverNew(model, finalTime); + const RCP > solver = solverNew(model, finalTime, "None"); const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); @@ -304,7 +319,7 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_NoDgDp_NoResponseSensitivity) new WeakenedModelEvaluator_NoDgDp(defaultModelNew())); const double finalTime = 0.0; - const RCP > solver = solverNew(model, finalTime); + const RCP > solver = solverNew(model, finalTime, "None"); const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); @@ -324,10 +339,9 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveInitialCondition) { const RCP > model = defaultModelNew(); const RCP > observer(new MockObserver); - //IKT, FIXME: set to 2.0 instead of 0.0 -- does not work correctly in this case - const double timeStamp = 0.0; + const double timeStamp = 2.0; - const RCP > solver = solverNew(model, timeStamp, timeStamp, observer); + const RCP > solver = solverNew(model, timeStamp, timeStamp, observer, "None"); const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); const Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); @@ -359,7 +373,7 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveFinalSolution) const double timeStepSize = 0.05; const RCP > solver = - solverNew(model, initialTime, finalTime, timeStepSize, observer); + solverNew(model, initialTime, finalTime, timeStepSize, observer, "None"); const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); @@ -379,5 +393,4 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveFinalSolution) TEST_FLOATING_EQUALITY(observer->lastStamp(), finalTime, tol); } - #endif /* HAVE_PIRO_TEMPUS */ diff --git a/packages/piro/test/SinCosModel.cpp b/packages/piro/test/SinCosModel.cpp new file mode 100644 index 000000000000..b8eed365ee13 --- /dev/null +++ b/packages/piro/test/SinCosModel.cpp @@ -0,0 +1,630 @@ +// @HEADER +// ************************************************************************ +// +// Piro: Strategy package for embedded analysis capabilitites +// Copyright (2010) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Andy Salinger (agsalin@sandia.gov), Sandia +// National Laboratories. +// +// ************************************************************************ +// @HEADER + +#include "SinCosModel.hpp" + +#include "Teuchos_StandardParameterEntryValidators.hpp" + +#include "Thyra_DefaultSpmdVectorSpace.hpp" +#include "Thyra_DetachedVectorView.hpp" +#include "Thyra_DetachedMultiVectorView.hpp" +#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp" +#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp" +#include "Thyra_DefaultLinearOpSource.hpp" +#include "Thyra_VectorStdOps.hpp" +#include "Thyra_MultiVectorStdOps.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" + +#include + +using Teuchos::RCP; +using Teuchos::rcp; + +SinCosModel:: +SinCosModel(Teuchos::RCP pList_) +{ + isInitialized_ = false; + dim_ = 2; + Np_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp) + np_ = 3; // Number of parameters in this vector (3) + Ng_ = 1; // Number of observation functions (1) + ng_ = dim_; // Number of elements in this observation function ( == x ) + acceptModelParams_ = false; + useDfDpAsTangent_ = false; + haveIC_ = true; + a_ = 0.0; + f_ = 1.0; + L_ = 1.0; + x0_ic_ = 0.0; + x1_ic_ = 1.0; + t0_ic_ = 0.0; + + // Create x_space and f_space + x_space_ = Thyra::defaultSpmdVectorSpace(dim_); + f_space_ = Thyra::defaultSpmdVectorSpace(dim_); + // Create p_space and g_space + p_space_ = Thyra::defaultSpmdVectorSpace(np_); + g_space_ = Thyra::defaultSpmdVectorSpace(ng_); + + setParameterList(pList_); + + // Create DxDp product space + DxDp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_); +} + +// +// 06/24/09 tscoffe: +// These are the exact sensitivities for the problem assuming the initial conditions are specified as: +// s(0) = [1;0] s(1) = [0;b/L] s(2) = [0;-b*f/(L*L)] +// sdot(0) = [0;0] sdot(1) = [0;-3*f*f*b/(L*L*L)] sdot(2) = [0;3*f*f*f*b/(L*L*L*L)] +// + + +Thyra::ModelEvaluatorBase::InArgs +SinCosModel:: +getExactSolution(double t) const +{ + TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error, + "Error, setupInOutArgs_ must be called first!\n"); + Thyra::ModelEvaluatorBase::InArgs inArgs = inArgs_; + double exact_t = t; + inArgs.set_t(exact_t); + Teuchos::RCP > exact_x = createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView exact_x_view(*exact_x); + exact_x_view[0] = a_+b_*sin((f_/L_)*t+phi_); + exact_x_view[1] = b_*(f_/L_)*cos((f_/L_)*t+phi_); + } + inArgs.set_x(exact_x); + Teuchos::RCP > exact_x_dot = createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView exact_x_dot_view(*exact_x_dot); + exact_x_dot_view[0] = b_*(f_/L_)*cos((f_/L_)*t+phi_); + exact_x_dot_view[1] = -b_*(f_/L_)*(f_/L_)*sin((f_/L_)*t+phi_); + } + inArgs.set_x_dot(exact_x_dot); + return(inArgs); +} + +// +// 06/24/09 tscoffe: +// These are the exact sensitivities for the problem assuming the initial conditions are specified as: +// s(0) = [1;0] s(1) = [0;b/L] s(2) = [0;-b*f/(L*L)] +// sdot(0) = [0;0] sdot(1) = [0;-3*f*f*b/(L*L*L)] sdot(2) = [0;3*f*f*f*b/(L*L*L*L)] +// +Thyra::ModelEvaluatorBase::InArgs +SinCosModel:: +getExactSensSolution(int j, double t) const +{ + TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error, + "Error, setupInOutArgs_ must be called first!\n"); + Thyra::ModelEvaluatorBase::InArgs inArgs = inArgs_; + if (!acceptModelParams_) { + return inArgs; + } + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, np_ ); + double exact_t = t; + double b = b_; + double f = f_; + double L = L_; + double phi = phi_; + inArgs.set_t(exact_t); + Teuchos::RCP > exact_s = createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView exact_s_view(*exact_s); + if (j == 0) { // dx/da + exact_s_view[0] = 1.0; + exact_s_view[1] = 0.0; + } else if (j == 1) { // dx/df + exact_s_view[0] = (b/L)*t*cos((f/L)*t+phi); + exact_s_view[1] = (b/L)*cos((f/L)*t+phi)-(b*f*t/(L*L))*sin((f/L)*t+phi); + } else { // dx/dL + exact_s_view[0] = -(b*f*t/(L*L))*cos((f/L)*t+phi); + exact_s_view[1] = -(b*f/(L*L))*cos((f/L)*t+phi)+(b*f*f*t/(L*L*L))*sin((f/L)*t+phi); + } + } + inArgs.set_x(exact_s); + Teuchos::RCP > exact_s_dot = createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView exact_s_dot_view(*exact_s_dot); + if (j == 0) { // dxdot/da + exact_s_dot_view[0] = 0.0; + exact_s_dot_view[1] = 0.0; + } else if (j == 1) { // dxdot/df + exact_s_dot_view[0] = (b/L)*cos((f/L)*t+phi)-(b*f*t/(L*L))*sin((f/L)*t+phi); + exact_s_dot_view[1] = -(2.0*b*f/(L*L))*sin((f/L)*t+phi)-(b*f*f*t/(L*L*L))*cos((f/L)*t+phi); + } else { // dxdot/dL + exact_s_dot_view[0] = -(b*f/(L*L))*cos((f/L)*t+phi)+(b*f*f*t/(L*L*L))*sin((f/L)*t+phi); + exact_s_dot_view[1] = (2.0*b*f*f/(L*L*L))*sin((f/L)*t+phi)+(b*f*f*f*t/(L*L*L*L))*cos((f/L)*t+phi); + } + } + inArgs.set_x_dot(exact_s_dot); + return(inArgs); +} + + + +Teuchos::RCP > +SinCosModel:: +get_x_space() const +{ + return x_space_; +} + +Teuchos::RCP > +SinCosModel:: +get_f_space() const +{ + return f_space_; +} + +Thyra::ModelEvaluatorBase::InArgs +SinCosModel:: +getNominalValues() const +{ + TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error, + "Error, setupInOutArgs_ must be called first!\n"); + return nominalValues_; +} + +Teuchos::RCP > +SinCosModel:: +create_W() const +{ + using Teuchos::RCP; + RCP > W_factory = this->get_W_factory(); + RCP > matrix = this->create_W_op(); + { + // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to + // linearOpWithSolve because it ends up factoring the matrix during + // initialization, which it really shouldn't do, or I'm doing something + // wrong here. The net effect is that I get exceptions thrown in + // optimized mode due to the matrix being rank deficient unless I do this. + RCP > multivec = Teuchos::rcp_dynamic_cast >(matrix,true); + { + RCP > vec = Thyra::createMember(x_space_); + { + Thyra::DetachedVectorView vec_view( *vec ); + vec_view[0] = 0.0; + vec_view[1] = 1.0; + } + V_V(multivec->col(0).ptr(),*vec); + { + Thyra::DetachedVectorView vec_view( *vec ); + vec_view[0] = 1.0; + vec_view[1] = 0.0; + } + V_V(multivec->col(1).ptr(),*vec); + } + } + RCP > W = + Thyra::linearOpWithSolve(*W_factory, matrix); + return W; +} + +Teuchos::RCP > +SinCosModel:: +create_W_op() const +{ + Teuchos::RCP > matrix = Thyra::createMembers(x_space_, dim_); + return(matrix); +} + +Teuchos::RCP > +SinCosModel:: +get_W_factory() const +{ + Teuchos::RCP > W_factory = + Thyra::defaultSerialDenseLinearOpWithSolveFactory(); + return W_factory; +} + + +Thyra::ModelEvaluatorBase::InArgs +SinCosModel:: +createInArgs() const +{ + setupInOutArgs_(); + return inArgs_; +} + +Thyra::ModelEvaluatorBase::OutArgs +SinCosModel:: +createOutArgsImpl() const +{ + setupInOutArgs_(); + return outArgs_; +} + +void +SinCosModel:: +evalModelImpl( const Thyra::ModelEvaluatorBase::InArgs &inArgs, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs) const +{ + typedef Thyra::DefaultMultiVectorProductVector DMVPV; + using Teuchos::RCP; + using Thyra::VectorBase; + using Thyra::MultiVectorBase; + using Teuchos::rcp_dynamic_cast; + TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error, + "Error, setupInOutArgs_ must be called first!\n"); + + const RCP > x_in = inArgs.get_x().assert_not_null(); + Thyra::ConstDetachedVectorView x_in_view( *x_in ); + + //double t = inArgs.get_t(); + double a = a_; + double f = f_; + double L = L_; + if (acceptModelParams_) { + const RCP > p_in = + inArgs.get_p(0).assert_not_null(); + Thyra::ConstDetachedVectorView p_in_view( *p_in ); + a = p_in_view[0]; + f = p_in_view[1]; + L = p_in_view[2]; + } + RCP > DxDp_in, DxdotDp_in; + if (acceptModelParams_) { + if (inArgs.get_p(1) != Teuchos::null) + DxDp_in = + rcp_dynamic_cast(inArgs.get_p(1))->getMultiVector(); + if (inArgs.get_p(2) != Teuchos::null) + DxdotDp_in = + rcp_dynamic_cast(inArgs.get_p(2))->getMultiVector(); + } + + double beta = inArgs.get_beta(); + + const RCP > f_out = outArgs.get_f(); + const RCP > W_out = outArgs.get_W_op(); + RCP > DfDp_out; + if (acceptModelParams_) { + Thyra::ModelEvaluatorBase::Derivative DfDp = outArgs.get_DfDp(0); + DfDp_out = DfDp.getMultiVector(); + } + if (inArgs.get_x_dot().is_null()) { + + // Evaluate the Explicit ODE f(x,t) [= 0] + if (!is_null(f_out)) { + Thyra::DetachedVectorView f_out_view( *f_out ); + f_out_view[0] = x_in_view[1]; + f_out_view[1] = (f/L)*(f/L)*(a-x_in_view[0]); + } + if (!is_null(W_out)) { + RCP > matrix = + Teuchos::rcp_dynamic_cast >(W_out,true); + Thyra::DetachedMultiVectorView matrix_view( *matrix ); + matrix_view(0,0) = 0.0; // d(f0)/d(x0_n) + matrix_view(0,1) = +beta; // d(f0)/d(x1_n) + matrix_view(1,0) = -beta*(f/L)*(f/L); // d(f1)/d(x0_n) + matrix_view(1,1) = 0.0; // d(f1)/d(x1_n) + // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n) + } + if (!is_null(DfDp_out)) { + Thyra::DetachedMultiVectorView DfDp_out_view( *DfDp_out ); + DfDp_out_view(0,0) = 0.0; + DfDp_out_view(0,1) = 0.0; + DfDp_out_view(0,2) = 0.0; + DfDp_out_view(1,0) = (f/L)*(f/L); + DfDp_out_view(1,1) = (2.0*f/(L*L))*(a-x_in_view[0]); + DfDp_out_view(1,2) = -(2.0*f*f/(L*L*L))*(a-x_in_view[0]); + // Compute df/dp + (df/dx) * (dx/dp) + if (useDfDpAsTangent_ && !is_null(DxDp_in)) { + Thyra::ConstDetachedMultiVectorView DxDp( *DxDp_in ); + DfDp_out_view(0,0) += DxDp(1,0); + DfDp_out_view(0,1) += DxDp(1,1); + DfDp_out_view(0,2) += DxDp(1,2); + DfDp_out_view(1,0) += -(f/L)*(f/L) * DxDp(0,0); + DfDp_out_view(1,1) += -(f/L)*(f/L) * DxDp(0,1); + DfDp_out_view(1,2) += -(f/L)*(f/L) * DxDp(0,2); + } + } + } else { + + // Evaluate the implicit ODE f(xdot, x, t) [= 0] + RCP > x_dot_in; + x_dot_in = inArgs.get_x_dot().assert_not_null(); + double alpha = inArgs.get_alpha(); + if (!is_null(f_out)) { + Thyra::DetachedVectorView f_out_view( *f_out ); + Thyra::ConstDetachedVectorView x_dot_in_view( *x_dot_in ); + f_out_view[0] = x_dot_in_view[0] - x_in_view[1]; + f_out_view[1] = x_dot_in_view[1] - (f/L)*(f/L)*(a-x_in_view[0]); + } + if (!is_null(W_out)) { + RCP > matrix = + Teuchos::rcp_dynamic_cast >(W_out,true); + Thyra::DetachedMultiVectorView matrix_view( *matrix ); + matrix_view(0,0) = alpha; // d(f0)/d(x0_n) + matrix_view(0,1) = -beta; // d(f0)/d(x1_n) + matrix_view(1,0) = +beta*(f/L)*(f/L); // d(f1)/d(x0_n) + matrix_view(1,1) = alpha; // d(f1)/d(x1_n) + // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n) + } + if (!is_null(DfDp_out)) { + Thyra::DetachedMultiVectorView DfDp_out_view( *DfDp_out ); + DfDp_out_view(0,0) = 0.0; + DfDp_out_view(0,1) = 0.0; + DfDp_out_view(0,2) = 0.0; + DfDp_out_view(1,0) = -(f/L)*(f/L); + DfDp_out_view(1,1) = -(2.0*f/(L*L))*(a-x_in_view[0]); + DfDp_out_view(1,2) = +(2.0*f*f/(L*L*L))*(a-x_in_view[0]); + + // Compute df/dp + (df/dx_dot) * (dx_dot/dp) + (df/dx) * (dx/dp) + if (useDfDpAsTangent_ && !is_null(DxdotDp_in)) { + Thyra::ConstDetachedMultiVectorView DxdotDp( *DxdotDp_in ); + DfDp_out_view(0,0) += DxdotDp(0,0); + DfDp_out_view(0,1) += DxdotDp(0,1); + DfDp_out_view(0,2) += DxdotDp(0,2); + DfDp_out_view(1,0) += DxdotDp(1,0); + DfDp_out_view(1,1) += DxdotDp(1,1); + DfDp_out_view(1,2) += DxdotDp(1,2); + } + if (useDfDpAsTangent_ && !is_null(DxDp_in)) { + Thyra::ConstDetachedMultiVectorView DxDp( *DxDp_in ); + DfDp_out_view(0,0) += -DxDp(1,0); + DfDp_out_view(0,1) += -DxDp(1,1); + DfDp_out_view(0,2) += -DxDp(1,2); + DfDp_out_view(1,0) += (f/L)*(f/L) * DxDp(0,0); + DfDp_out_view(1,1) += (f/L)*(f/L) * DxDp(0,1); + DfDp_out_view(1,2) += (f/L)*(f/L) * DxDp(0,2); + } + } + } + + // Responses: g = x + if (acceptModelParams_) { + RCP > g_out = outArgs.get_g(0); + if (g_out != Teuchos::null) + Thyra::assign(g_out.ptr(), *x_in); + + RCP > DgDp_out = + outArgs.get_DgDp(0,0).getMultiVector(); + if (DgDp_out != Teuchos::null) + Thyra::assign(DgDp_out.ptr(), double(0.0)); + + RCP > DgDx_out = + outArgs.get_DgDx(0).getMultiVector(); + if (DgDx_out != Teuchos::null) { + Thyra::DetachedMultiVectorView DgDx_out_view( *DgDx_out ); + DgDx_out_view(0,0) = 1.0; + DgDx_out_view(0,1) = 0.0; + DgDx_out_view(1,0) = 0.0; + DgDx_out_view(1,1) = 1.0; + } + } +} + + +Teuchos::RCP > +SinCosModel:: +get_p_space(int l) const +{ + if (!acceptModelParams_) { + return Teuchos::null; + } + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ ); + if (l == 0) + return p_space_; + else if (l == 1 || l == 2) + return DxDp_space_; + return Teuchos::null; +} + +Teuchos::RCP > +SinCosModel:: +get_p_names(int l) const +{ + if (!acceptModelParams_) { + return Teuchos::null; + } + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ ); + Teuchos::RCP > p_strings = + Teuchos::rcp(new Teuchos::Array()); + if (l == 0) { + p_strings->push_back("Model Coefficient: a"); + p_strings->push_back("Model Coefficient: f"); + p_strings->push_back("Model Coefficient: L"); + } + else if (l == 1) + p_strings->push_back("DxDp"); + else if (l == 2) + p_strings->push_back("Dx_dotDp"); + return p_strings; +} + +Teuchos::RCP > +SinCosModel:: +get_g_space(int j) const +{ + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ ); + return g_space_; +} + +void +SinCosModel:: +setupInOutArgs_() const +{ + if (isInitialized_) { + return; + } + + using Teuchos::RCP; + typedef Thyra::ModelEvaluatorBase MEB; + { + // Set up prototypical InArgs + MEB::InArgsSetup inArgs; + inArgs.setModelEvalDescription(this->description()); + inArgs.setSupports( MEB::IN_ARG_t ); + inArgs.setSupports( MEB::IN_ARG_x ); + inArgs.setSupports( MEB::IN_ARG_beta ); + inArgs.setSupports( MEB::IN_ARG_x_dot ); + inArgs.setSupports( MEB::IN_ARG_alpha ); + if (acceptModelParams_) { + inArgs.set_Np(Np_); + } + inArgs_ = inArgs; + } + + { + // Set up prototypical OutArgs + MEB::OutArgsSetup outArgs; + outArgs.setModelEvalDescription(this->description()); + outArgs.setSupports( MEB::OUT_ARG_f ); + outArgs.setSupports( MEB::OUT_ARG_W_op ); + if (acceptModelParams_) { + outArgs.set_Np_Ng(Np_,Ng_); + outArgs.setSupports( MEB::OUT_ARG_DfDp,0, + MEB::DERIV_MV_JACOBIAN_FORM ); + outArgs.setSupports( MEB::OUT_ARG_DgDp,0,0, + MEB::DERIV_MV_JACOBIAN_FORM ); + outArgs.setSupports( MEB::OUT_ARG_DgDx,0, + MEB::DERIV_MV_GRADIENT_FORM ); + } + outArgs_ = outArgs; + } + // Set up nominal values + nominalValues_ = inArgs_; + if (haveIC_) + { + nominalValues_.set_t(t0_ic_); + const RCP > x_ic = createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView x_ic_view( *x_ic ); + x_ic_view[0] = a_+b_*sin((f_/L_)*t0_ic_+phi_); + x_ic_view[1] = b_*(f_/L_)*cos((f_/L_)*t0_ic_+phi_); + } + nominalValues_.set_x(x_ic); + if (acceptModelParams_) { + const RCP > p_ic = createMember(p_space_); + { + Thyra::DetachedVectorView p_ic_view( *p_ic ); + p_ic_view[0] = a_; + p_ic_view[1] = f_; + p_ic_view[2] = L_; + } + nominalValues_.set_p(0,p_ic); + } + const RCP > x_dot_ic = createMember(x_space_); + { // scope to delete DetachedVectorView + Thyra::DetachedVectorView x_dot_ic_view( *x_dot_ic ); + x_dot_ic_view[0] = b_*(f_/L_)*cos((f_/L_)*t0_ic_+phi_); + x_dot_ic_view[1] = -b_*(f_/L_)*(f_/L_)*sin((f_/L_)*t0_ic_+phi_); + } + nominalValues_.set_x_dot(x_dot_ic); + } + + isInitialized_ = true; + +} + +void +SinCosModel:: +setParameterList(Teuchos::RCP const& paramList) +{ + using Teuchos::get; + using Teuchos::ParameterList; + Teuchos::RCP tmpPL = Teuchos::rcp(new ParameterList("SinCosModel")); + if (paramList != Teuchos::null) tmpPL = paramList; + tmpPL->validateParametersAndSetDefaults(*this->getValidParameters()); + this->setMyParamList(tmpPL); + Teuchos::RCP pl = this->getMyNonconstParamList(); + bool acceptModelParams = get(*pl,"Accept model parameters"); + bool haveIC = get(*pl,"Provide nominal values"); + bool useDfDpAsTangent = get(*pl, "Use DfDp as Tangent"); + if ( (acceptModelParams != acceptModelParams_) || + (haveIC != haveIC_) + ) { + isInitialized_ = false; + } + acceptModelParams_ = acceptModelParams; + haveIC_ = haveIC; + useDfDpAsTangent_ = useDfDpAsTangent; + a_ = get(*pl,"Coeff a"); + f_ = get(*pl,"Coeff f"); + L_ = get(*pl,"Coeff L"); + x0_ic_ = get(*pl,"IC x0"); + x1_ic_ = get(*pl,"IC x1"); + t0_ic_ = get(*pl,"IC t0"); + calculateCoeffFromIC_(); + setupInOutArgs_(); +} + +Teuchos::RCP +SinCosModel:: +getValidParameters() const +{ + static Teuchos::RCP validPL; + if (is_null(validPL)) { + Teuchos::RCP pl = Teuchos::parameterList(); + pl->set("Accept model parameters", false); + pl->set("Provide nominal values", true); + pl->set("Use DfDp as Tangent", false); + pl->set("Output File Name", "Tempus_BDF2_SinCos"); + Teuchos::setDoubleParameter( + "Coeff a", 0.0, "Coefficient a in model", &*pl); + Teuchos::setDoubleParameter( + "Coeff f", 1.0, "Coefficient f in model", &*pl); + Teuchos::setDoubleParameter( + "Coeff L", 1.0, "Coefficient L in model", &*pl); + Teuchos::setDoubleParameter( + "IC x0", 0.0, "Initial Condition for x0", &*pl); + Teuchos::setDoubleParameter( + "IC x1", 1.0, "Initial Condition for x1", &*pl); + Teuchos::setDoubleParameter( + "IC t0", 0.0, "Initial time t0", &*pl); + Teuchos::setIntParameter( + "Number of Time Step Sizes", 1, "Number time step sizes for convergence study", &*pl); + validPL = pl; + } + return validPL; +} + +void +SinCosModel:: +calculateCoeffFromIC_() +{ + phi_ = atan(((f_/L_)/x1_ic_)*(x0_ic_-a_))-(f_/L_)*t0_ic_; + b_ = x1_ic_/((f_/L_)*cos((f_/L_)*t0_ic_+phi_)); +} + diff --git a/packages/piro/test/SinCosModel.hpp b/packages/piro/test/SinCosModel.hpp new file mode 100644 index 000000000000..32dfe1ee53ba --- /dev/null +++ b/packages/piro/test/SinCosModel.hpp @@ -0,0 +1,236 @@ +// @HEADER +// ************************************************************************ +// +// Piro: Strategy package for embedded analysis capabilitites +// Copyright (2010) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Andy Salinger (agsalin@sandia.gov), Sandia +// National Laboratories. +// +// ************************************************************************ +// @HEADER + +#ifndef SINCOSMODEL_H +#define SINCOSMODEL_H + +#include "Teuchos_Assert.hpp" +#include "Teuchos_RCP.hpp" +#include "Tpetra_MultiVector.hpp" +#include "Tpetra_CrsMatrix.hpp" +#include "Thyra_TpetraThyraWrappers.hpp" +#include "Thyra_ModelEvaluator.hpp" // Interface +#include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation + +#include "Teuchos_ParameterListAcceptorDefaultBase.hpp" +#include "Teuchos_ParameterList.hpp" + + +/** \brief Sine-Cosine model problem from Tempus. + * This is a canonical Sine-Cosine differential equation + * \f[ + * \mathbf{\ddot{x}}=-\mathbf{x} + * \f] + * with a few enhancements. We start with the exact solution to the + * differential equation + * \f{eqnarray*}{ + * x_{0}(t) & = & a+b*\sin((f/L)*t+\phi)\\ + * x_{1}(t) & = & b*(f/L)*\cos((f/L)*t+\phi) + * \f} + * then the form of the model is + * \f{eqnarray*}{ + * \frac{d}{dt}x_{0}(t) & = & x_{1}(t)\\ + * \frac{d}{dt}x_{1}(t) & = & \left(\frac{f}{L}\right)^{2}(a-x_{0}(t)) + * \f} + * where the default parameter values are \f$a=0\f$, \f$f=1\f$, and \f$L=1\f$, + * and the initial conditions + * \f{eqnarray*}{ + * x_{0}(t_{0}=0) & = & \gamma_{0}[=0]\\ + * x_{1}(t_{0}=0) & = & \gamma_{1}[=1] + * \f} + * determine the remaining coefficients + * \f{eqnarray*}{ + * \phi & = & \arctan(((f/L)/\gamma_{1})*(\gamma_{0}-a))-(f/L)*t_{0}[=0]\\ + * b & = & \gamma_{1}/((f/L)*cos((f/L)*t_{0}+\phi))[=1] + * \f} + + * Therefore this model has three model parameters and two initial conditions + * which effect the exact solution as above. + * \f[ + * \mathbf{p}=(a,f,L) + * \f] + * \f[ + * \dot{\mathbf{x}}=\mathbf{F}(\mathbf{x},t,\mathbf{p}) + * \f] + * where + * \f{eqnarray*}{ + * F_{0} & = & x_{1}\\ + * F_{1} & = & \left(\frac{f}{L}\right)^{2}(a-x_{0}) + * \f} + + * The exact sensitivities, \f$\mathbf{s}=\partial\mathbf{x}/\partial\mathbf{p}\f$, + * for the problem are specified as + * \f[ + * \mathbf{s}(t)=\left[\begin{array}{cc} + * 1 & 0\\ + * \left(\frac{b}{L}\right)t\,\cos\left(\left(\frac{f}{L}\right)t+\phi\right) & \left(\frac{b}{L}\right)\cos\left(\left(\frac{f}{L}\right)t+\phi\right)-\frac{b\, f\, t}{L^{2}}\sin\left(\left(\frac{f}{L}\right)t+\phi\right)\\ + * -\frac{b\, f\, t}{L^{2}}\cos\left(\left(\frac{f}{L}\right)t+\phi\right) & -\frac{b\, f}{L^{2}}\cos\left(\left(\frac{f}{L}\right)t+\phi\right)+\frac{b\, f^{2}\, t}{L^{3}}\sin\left(\left(\frac{f}{L}\right)t+\phi\right) + * \end{array}\right] + * \f] + * and for the default initial conditions, \f$\phi=0\f$ and \f$b=1\f$ + * \f[ + * \mathbf{s}(t=0)=\left[\begin{array}{cc} + * 1 & 0\\ + * 0 & \frac{b}{L}\\ + * 0 & -\frac{f}{L^{2}} + * \end{array}\right] + * \f] + * The time differentiated sensitivities, \f$\dot{\mathbf{s}}=\partial\mathbf{s}/\partial t=\partial/\partial t(\partial\mathbf{x}/\partial\mathbf{p})=\partial/\partial\mathbf{p}(\partial\mathbf{x}/\partial t)\f$ + * are + * \f[ + * \dot{\mathbf{s}}(t)=\left[\begin{array}{cc} + * 0 & 0\\ + * \left(\frac{b}{L}\right)\cos\left(\left(\frac{f}{L}\right)t+\phi\right)-\frac{b\, f\, t}{L^{2}}\sin\left(\left(\frac{f}{L}\right)t+\phi\right) & -\frac{2b\, f}{L^{2}}\sin\left(\left(\frac{f}{L}\right)t+\phi\right)\left(\frac{b}{L}\right)-\frac{b\, f^{2}\, t}{L^{3}}\cos\left(\left(\frac{f}{L}\right)t+\phi\right)\\ + * -\frac{b\, f}{L^{2}}\cos\left(\left(\frac{f}{L}\right)t+\phi\right)+\frac{b\, f^{2}\, t}{L^{3}}\sin\left(\left(\frac{f}{L}\right)t+\phi\right) & \frac{2b\, f^{2}}{L^{3}}\sin\left(\left(\frac{f}{L}\right)t+\phi\right)+\frac{b\, f^{3}\, t}{L^{4}}\cos\left(\left(\frac{f}{L}\right)t+\phi\right) + * \end{array}\right] + * \f] + */ + + +using LO = Tpetra::Map<>::local_ordinal_type; +using GO = Tpetra::Map<>::global_ordinal_type; +typedef Tpetra::Map Tpetra_Map; +typedef Tpetra::Vector Tpetra_Vector; +typedef Tpetra::MultiVector Tpetra_MultiVector; +typedef Tpetra::Operator Tpetra_Operator; +typedef Tpetra::CrsGraph Tpetra_CrsGraph; +typedef Tpetra::CrsMatrix Tpetra_CrsMatrix; +typedef Thyra::TpetraOperatorVectorExtraction ConverterT; + + + +class SinCosModel + : public Thyra::StateFuncModelEvaluatorBase, + public Teuchos::ParameterListAcceptorDefaultBase + +{ + public: + + // Constructor + SinCosModel(Teuchos::RCP pList = Teuchos::null); + + // Exact solution + Thyra::ModelEvaluatorBase::InArgs getExactSolution(double t) const; + + // Exact sensitivity solution + Thyra::ModelEvaluatorBase::InArgs getExactSensSolution(int j, double t) const; + + /** \name Public functions overridden from ModelEvaluator. */ + //@{ + + Teuchos::RCP > get_x_space() const; + Teuchos::RCP > get_f_space() const; + Thyra::ModelEvaluatorBase::InArgs getNominalValues() const; + Teuchos::RCP > create_W() const; + Teuchos::RCP > create_W_op() const; + Teuchos::RCP > get_W_factory() const; + Thyra::ModelEvaluatorBase::InArgs createInArgs() const; + + Teuchos::RCP > get_p_space(int l) const; + Teuchos::RCP > get_p_names(int l) const; + Teuchos::RCP > get_g_space(int j) const; + + //@} + + /** \name Public functions overridden from ParameterListAcceptor. */ + //@{ + void setParameterList(Teuchos::RCP const& paramList); + Teuchos::RCP getValidParameters() const; + //@} + +private: + + void setupInOutArgs_() const; + + /** \name Private functions overridden from ModelEvaluatorDefaultBase. */ + //@{ + Thyra::ModelEvaluatorBase::OutArgs createOutArgsImpl() const; + void evalModelImpl( + const Thyra::ModelEvaluatorBase::InArgs &inArgs_bar, + const Thyra::ModelEvaluatorBase::OutArgs &outArgs_bar + ) const; + //@} + + void calculateCoeffFromIC_(); + +private: + int dim_; ///< Number of state unknowns (2) + int Np_; ///< Number of parameter vectors (1) + int np_; ///< Number of parameters in this vector (2) + int Ng_; ///< Number of observation functions (1) + int ng_; ///< Number of elements in this observation function (1) + bool haveIC_; ///< false => no nominal values are provided (default=true) + bool acceptModelParams_; ///< Changes inArgs to require parameters + bool useDfDpAsTangent_; ///< Treat DfDp OutArg as tangent (df/dx*dx/dp+df/dp) + mutable bool isInitialized_; + mutable Thyra::ModelEvaluatorBase::InArgs inArgs_; + mutable Thyra::ModelEvaluatorBase::OutArgs outArgs_; + mutable Thyra::ModelEvaluatorBase::InArgs nominalValues_; + Teuchos::RCP > x_space_; + Teuchos::RCP > f_space_; + Teuchos::RCP > p_space_; + Teuchos::RCP > g_space_; + Teuchos::RCP > DxDp_space_; + + // Parameters for the model: x_0(t) = a + b*sin(f*t+phi) + // x_1(t) = b*f*cos(f*t+phi) + double a_; ///< Model parameter + double f_; ///< Model parameter + double L_; ///< Model parameter + double phi_; ///< Parameter determined from the IC + double b_; ///< Parameter determined from the IC + double t0_ic_; ///< Time value where the initial condition is specified + double x0_ic_; ///< Initial condition for x0 + double x1_ic_; ///< Initial condition for x1 +}; + + +/// Non-member constructor +//Teuchos::RCP sineCosineModel( +// Teuchos::RCP pList_) +//{ +// Teuchos::RCP model = rcp(new SinCosModel(pList_)); +// return(model); +//} + + +#endif // SINCOSMODEL diff --git a/packages/piro/test/_input_Tempus_BackwardEuler_SinCos.xml b/packages/piro/test/_input_Tempus_BackwardEuler_SinCos.xml new file mode 100644 index 000000000000..126dec31bdd4 --- /dev/null +++ b/packages/piro/test/_input_Tempus_BackwardEuler_SinCos.xml @@ -0,0 +1,110 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/packages/tempus/test/TestModels/SinCosModel_impl.hpp b/packages/tempus/test/TestModels/SinCosModel_impl.hpp index 9a0626b615cb..9a5335eb321b 100644 --- a/packages/tempus/test/TestModels/SinCosModel_impl.hpp +++ b/packages/tempus/test/TestModels/SinCosModel_impl.hpp @@ -524,11 +524,11 @@ setupInOutArgs_() const if (acceptModelParams_) { outArgs.set_Np_Ng(Np_,Ng_); outArgs.setSupports( MEB::OUT_ARG_DfDp,0, - MEB::DERIV_MV_BY_COL ); + MEB::DERIV_MV_JACOBIAN_FORM ); outArgs.setSupports( MEB::OUT_ARG_DgDp,0,0, - MEB::DERIV_MV_BY_COL ); + MEB::DERIV_MV_JACOBIAN_FORM ); outArgs.setSupports( MEB::OUT_ARG_DgDx,0, - MEB::DERIV_TRANS_MV_BY_ROW ); + MEB::DERIV_MV_GRADIENT_FORM ); } outArgs_ = outArgs; } diff --git a/packages/tempus/test/TestModels/VanDerPolModel_impl.hpp b/packages/tempus/test/TestModels/VanDerPolModel_impl.hpp index 81193259316a..7c0f4bdaeb49 100644 --- a/packages/tempus/test/TestModels/VanDerPolModel_impl.hpp +++ b/packages/tempus/test/TestModels/VanDerPolModel_impl.hpp @@ -351,7 +351,7 @@ setupInOutArgs_() const if (acceptModelParams_) { outArgs.set_Np_Ng(Np_,Ng_); outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0, - Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL ); + Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM ); } outArgs_ = outArgs; } diff --git a/packages/tempus/test/TestModels/VanDerPol_IMEXPart_ImplicitModel_impl.hpp b/packages/tempus/test/TestModels/VanDerPol_IMEXPart_ImplicitModel_impl.hpp index b9a230669ac9..487a77609225 100644 --- a/packages/tempus/test/TestModels/VanDerPol_IMEXPart_ImplicitModel_impl.hpp +++ b/packages/tempus/test/TestModels/VanDerPol_IMEXPart_ImplicitModel_impl.hpp @@ -368,9 +368,9 @@ setupInOutArgs_() const outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op ); outArgs.set_Np_Ng(Np_,Ng_); outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0, - Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL ); + Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM ); outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,4, - Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL ); + Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM ); outArgs_ = outArgs; } diff --git a/packages/tempus/test/TestModels/VanDerPol_IMEX_ExplicitModel_impl.hpp b/packages/tempus/test/TestModels/VanDerPol_IMEX_ExplicitModel_impl.hpp index 11b3e55aba59..88996bdf1af3 100644 --- a/packages/tempus/test/TestModels/VanDerPol_IMEX_ExplicitModel_impl.hpp +++ b/packages/tempus/test/TestModels/VanDerPol_IMEX_ExplicitModel_impl.hpp @@ -174,7 +174,7 @@ setupInOutArgs_() const if (acceptModelParams_) { outArgs.set_Np_Ng(Np_,Ng_); outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0, - Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL ); + Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM ); } outArgs_ = outArgs; } diff --git a/packages/tempus/test/TestModels/VanDerPol_IMEX_ImplicitModel_impl.hpp b/packages/tempus/test/TestModels/VanDerPol_IMEX_ImplicitModel_impl.hpp index 794ce1d8659b..1bdd9ccf2a97 100644 --- a/packages/tempus/test/TestModels/VanDerPol_IMEX_ImplicitModel_impl.hpp +++ b/packages/tempus/test/TestModels/VanDerPol_IMEX_ImplicitModel_impl.hpp @@ -369,7 +369,7 @@ setupInOutArgs_() const if (acceptModelParams_) { outArgs.set_Np_Ng(Np_,Ng_); outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0, - Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL ); + Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM ); } outArgs_ = outArgs; }