@acer Thank you.   I don't know how I miss that but I did.  No excuse.   

@Rouben Rostamian  No problem.  Thanks for the help.

@Carl Love Many thanks.  Below is a piece of code without alias that sends out a result.  It is not elegant, does not use the Physics package (the notation is beyond my understanding)  and the result is not in the form I want (z2=Grad(F_1).F).

restart:
#alias(x=x(t),y=y(t),f1=f1(x(t),y(t)),f2=f2(x(t),y(t))):
F:=[f1(x(t),y(t)),f2(x(t),y(t))];
phi:=[x(t),phi1(x(t),y(t))];
sys:=[diff(x(t),t)=F[1],diff(y(t),t)=F[2]];
vars:=[diff(x(t),t),diff(y(t),t)];
             F := [f1(x(t), y(t)), f2(x(t), y(t))]

                phi := [x(t), phi1(x(t), y(t))]

        [ d                          d                       ]
 sys := [--- x(t) = f1(x(t), y(t)), --- y(t) = f2(x(t), y(t))]
        [ dt                         dt                      ]

                          [ d         d      ]
                  vars := [--- x(t), --- y(t)]
                          [ dt        dt     ]

s1:=convert(op(solve(sys,vars)),list):map(print,s1):
                    d                       
                   --- x(t) = f1(x(t), y(t))
                    dt                      

                    d                       
                   --- y(t) = f2(x(t), y(t))
                    dt                      

s2aux:=map(w->diff(w,t),s1):
s2:=map(w->lhs(w)=simplify(subs(s1,rhs(w))),s2aux):map(print,s2):
          2                                            
         d                                             
        ---- x(t) = D[1](f1)(x(t), y(t)) f1(x(t), y(t))
           2                                           
         dt                                            

           + D[2](f1)(x(t), y(t)) f2(x(t), y(t))


          2                                            
         d                                             
        ---- y(t) = D[1](f2)(x(t), y(t)) f1(x(t), y(t))
           2                                           
         dt                                            

           + D[2](f2)(x(t), y(t)) f2(x(t), y(t))


z__1:=rhs(s1[1]);   #  d(Lf(x))/dt = F[1]   
z__2:=rhs(s2[1]);   #  d(Lf(Lf(x)))/dt = Grad(F[1]).F
                     z__1 := f1(x(t), y(t))

          z__2 := D[1](f1)(x(t), y(t)) f1(x(t), y(t))

             + D[2](f1)(x(t), y(t)) f2(x(t), y(t))

@Rouben Rostamian  Many thanks for the quick answer. 

I need further help to understand how the package works.  In the example below (pseudocode) I did all the calculations by going back or forward to avoid the time dependence and I wonder whether the Physics package will do the same elegantly. 

f:=[f1(x(t),y(t),f2(x(t),y(t)];
sys:={diff(x,t)=f[1],diff(y,t)=f[2]);
newsys:=subs({x=phi[1],y=phi[2]},sys);  # A transformation of Coordinates - Can we use the Physics package here?
# And another transformation of coordinates
z__1:=x;
z__2   # Lf(z__1) - Lie derivative of z__1 in f
z__3  # Lf(z__2) 
.....

Many thanks

@Carl Love and @mmcdara Many thanks.   

@Carl Love I remember reading in one of your answers only one of the time outputs of Usage matters when using Grid or Threads but I don't remember which one.  Could you remind me?

I did not have time to run a memory demanding example but a small (not that small) example showed me that Carl's solution is blazing fast.   

@Carl Love Since your solution is very concise, could you be so kind to open it up to show the details (curry, splitranks,..)?  

@Carl Love Many thanks.   Before reading your reply I tried exactly what you suggested (first option) and it worked.   

@tomleslie Many thanks.   Grid:-Set(UtilsIORelation) did the trick for removing the error message.  Now I have to check how to deal with the other issues.  

@Carl Love Many thanks.  I have just sent you an email with a link to the file with only two variables (parms10 and parms11).  

The Grid version of my procedure is running so far without any problem.   I hope Grid does not choke (as it does in other calculations for no reason).  

@Carl Love Many thanks.  Yes, "kernel panic" meant "maple: fatal error, lost connection to kernel". 

By changing the line, I got a different message - "Execution stopped: Stack limit reached."  - The limit was changed but the same msg came out.  kernelopts adjustments made no difference.  

@Carl Love Many thanks.  I am using Maple 2020.2 on both machines:  the server and the MacBook pro.  I sent the examples from the MacBook pro because I was running a much larger example on the server using your code without any modification (my attempt to modify your not-threaded code is not worth spending time). The outcome was "kernel panic". 

parse/sprintf were used to select the examples without having to modify the worksheet (I use batch files to run mpl files). As you know I had to save the intermediate results in files due to the size of the computations.  

As I said the use of the threaded version ended up with a kernel panic on the server.   

ProbSize: [[1122436, 1639162], [11224360, 17995320], [29, 29]]
 

@Carl Love Many thanks.  Here are the results for a less computationally demanding case.    SubsetPairsSparse2 is a procedure with the last line changed to [Threads:-Seq](`intersect`(seq(T[x], x= a)),a in A).

Questions:

1) Why does Maple return numcpus when kernelopts is called for the first time?

2) Why does Maple return 2 instead of 4 in some of the subsequent calls of kernelopts?  

@Carl Love The result of issuing the command shown above and running SubsetPairsSparse is

memory used=16.37GiB, alloc change=0 bytes, cpu time=5.53h, real time=25.90m, gc time=5.48h
 

There is certainly an improvement over the last calculation but did not reach 17 minutes when using Grid (as I said Grid is not a good option because I don't know when it is going to hang for no reason and memory demand).  

@vv I am not sure what you meant by all computations. No, it is no a CPU test, although if I use Grid:-Map, it shows that there is a problem when Grid is used for large computations (reported in another post to MaplePrimes).   Yes, I could use C(or Fortran, sorry but I am old school) but I am trying to have everything under one platform, that is, Maple.  I used to have parts in Mathematica, Matlab, R, and Fortran.   

@Carl Love Many thanks.  Sorry, I meant nops (old matlab habit). 

For the largest example (not discussed here), the output of Probsize is: [[1905271, 1891223], [34294878, 32150791], [30, 30]]

For the example shown in Parallelization: [[52190, 121802], [365330, 974416], [29, 29]]

I run SubsetPairsSparse with the same example and obtained the following metrics:

memory used=16.37GiB, alloc change=66.22MiB, cpu time=7.42h, real time=33.65m, gc time=7.35

approximately 7 times faster than SubsetPairs1a. If you can, please send the threaded version of it.

Since I am not sure if the procedure will work for larger problems (three trillion - a bunch of yikes), I was thinking that perhaps the following problem is feasible: given lists of sets A and B return a list L of sets in B that does not contain any set of A. 

@Carl Love @vv  @tomleslie Thank you for the procedures.  I am testing them with L1 (length=52190) and L2(lenth=121802).   Some of the procedures are taking more than 3 hours to finish.   I have modified my procedure to use Grid:-Map (Unfortunately, I don't know how to use infolevel[Threads:-Map] to check the progress, and therefore I left the use of Threads for later), and by doing that, the time is down to 17 minutes.   I guess all the procedures will need some tuning to use parallelization tools.   I can't imagine how long it will take to run the procedure when the lengths are over one million.   

@Carl Love Many, many thanks for the implementation of the bipartite graph.  I would appreciate if you could show me how to use embellishment options.  

Updates:

a) @tomleslie 's second routine

memory used=170.20GiB, alloc change=4.80MiB, cpu time=77.26h, real time=3.04h, gc time=76.90h
 

b) @Carl Love 's procedure 

memory used=0.65TiB, alloc change=0 bytes, cpu time=20.05h, real time=3.26h, gc time=17.38h
 

c) @vv 's procedure

memory used=0.74TiB, alloc change=23.93GiB, cpu time=57.76h, real time=5.09h, gc time=54.58h
 

Bipartite option returns "Error, (in GraphTheory:-DrawGraph) input graph is not bipartite" using the example shown in the post above.